sysdeps/ieee754/ldbl-96/e_hypotl.c

   1 /* e_hypotl.c -- long double version of e_hypot.c.
   2  * Conversion to long double by Ulrich Drepper,
   3  * Cygnus Support, drepper@cygnus.com.
   4  */
   5
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  10  * Developed at SunPro, a Sun Microsystems, Inc. business.
  11  * Permission to use, copy, modify, and distribute this
  12  * software is freely granted, provided that this notice
  13  * is preserved.
  14  * ====================================================
  15  */
  16
  17 /* __ieee754_hypotl(x,y)
  18  *
  19  * Method :
  20  *      If (assume round-to-nearest) z=x*x+y*y
  21  *      has error less than sqrt(2)/2 ulp, than
  22  *      sqrt(z) has error less than 1 ulp (exercise).
  23  *
  24  *      So, compute sqrt(x*x+y*y) with some care as
  25  *      follows to get the error below 1 ulp:
  26  *
  27  *      Assume x>y>0;
  28  *      (if possible, set rounding to round-to-nearest)
  29  *      1. if x > 2y  use
  30  *              x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
  31  *      where x1 = x with lower 32 bits cleared, x2 = x-x1; else
  32  *      2. if x <= 2y use
  33  *              t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
  34  *      where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
  35  *      y1= y with lower 32 bits chopped, y2 = y-y1.
  36  *
  37  *      NOTE: scaling may be necessary if some argument is too
  38  *            large or too tiny
  39  *
  40  * Special cases:
  41  *      hypot(x,y) is INF if x or y is +INF or -INF; else
  42  *      hypot(x,y) is NAN if x or y is NAN.
  43  *
  44  * Accuracy:
  45  *      hypot(x,y) returns sqrt(x^2+y^2) with error less
  46  *      than 1 ulps (units in the last place)
  47  */
  48
  49 #include <math.h>
  50 #include <math_private.h>
  51
  52 long double __ieee754_hypotl(long double x, long double y)
  53 {
  54         long double a,b,t1,t2,y1,y2,w;
  55         u_int32_t j,k,ea,eb;
  56
  57         GET_LDOUBLE_EXP(ea,x);
  58         ea &= 0x7fff;
  59         GET_LDOUBLE_EXP(eb,y);
  60         eb &= 0x7fff;
  61         if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
  62         SET_LDOUBLE_EXP(a,ea);  /* a <- |a| */
  63         SET_LDOUBLE_EXP(b,eb);  /* b <- |b| */
  64         if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */
  65         k=0;
  66         if(__builtin_expect(ea > 0x5f3f,0)) {   /* a>2**8000 */
  67            if(ea == 0x7fff) {   /* Inf or NaN */
  68                u_int32_t exp __attribute__ ((unused));
  69                u_int32_t high,low;
  70                w = a+b;                 /* for sNaN */
  71                GET_LDOUBLE_WORDS(exp,high,low,a);
  72                if(((high&0x7fffffff)|low)==0) w = a;
  73                GET_LDOUBLE_WORDS(exp,high,low,b);
  74                if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
  75                return w;
  76            }
  77            /* scale a and b by 2**-9600 */
  78            ea -= 0x2580; eb -= 0x2580;  k += 9600;
  79            SET_LDOUBLE_EXP(a,ea);
  80            SET_LDOUBLE_EXP(b,eb);
  81         }
  82         if(__builtin_expect(eb < 0x20bf, 0)) {  /* b < 2**-8000 */
  83             if(eb == 0) {       /* subnormal b or 0 */
  84                 u_int32_t exp __attribute__ ((unused));
  85                 u_int32_t high,low;
  86                 GET_LDOUBLE_WORDS(exp,high,low,b);
  87                 if((high|low)==0) return a;
  88                 SET_LDOUBLE_WORDS(t1, 0x7ffd, 0x80000000, 0); /* t1=2^16382 */
  89                 b *= t1;
  90                 a *= t1;
  91                 k -= 16382;
  92                 GET_LDOUBLE_EXP (ea, a);
  93                 GET_LDOUBLE_EXP (eb, b);
  94                 if (eb > ea)
  95                   {
  96                     t1 = a;
  97                     a = b;
  98                     b = t1;
  99                     j = ea;
 100                     ea = eb;
 101                     eb = j;
 102                   }
 103             } else {            /* scale a and b by 2^9600 */
 104                 ea += 0x2580;   /* a *= 2^9600 */
 105                 eb += 0x2580;   /* b *= 2^9600 */
 106                 k -= 9600;
 107                 SET_LDOUBLE_EXP(a,ea);
 108                 SET_LDOUBLE_EXP(b,eb);
 109             }
 110         }
 111     /* medium size a and b */
 112         w = a-b;
 113         if (w>b) {
 114             u_int32_t high;
 115             GET_LDOUBLE_MSW(high,a);
 116             SET_LDOUBLE_WORDS(t1,ea,high,0);
 117             t2 = a-t1;
 118             w  = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
 119         } else {
 120             u_int32_t high;
 121             GET_LDOUBLE_MSW(high,b);
 122             a  = a+a;
 123             SET_LDOUBLE_WORDS(y1,eb,high,0);
 124             y2 = b - y1;
 125             GET_LDOUBLE_MSW(high,a);
 126             SET_LDOUBLE_WORDS(t1,ea+1,high,0);
 127             t2 = a - t1;
 128             w  = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 129         }
 130         if(k!=0) {
 131             u_int32_t exp;
 132             t1 = 1.0;
 133             GET_LDOUBLE_EXP(exp,t1);
 134             SET_LDOUBLE_EXP(t1,exp+k);
 135             return t1*w;
 136         } else return w;
 137 }
 138 strong_alias (__ieee754_hypotl, __hypotl_finite)