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 Ulrich Drepper,
00003  * Cygnus Support, drepper@cygnus.com.
00004  */
00005 
00006 /*
00007  * ====================================================
00008  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00009  *
00010  * Developed at SunPro, a Sun Microsystems, Inc. business.
00011  * Permission to use, copy, modify, and distribute this
00012  * software is freely granted, provided that this notice
00013  * is preserved.
00014  * ====================================================
00015  */
00016 
00017 #if defined(LIBM_SCCS) && !defined(lint)
00018 static char rcsid[] = "$NetBSD: $";
00019 #endif
00020 
00021 /* __ieee754_hypotl(x,y)
00022  *
00023  * Method :
00024  *     If (assume round-to-nearest) z=x*x+y*y
00025  *     has error less than sqrt(2)/2 ulp, than
00026  *     sqrt(z) has error less than 1 ulp (exercise).
00027  *
00028  *     So, compute sqrt(x*x+y*y) with some care as
00029  *     follows to get the error below 1 ulp:
00030  *
00031  *     Assume x>y>0;
00032  *     (if possible, set rounding to round-to-nearest)
00033  *     1. if x > 2y  use
00034  *            x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
00035  *     where x1 = x with lower 32 bits cleared, x2 = x-x1; else
00036  *     2. if x <= 2y use
00037  *            t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
00038  *     where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
00039  *     y1= y with lower 32 bits chopped, y2 = y-y1.
00040  *
00041  *     NOTE: scaling may be necessary if some argument is too
00042  *           large or too tiny
00043  *
00044  * Special cases:
00045  *     hypot(x,y) is INF if x or y is +INF or -INF; else
00046  *     hypot(x,y) is NAN if x or y is NAN.
00047  *
00048  * Accuracy:
00049  *     hypot(x,y) returns sqrt(x^2+y^2) with error less
00050  *     than 1 ulps (units in the last place)
00051  */
00052 
00053 #include "math.h"
00054 #include "math_private.h"
00055 
00056 #ifdef __STDC__
00057        long double __ieee754_hypotl(long double x, long double y)
00058 #else
00059        long double __ieee754_hypotl(x,y)
00060        long double x, y;
00061 #endif
00062 {
00063        long double a,b,t1,t2,y1,y2,w;
00064        u_int32_t j,k,ea,eb;
00065 
00066        GET_LDOUBLE_EXP(ea,x);
00067        ea &= 0x7fff;
00068        GET_LDOUBLE_EXP(eb,y);
00069        eb &= 0x7fff;
00070        if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
00071        SET_LDOUBLE_EXP(a,ea);      /* a <- |a| */
00072        SET_LDOUBLE_EXP(b,eb);      /* b <- |b| */
00073        if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */
00074        k=0;
00075        if(ea > 0x5f3f) {    /* a>2**8000 */
00076           if(ea == 0x7fff) {       /* Inf or NaN */
00077               u_int32_t exp,high,low;
00078               w = a+b;                    /* for sNaN */
00079               GET_LDOUBLE_WORDS(exp,high,low,a);
00080               if(((high&0x7fffffff)|low)==0) w = a;
00081               GET_LDOUBLE_WORDS(exp,high,low,b);
00082               if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
00083               return w;
00084           }
00085           /* scale a and b by 2**-9600 */
00086           ea -= 0x2580; eb -= 0x2580;     k += 9600;
00087           SET_LDOUBLE_EXP(a,ea);
00088           SET_LDOUBLE_EXP(b,eb);
00089        }
00090        if(eb < 0x20bf) {    /* b < 2**-8000 */
00091            if(eb == 0) {    /* subnormal b or 0 */
00092                u_int32_t exp,high,low;
00093               GET_LDOUBLE_WORDS(exp,high,low,b);
00094               if((high|low)==0) return a;
00095               SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0);      /* t1=2^16382 */
00096               b *= t1;
00097               a *= t1;
00098               k -= 16382;
00099            } else {         /* scale a and b by 2^9600 */
00100                ea += 0x2580;       /* a *= 2^9600 */
00101               eb += 0x2580; /* b *= 2^9600 */
00102               k -= 9600;
00103               SET_LDOUBLE_EXP(a,ea);
00104               SET_LDOUBLE_EXP(b,eb);
00105            }
00106        }
00107     /* medium size a and b */
00108        w = a-b;
00109        if (w>b) {
00110            u_int32_t high;
00111            GET_LDOUBLE_MSW(high,a);
00112            SET_LDOUBLE_WORDS(t1,ea,high,0);
00113            t2 = a-t1;
00114            w  = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
00115        } else {
00116            u_int32_t high;
00117            GET_LDOUBLE_MSW(high,b);
00118            a  = a+a;
00119            SET_LDOUBLE_WORDS(y1,eb,high,0);
00120            y2 = b - y1;
00121            GET_LDOUBLE_MSW(high,a);
00122            SET_LDOUBLE_WORDS(t1,ea+1,high,0);
00123            t2 = a - t1;
00124            w  = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
00125        }
00126        if(k!=0) {
00127            u_int32_t exp;
00128            t1 = 1.0;
00129            GET_LDOUBLE_EXP(exp,t1);
00130            SET_LDOUBLE_EXP(t1,exp+k);
00131            return t1*w;
00132        } else return w;
00133 }