]>
Commit | Line | Data |
---|---|---|
f964490f RM |
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 | ||
f964490f RM |
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: | |
0ac5ae23 UD |
41 | * hypotl(x,y) returns sqrtl(x^2+y^2) with error less |
42 | * than 1 ulps (units in the last place) | |
f964490f RM |
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 | ||
0ac5ae23 UD |
51 | long double |
52 | __ieee754_hypotl(long double x, long double y) | |
f964490f RM |
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 */ | |
0ac5ae23 | 88 | u_int64_t low; |
f964490f RM |
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 */ | |
0ac5ae23 | 97 | ha += 0x2580000000000000LL; /* a *= 2^600 */ |
f964490f RM |
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 | } | |
0ac5ae23 | 124 | strong_alias (__ieee754_hypotl, __hypotl_finite) |