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 #if defined(LIBM_SCCS) && !defined(lint)
  14 static char rcsid[] = "$NetBSD: e_hypotl.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
  15 #endif
  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 sqrtl(2)/2 ulp, than
  22  *      sqrtl(z) has error less than 1 ulp (exercise).
  23  *
  24  *      So, compute sqrtl(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 53 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 53 bits cleared, t2 = 2x-t1,
  35  *      y1= y with lower 53 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  *      hypotl(x,y) is INF if x or y is +INF or -INF; else
  42  *      hypotl(x,y) is NAN if x or y is NAN.
  43  *
  44  * Accuracy:
  45  *      hypotl(x,y) returns sqrtl(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 static const long double two600 = 0x1.0p+600L;
  53 static const long double two1022 = 0x1.0p+1022L;
  54
  55 #ifdef __STDC__
  56         long double __ieee754_hypotl(long double x, long double y)
  57 #else
  58         long double __ieee754_hypotl(x,y)
  59         long double x, y;
  60 #endif
  61 {
  62         long double a,b,t1,t2,y1,y2,w,kld;
  63         int64_t j,k,ha,hb;
  64
  65         GET_LDOUBLE_MSW64(ha,x);
  66         ha &= 0x7fffffffffffffffLL;
  67         GET_LDOUBLE_MSW64(hb,y);
  68         hb &= 0x7fffffffffffffffLL;
  69         if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
  70         a = fabsl(a);   /* a <- |a| */
  71         b = fabsl(b);   /* b <- |b| */
  72         if((ha-hb)>0x3c0000000000000LL) {return a+b;} /* x/y > 2**60 */
  73         k=0;
  74         kld = 1.0L;
  75         if(ha > 0x5f30000000000000LL) { /* a>2**500 */
  76            if(ha >= 0x7ff0000000000000LL) {     /* Inf or NaN */
  77                u_int64_t low;
  78                w = a+b;                 /* for sNaN */
  79                GET_LDOUBLE_LSW64(low,a);
  80                if(((ha&0xfffffffffffffLL)|(low&0x7fffffffffffffffLL))==0)
  81                  w = a;
  82                GET_LDOUBLE_LSW64(low,b);
  83                if(((hb^0x7ff0000000000000LL)|(low&0x7fffffffffffffffLL))==0)
  84                  w = b;
  85                return w;
  86            }
  87            /* scale a and b by 2**-600 */
  88            ha -= 0x2580000000000000LL; hb -= 0x2580000000000000LL; k += 600;
  89            a /= two600;
  90            b /= two600;
  91            k += 600;
  92            kld = two600;
  93         }
  94         if(hb < 0x20b0000000000000LL) { /* b < 2**-500 */
  95             if(hb <= 0x000fffffffffffffLL) {    /* subnormal b or 0 */
  96                 u_int64_t low;
  97                 GET_LDOUBLE_LSW64(low,b);
  98                 if((hb|(low&0x7fffffffffffffffLL))==0) return a;
  99                 t1=two1022;     /* t1=2^1022 */
 100                 b *= t1;
 101                 a *= t1;
 102                 k -= 1022;
 103                 kld = kld / two1022;
 104             } else {            /* scale a and b by 2^600 */
 105                 ha += 0x2580000000000000LL;     /* a *= 2^600 */
 106                 hb += 0x2580000000000000LL;     /* b *= 2^600 */
 107                 k -= 600;
 108                 a *= two600;
 109                 b *= two600;
 110                 kld = kld / two600;
 111             }
 112         }
 113     /* medium size a and b */
 114         w = a-b;
 115         if (w>b) {
 116             SET_LDOUBLE_WORDS64(t1,ha,0);
 117             t2 = a-t1;
 118             w  = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
 119         } else {
 120             a  = a+a;
 121             SET_LDOUBLE_WORDS64(y1,hb,0);
 122             y2 = b - y1;
 123             SET_LDOUBLE_WORDS64(t1,ha+0x0010000000000000LL,0);
 124             t2 = a - t1;
 125             w  = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 126         }
 127         if(k!=0)
 128             return w*kld;
 129         else
 130             return w;
 131 }