sysdeps/ieee754/ldbl-128ibm/e_hypotl.c

   1 /* @(#)e_hypotl.c 5.1 93/09/24 */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 /* __ieee754_hypotl(x,y)
  14  *
  15  * Method :
  16  *      If (assume round-to-nearest) z=x*x+y*y
  17  *      has error less than sqrtl(2)/2 ulp, than
  18  *      sqrtl(z) has error less than 1 ulp (exercise).
  19  *
  20  *      So, compute sqrtl(x*x+y*y) with some care as
  21  *      follows to get the error below 1 ulp:
  22  *
  23  *      Assume x>y>0;
  24  *      (if possible, set rounding to round-to-nearest)
  25  *      1. if x > 2y  use
  26  *              x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
  27  *      where x1 = x with lower 53 bits cleared, x2 = x-x1; else
  28  *      2. if x <= 2y use
  29  *              t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
  30  *      where t1 = 2x with lower 53 bits cleared, t2 = 2x-t1,
  31  *      y1= y with lower 53 bits chopped, y2 = y-y1.
  32  *
  33  *      NOTE: scaling may be necessary if some argument is too
  34  *            large or too tiny
  35  *
  36  * Special cases:
  37  *      hypotl(x,y) is INF if x or y is +INF or -INF; else
  38  *      hypotl(x,y) is NAN if x or y is NAN.
  39  *
  40  * Accuracy:
  41  *      hypotl(x,y) returns sqrtl(x^2+y^2) with error less
  42  *      than 1 ulps (units in the last place)
  43  */
  44
  45 #include "math.h"
  46 #include "math_private.h"
  47
  48 static const long double two600 = 0x1.0p+600L;
  49 static const long double two1022 = 0x1.0p+1022L;
  50
  51 long double
  52 __ieee754_hypotl(long double x, long double y)
  53 {
  54         long double a,b,t1,t2,y1,y2,w,kld;
  55         int64_t j,k,ha,hb;
  56
  57         GET_LDOUBLE_MSW64(ha,x);
  58         ha &= 0x7fffffffffffffffLL;
  59         GET_LDOUBLE_MSW64(hb,y);
  60         hb &= 0x7fffffffffffffffLL;
  61         if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
  62         a = fabsl(a);   /* a <- |a| */
  63         b = fabsl(b);   /* b <- |b| */
  64         if((ha-hb)>0x3c0000000000000LL) {return a+b;} /* x/y > 2**60 */
  65         k=0;
  66         kld = 1.0L;
  67         if(ha > 0x5f30000000000000LL) { /* a>2**500 */
  68            if(ha >= 0x7ff0000000000000LL) {     /* Inf or NaN */
  69                u_int64_t low;
  70                w = a+b;                 /* for sNaN */
  71                GET_LDOUBLE_LSW64(low,a);
  72                if(((ha&0xfffffffffffffLL)|(low&0x7fffffffffffffffLL))==0)
  73                  w = a;
  74                GET_LDOUBLE_LSW64(low,b);
  75                if(((hb^0x7ff0000000000000LL)|(low&0x7fffffffffffffffLL))==0)
  76                  w = b;
  77                return w;
  78            }
  79            /* scale a and b by 2**-600 */
  80            ha -= 0x2580000000000000LL; hb -= 0x2580000000000000LL; k += 600;
  81            a /= two600;
  82            b /= two600;
  83            k += 600;
  84            kld = two600;
  85         }
  86         if(hb < 0x20b0000000000000LL) { /* b < 2**-500 */
  87             if(hb <= 0x000fffffffffffffLL) {    /* subnormal b or 0 */
  88                 u_int64_t low;
  89                 GET_LDOUBLE_LSW64(low,b);
  90                 if((hb|(low&0x7fffffffffffffffLL))==0) return a;
  91                 t1=two1022;     /* t1=2^1022 */
  92                 b *= t1;
  93                 a *= t1;
  94                 k -= 1022;
  95                 kld = kld / two1022;
  96             } else {            /* scale a and b by 2^600 */
  97                 ha += 0x2580000000000000LL;     /* a *= 2^600 */
  98                 hb += 0x2580000000000000LL;     /* b *= 2^600 */
  99                 k -= 600;
 100                 a *= two600;
 101                 b *= two600;
 102                 kld = kld / two600;
 103             }
 104         }
 105     /* medium size a and b */
 106         w = a-b;
 107         if (w>b) {
 108             SET_LDOUBLE_WORDS64(t1,ha,0);
 109             t2 = a-t1;
 110             w  = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
 111         } else {
 112             a  = a+a;
 113             SET_LDOUBLE_WORDS64(y1,hb,0);
 114             y2 = b - y1;
 115             SET_LDOUBLE_WORDS64(t1,ha+0x0010000000000000LL,0);
 116             t2 = a - t1;
 117             w  = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 118         }
 119         if(k!=0)
 120             return w*kld;
 121         else
 122             return w;
 123 }
 124 strong_alias (__ieee754_hypotl, __hypotl_finite)