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