Back to index

glibc  2.9
e_hypotl.c
Go to the documentation of this file.
00001 /* e_hypotl.c -- long double version of e_hypot.c.
00002  * Conversion to long double by Jakub Jelinek, jakub@redhat.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_hypotl.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
00018 #endif
00019 
00020 /* __ieee754_hypotl(x,y)
00021  *
00022  * Method :
00023  *     If (assume round-to-nearest) z=x*x+y*y
00024  *     has error less than sqrtl(2)/2 ulp, than
00025  *     sqrtl(z) has error less than 1 ulp (exercise).
00026  *
00027  *     So, compute sqrtl(x*x+y*y) with some care as
00028  *     follows to get the error below 1 ulp:
00029  *
00030  *     Assume x>y>0;
00031  *     (if possible, set rounding to round-to-nearest)
00032  *     1. if x > 2y  use
00033  *            x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
00034  *     where x1 = x with lower 64 bits cleared, x2 = x-x1; else
00035  *     2. if x <= 2y use
00036  *            t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
00037  *     where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
00038  *     y1= y with lower 64 bits chopped, y2 = y-y1.
00039  *
00040  *     NOTE: scaling may be necessary if some argument is too
00041  *           large or too tiny
00042  *
00043  * Special cases:
00044  *     hypotl(x,y) is INF if x or y is +INF or -INF; else
00045  *     hypotl(x,y) is NAN if x or y is NAN.
00046  *
00047  * Accuracy:
00048  *     hypotl(x,y) returns sqrtl(x^2+y^2) with error less
00049  *     than 1 ulps (units in the last place)
00050  */
00051 
00052 #include "math.h"
00053 #include "math_private.h"
00054 
00055 #ifdef __STDC__
00056        long double __ieee754_hypotl(long double x, long double y)
00057 #else
00058        long double __ieee754_hypotl(x,y)
00059        long double x, y;
00060 #endif
00061 {
00062        long double a,b,t1,t2,y1,y2,w;
00063        int64_t j,k,ha,hb;
00064 
00065        GET_LDOUBLE_MSW64(ha,x);
00066        ha &= 0x7fffffffffffffffLL;
00067        GET_LDOUBLE_MSW64(hb,y);
00068        hb &= 0x7fffffffffffffffLL;
00069        if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
00070        SET_LDOUBLE_MSW64(a,ha);    /* a <- |a| */
00071        SET_LDOUBLE_MSW64(b,hb);    /* b <- |b| */
00072        if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
00073        k=0;
00074        if(ha > 0x5f3f000000000000LL) {    /* a>2**8000 */
00075           if(ha >= 0x7fff000000000000LL) {       /* Inf or NaN */
00076               u_int64_t low;
00077               w = a+b;                    /* for sNaN */
00078               GET_LDOUBLE_LSW64(low,a);
00079               if(((ha&0xffffffffffffLL)|low)==0) w = a;
00080               GET_LDOUBLE_LSW64(low,b);
00081               if(((hb^0x7fff000000000000LL)|low)==0) w = b;
00082               return w;
00083           }
00084           /* scale a and b by 2**-9600 */
00085           ha -= 0x2580000000000000LL;
00086           hb -= 0x2580000000000000LL;     k += 9600;
00087           SET_LDOUBLE_MSW64(a,ha);
00088           SET_LDOUBLE_MSW64(b,hb);
00089        }
00090        if(hb < 0x20bf000000000000LL) {    /* b < 2**-8000 */
00091            if(hb <= 0x0000ffffffffffffLL) {      /* subnormal b or 0 */
00092                u_int64_t low;
00093               GET_LDOUBLE_LSW64(low,b);
00094               if((hb|low)==0) return a;
00095               t1=0;
00096               SET_LDOUBLE_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
00097               b *= t1;
00098               a *= t1;
00099               k -= 16382;
00100            } else {         /* scale a and b by 2^9600 */
00101                ha += 0x2580000000000000LL;       /* a *= 2^9600 */
00102               hb += 0x2580000000000000LL; /* b *= 2^9600 */
00103               k -= 9600;
00104               SET_LDOUBLE_MSW64(a,ha);
00105               SET_LDOUBLE_MSW64(b,hb);
00106            }
00107        }
00108     /* medium size a and b */
00109        w = a-b;
00110        if (w>b) {
00111            t1 = 0;
00112            SET_LDOUBLE_MSW64(t1,ha);
00113            t2 = a-t1;
00114            w  = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
00115        } else {
00116            a  = a+a;
00117            y1 = 0;
00118            SET_LDOUBLE_MSW64(y1,hb);
00119            y2 = b - y1;
00120            t1 = 0;
00121            SET_LDOUBLE_MSW64(t1,ha+0x0001000000000000LL);
00122            t2 = a - t1;
00123            w  = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
00124        }
00125        if(k!=0) {
00126            u_int64_t high;
00127            t1 = 1.0L;
00128            GET_LDOUBLE_MSW64(high,t1);
00129            SET_LDOUBLE_MSW64(t1,high+(k<<48));
00130            return t1*w;
00131        } else return w;
00132 }