1 | /* @(#)e_hypot.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 | /* __ieee754_hypot(x,y) |
14 | * |
15 | * Method : |
16 | * If (assume round-to-nearest) z=x*x+y*y |
17 | * has error less than sqrt(2)/2 ulp, than |
18 | * sqrt(z) has error less than 1 ulp (exercise). |
19 | * |
20 | * So, compute sqrt(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 32 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 32 bits cleared, t2 = 2x-t1, |
31 | * y1= y with lower 32 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 | * hypot(x,y) is INF if x or y is +INF or -INF; else |
38 | * hypot(x,y) is NAN if x or y is NAN. |
39 | * |
40 | * Accuracy: |
41 | * hypot(x,y) returns sqrt(x^2+y^2) with error less |
42 | * than 1 ulps (units in the last place) |
43 | */ |
44 | |
45 | #include <math.h> |
46 | #include <math_private.h> |
47 | |
48 | double |
49 | __ieee754_hypot (double x, double y) |
50 | { |
51 | double a, b, t1, t2, y1, y2, w; |
52 | int32_t j, k, ha, hb; |
53 | |
54 | GET_HIGH_WORD (ha, x); |
55 | ha &= 0x7fffffff; |
56 | GET_HIGH_WORD (hb, y); |
57 | hb &= 0x7fffffff; |
58 | if (hb > ha) |
59 | { |
60 | a = y; b = x; j = ha; ha = hb; hb = j; |
61 | } |
62 | else |
63 | { |
64 | a = x; b = y; |
65 | } |
66 | SET_HIGH_WORD (a, ha); /* a <- |a| */ |
67 | SET_HIGH_WORD (b, hb); /* b <- |b| */ |
68 | if ((ha - hb) > 0x3c00000) |
69 | { |
70 | return a + b; |
71 | } /* x/y > 2**60 */ |
72 | k = 0; |
73 | if (__glibc_unlikely (ha > 0x5f300000)) /* a>2**500 */ |
74 | { |
75 | if (ha >= 0x7ff00000) /* Inf or NaN */ |
76 | { |
77 | uint32_t low; |
78 | w = a + b; /* for sNaN */ |
79 | if (issignaling (a) || issignaling (b)) |
80 | return w; |
81 | GET_LOW_WORD (low, a); |
82 | if (((ha & 0xfffff) | low) == 0) |
83 | w = a; |
84 | GET_LOW_WORD (low, b); |
85 | if (((hb ^ 0x7ff00000) | low) == 0) |
86 | w = b; |
87 | return w; |
88 | } |
89 | /* scale a and b by 2**-600 */ |
90 | ha -= 0x25800000; hb -= 0x25800000; k += 600; |
91 | SET_HIGH_WORD (a, ha); |
92 | SET_HIGH_WORD (b, hb); |
93 | } |
94 | if (__builtin_expect (hb < 0x23d00000, 0)) /* b < 2**-450 */ |
95 | { |
96 | if (hb <= 0x000fffff) /* subnormal b or 0 */ |
97 | { |
98 | uint32_t low; |
99 | GET_LOW_WORD (low, b); |
100 | if ((hb | low) == 0) |
101 | return a; |
102 | t1 = 0; |
103 | SET_HIGH_WORD (t1, 0x7fd00000); /* t1=2^1022 */ |
104 | b *= t1; |
105 | a *= t1; |
106 | k -= 1022; |
107 | GET_HIGH_WORD (ha, a); |
108 | GET_HIGH_WORD (hb, b); |
109 | if (hb > ha) |
110 | { |
111 | t1 = a; |
112 | a = b; |
113 | b = t1; |
114 | j = ha; |
115 | ha = hb; |
116 | hb = j; |
117 | } |
118 | } |
119 | else /* scale a and b by 2^600 */ |
120 | { |
121 | ha += 0x25800000; /* a *= 2^600 */ |
122 | hb += 0x25800000; /* b *= 2^600 */ |
123 | k -= 600; |
124 | SET_HIGH_WORD (a, ha); |
125 | SET_HIGH_WORD (b, hb); |
126 | } |
127 | } |
128 | /* medium size a and b */ |
129 | w = a - b; |
130 | if (w > b) |
131 | { |
132 | t1 = 0; |
133 | SET_HIGH_WORD (t1, ha); |
134 | t2 = a - t1; |
135 | w = __ieee754_sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1))); |
136 | } |
137 | else |
138 | { |
139 | a = a + a; |
140 | y1 = 0; |
141 | SET_HIGH_WORD (y1, hb); |
142 | y2 = b - y1; |
143 | t1 = 0; |
144 | SET_HIGH_WORD (t1, ha + 0x00100000); |
145 | t2 = a - t1; |
146 | w = __ieee754_sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b))); |
147 | } |
148 | if (k != 0) |
149 | { |
150 | uint32_t high; |
151 | t1 = 1.0; |
152 | GET_HIGH_WORD (high, t1); |
153 | SET_HIGH_WORD (t1, high + (k << 20)); |
154 | w *= t1; |
155 | math_check_force_underflow_nonneg (w); |
156 | return w; |
157 | } |
158 | else |
159 | return w; |
160 | } |
161 | strong_alias (__ieee754_hypot, __hypot_finite) |
162 | |