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 | #include <math-underflow.h> |
48 | #include <libm-alias-finite.h> |
49 | |
50 | double |
51 | __ieee754_hypot (double x, double y) |
52 | { |
53 | double a, b, t1, t2, y1, y2, w; |
54 | int32_t j, k, ha, hb; |
55 | |
56 | GET_HIGH_WORD (ha, x); |
57 | ha &= 0x7fffffff; |
58 | GET_HIGH_WORD (hb, y); |
59 | hb &= 0x7fffffff; |
60 | if (hb > ha) |
61 | { |
62 | a = y; b = x; j = ha; ha = hb; hb = j; |
63 | } |
64 | else |
65 | { |
66 | a = x; b = y; |
67 | } |
68 | SET_HIGH_WORD (a, ha); /* a <- |a| */ |
69 | SET_HIGH_WORD (b, hb); /* b <- |b| */ |
70 | if ((ha - hb) > 0x3c00000) |
71 | { |
72 | return a + b; |
73 | } /* x/y > 2**60 */ |
74 | k = 0; |
75 | if (__glibc_unlikely (ha > 0x5f300000)) /* a>2**500 */ |
76 | { |
77 | if (ha >= 0x7ff00000) /* Inf or NaN */ |
78 | { |
79 | uint32_t low; |
80 | w = a + b; /* for sNaN */ |
81 | if (issignaling (a) || issignaling (b)) |
82 | return w; |
83 | GET_LOW_WORD (low, a); |
84 | if (((ha & 0xfffff) | low) == 0) |
85 | w = a; |
86 | GET_LOW_WORD (low, b); |
87 | if (((hb ^ 0x7ff00000) | low) == 0) |
88 | w = b; |
89 | return w; |
90 | } |
91 | /* scale a and b by 2**-600 */ |
92 | ha -= 0x25800000; hb -= 0x25800000; k += 600; |
93 | SET_HIGH_WORD (a, ha); |
94 | SET_HIGH_WORD (b, hb); |
95 | } |
96 | if (__builtin_expect (hb < 0x23d00000, 0)) /* b < 2**-450 */ |
97 | { |
98 | if (hb <= 0x000fffff) /* subnormal b or 0 */ |
99 | { |
100 | uint32_t low; |
101 | GET_LOW_WORD (low, b); |
102 | if ((hb | low) == 0) |
103 | return a; |
104 | t1 = 0; |
105 | SET_HIGH_WORD (t1, 0x7fd00000); /* t1=2^1022 */ |
106 | b *= t1; |
107 | a *= t1; |
108 | k -= 1022; |
109 | GET_HIGH_WORD (ha, a); |
110 | GET_HIGH_WORD (hb, b); |
111 | if (hb > ha) |
112 | { |
113 | t1 = a; |
114 | a = b; |
115 | b = t1; |
116 | j = ha; |
117 | ha = hb; |
118 | hb = j; |
119 | } |
120 | } |
121 | else /* scale a and b by 2^600 */ |
122 | { |
123 | ha += 0x25800000; /* a *= 2^600 */ |
124 | hb += 0x25800000; /* b *= 2^600 */ |
125 | k -= 600; |
126 | SET_HIGH_WORD (a, ha); |
127 | SET_HIGH_WORD (b, hb); |
128 | } |
129 | } |
130 | /* medium size a and b */ |
131 | w = a - b; |
132 | if (w > b) |
133 | { |
134 | t1 = 0; |
135 | SET_HIGH_WORD (t1, ha); |
136 | t2 = a - t1; |
137 | w = sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1))); |
138 | } |
139 | else |
140 | { |
141 | a = a + a; |
142 | y1 = 0; |
143 | SET_HIGH_WORD (y1, hb); |
144 | y2 = b - y1; |
145 | t1 = 0; |
146 | SET_HIGH_WORD (t1, ha + 0x00100000); |
147 | t2 = a - t1; |
148 | w = sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b))); |
149 | } |
150 | if (k != 0) |
151 | { |
152 | uint32_t high; |
153 | t1 = 1.0; |
154 | GET_HIGH_WORD (high, t1); |
155 | SET_HIGH_WORD (t1, high + (k << 20)); |
156 | w *= t1; |
157 | math_check_force_underflow_nonneg (w); |
158 | return w; |
159 | } |
160 | else |
161 | return w; |
162 | } |
163 | #ifndef __ieee754_hypot |
164 | libm_alias_finite (__ieee754_hypot, __hypot) |
165 | #endif |
166 | |