149393c00Smartynas /* @(#)e_hypot.c 5.1 93/09/24 */
249393c00Smartynas /*
349393c00Smartynas * ====================================================
449393c00Smartynas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
549393c00Smartynas *
649393c00Smartynas * Developed at SunPro, a Sun Microsystems, Inc. business.
749393c00Smartynas * Permission to use, copy, modify, and distribute this
849393c00Smartynas * software is freely granted, provided that this notice
949393c00Smartynas * is preserved.
1049393c00Smartynas * ====================================================
1149393c00Smartynas */
1249393c00Smartynas
1349393c00Smartynas /* hypotl(x,y)
1449393c00Smartynas *
1549393c00Smartynas * Method :
1649393c00Smartynas * If (assume round-to-nearest) z=x*x+y*y
1749393c00Smartynas * has error less than sqrtl(2)/2 ulp, than
1849393c00Smartynas * sqrtl(z) has error less than 1 ulp (exercise).
1949393c00Smartynas *
2049393c00Smartynas * So, compute sqrtl(x*x+y*y) with some care as
2149393c00Smartynas * follows to get the error below 1 ulp:
2249393c00Smartynas *
2349393c00Smartynas * Assume x>y>0;
2449393c00Smartynas * (if possible, set rounding to round-to-nearest)
2549393c00Smartynas * 1. if x > 2y use
2649393c00Smartynas * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
2749393c00Smartynas * where x1 = x with lower 64 bits cleared, x2 = x-x1; else
2849393c00Smartynas * 2. if x <= 2y use
2949393c00Smartynas * t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y))
3049393c00Smartynas * where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
3149393c00Smartynas * yy1= y with lower 64 bits chopped, y2 = y-yy1.
3249393c00Smartynas *
3349393c00Smartynas * NOTE: scaling may be necessary if some argument is too
3449393c00Smartynas * large or too tiny
3549393c00Smartynas *
3649393c00Smartynas * Special cases:
3749393c00Smartynas * hypotl(x,y) is INF if x or y is +INF or -INF; else
3849393c00Smartynas * hypotl(x,y) is NAN if x or y is NAN.
3949393c00Smartynas *
4049393c00Smartynas * Accuracy:
4149393c00Smartynas * hypotl(x,y) returns sqrtl(x^2+y^2) with error less
4249393c00Smartynas * than 1 ulps (units in the last place)
4349393c00Smartynas */
4449393c00Smartynas
4549393c00Smartynas #include <math.h>
4649393c00Smartynas
4749393c00Smartynas #include "math_private.h"
4849393c00Smartynas
4949393c00Smartynas long double
hypotl(long double x,long double y)5049393c00Smartynas hypotl(long double x, long double y)
5149393c00Smartynas {
5249393c00Smartynas long double a,b,t1,t2,yy1,y2,w;
5349393c00Smartynas int64_t j,k,ha,hb;
5449393c00Smartynas
5549393c00Smartynas GET_LDOUBLE_MSW64(ha,x);
5649393c00Smartynas ha &= 0x7fffffffffffffffLL;
5749393c00Smartynas GET_LDOUBLE_MSW64(hb,y);
5849393c00Smartynas hb &= 0x7fffffffffffffffLL;
5949393c00Smartynas if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
6049393c00Smartynas SET_LDOUBLE_MSW64(a,ha); /* a <- |a| */
6149393c00Smartynas SET_LDOUBLE_MSW64(b,hb); /* b <- |b| */
6249393c00Smartynas if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
6349393c00Smartynas k=0;
6449393c00Smartynas if(ha > 0x5f3f000000000000LL) { /* a>2**8000 */
6549393c00Smartynas if(ha >= 0x7fff000000000000LL) { /* Inf or NaN */
6649393c00Smartynas u_int64_t low;
6749393c00Smartynas w = a+b; /* for sNaN */
6849393c00Smartynas GET_LDOUBLE_LSW64(low,a);
6949393c00Smartynas if(((ha&0xffffffffffffLL)|low)==0) w = a;
7049393c00Smartynas GET_LDOUBLE_LSW64(low,b);
7149393c00Smartynas if(((hb^0x7fff000000000000LL)|low)==0) w = b;
7249393c00Smartynas return w;
7349393c00Smartynas }
7449393c00Smartynas /* scale a and b by 2**-9600 */
7549393c00Smartynas ha -= 0x2580000000000000LL;
7649393c00Smartynas hb -= 0x2580000000000000LL; k += 9600;
7749393c00Smartynas SET_LDOUBLE_MSW64(a,ha);
7849393c00Smartynas SET_LDOUBLE_MSW64(b,hb);
7949393c00Smartynas }
8049393c00Smartynas if(hb < 0x20bf000000000000LL) { /* b < 2**-8000 */
8149393c00Smartynas if(hb <= 0x0000ffffffffffffLL) { /* subnormal b or 0 */
8249393c00Smartynas u_int64_t low;
8349393c00Smartynas GET_LDOUBLE_LSW64(low,b);
8449393c00Smartynas if((hb|low)==0) return a;
8549393c00Smartynas t1=0;
8649393c00Smartynas SET_LDOUBLE_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
8749393c00Smartynas b *= t1;
8849393c00Smartynas a *= t1;
8949393c00Smartynas k -= 16382;
9049393c00Smartynas } else { /* scale a and b by 2^9600 */
9149393c00Smartynas ha += 0x2580000000000000LL; /* a *= 2^9600 */
9249393c00Smartynas hb += 0x2580000000000000LL; /* b *= 2^9600 */
9349393c00Smartynas k -= 9600;
9449393c00Smartynas SET_LDOUBLE_MSW64(a,ha);
9549393c00Smartynas SET_LDOUBLE_MSW64(b,hb);
9649393c00Smartynas }
9749393c00Smartynas }
9849393c00Smartynas /* medium size a and b */
9949393c00Smartynas w = a-b;
10049393c00Smartynas if (w>b) {
10149393c00Smartynas t1 = 0;
10249393c00Smartynas SET_LDOUBLE_MSW64(t1,ha);
10349393c00Smartynas t2 = a-t1;
10449393c00Smartynas w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
10549393c00Smartynas } else {
10649393c00Smartynas a = a+a;
10749393c00Smartynas yy1 = 0;
10849393c00Smartynas SET_LDOUBLE_MSW64(yy1,hb);
10949393c00Smartynas y2 = b - yy1;
11049393c00Smartynas t1 = 0;
11149393c00Smartynas SET_LDOUBLE_MSW64(t1,ha+0x0001000000000000LL);
11249393c00Smartynas t2 = a - t1;
11349393c00Smartynas w = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
11449393c00Smartynas }
11549393c00Smartynas if(k!=0) {
11649393c00Smartynas u_int64_t high;
11749393c00Smartynas t1 = 1.0L;
11849393c00Smartynas GET_LDOUBLE_MSW64(high,t1);
11949393c00Smartynas SET_LDOUBLE_MSW64(t1,high+(k<<48));
12049393c00Smartynas return t1*w;
12149393c00Smartynas } else return w;
12249393c00Smartynas }
123*2f2c0062Sguenther DEF_STD(hypotl);
124