]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
688903eb | 4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
f7eac6eb | 5 | * |
e4d82761 UD |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
f7eac6eb | 10 | * |
e4d82761 UD |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
f7eac6eb | 15 | * |
e4d82761 | 16 | * You should have received a copy of the GNU Lesser General Public License |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
19 | /***************************************************************************/ |
20 | /* MODULE_NAME: upow.c */ | |
21 | /* */ | |
22 | /* FUNCTIONS: upow */ | |
e4d82761 UD |
23 | /* log1 */ |
24 | /* checkint */ | |
25 | /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */ | |
e4d82761 UD |
26 | /* root.tbl uexp.tbl upow.tbl */ |
27 | /* An ultimate power routine. Given two IEEE double machine numbers y,x */ | |
28 | /* it computes the correctly rounded (to nearest) value of x^y. */ | |
29 | /* Assumption: Machine arithmetic operations are performed in */ | |
30 | /* round to nearest mode of IEEE 754 standard. */ | |
31 | /* */ | |
32 | /***************************************************************************/ | |
4da6db51 | 33 | #include <math.h> |
e4d82761 UD |
34 | #include "endian.h" |
35 | #include "upow.h" | |
c8b3296b | 36 | #include <dla.h> |
e4d82761 UD |
37 | #include "mydefs.h" |
38 | #include "MathLib.h" | |
39 | #include "upow.tbl" | |
1ed0291c | 40 | #include <math_private.h> |
8f5b00d3 | 41 | #include <math-underflow.h> |
b7cd39e8 | 42 | #include <fenv.h> |
f7eac6eb | 43 | |
31d3cc00 UD |
44 | #ifndef SECTION |
45 | # define SECTION | |
46 | #endif | |
47 | ||
d7dd9453 | 48 | static const double huge = 1.0e300, tiny = 1.0e-300; |
f7eac6eb | 49 | |
c3d466cb WD |
50 | double __exp1 (double x, double xx); |
51 | static double log1 (double x, double *delta); | |
88576635 | 52 | static int checkint (double x); |
f7eac6eb | 53 | |
88576635 SP |
54 | /* An ultimate power routine. Given two IEEE double machine numbers y, x it |
55 | computes the correctly rounded (to nearest) value of X^y. */ | |
31d3cc00 UD |
56 | double |
57 | SECTION | |
88576635 SP |
58 | __ieee754_pow (double x, double y) |
59 | { | |
c3d466cb | 60 | double z, a, aa, t, a1, a2, y1, y2; |
88576635 | 61 | mynumber u, v; |
e4d82761 | 62 | int k; |
88576635 SP |
63 | int4 qx, qy; |
64 | v.x = y; | |
65 | u.x = x; | |
66 | if (v.i[LOW_HALF] == 0) | |
67 | { /* of y */ | |
68 | qx = u.i[HIGH_HALF] & 0x7fffffff; | |
69 | /* Is x a NaN? */ | |
72d839a4 | 70 | if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000)) |
90ab295a JM |
71 | && (y != 0 || issignaling (x))) |
72 | return x + x; | |
88576635 SP |
73 | if (y == 1.0) |
74 | return x; | |
75 | if (y == 2.0) | |
76 | return x * x; | |
77 | if (y == -1.0) | |
78 | return 1.0 / x; | |
79 | if (y == 0) | |
80 | return 1.0; | |
81 | } | |
e4d82761 | 82 | /* else */ |
88576635 SP |
83 | if (((u.i[HIGH_HALF] > 0 && u.i[HIGH_HALF] < 0x7ff00000) || /* x>0 and not x->0 */ |
84 | (u.i[HIGH_HALF] == 0 && u.i[LOW_HALF] != 0)) && | |
85 | /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */ | |
86 | (v.i[HIGH_HALF] & 0x7fffffff) < 0x4ff00000) | |
87 | { /* if y<-1 or y>1 */ | |
88 | double retval; | |
b7cd39e8 | 89 | |
4da6db51 JM |
90 | { |
91 | SET_RESTORE_ROUND (FE_TONEAREST); | |
b7cd39e8 | 92 | |
4da6db51 JM |
93 | /* Avoid internal underflow for tiny y. The exact value of y does |
94 | not matter if |y| <= 2**-64. */ | |
0e9be4db | 95 | if (fabs (y) < 0x1p-64) |
4da6db51 | 96 | y = y < 0 ? -0x1p-64 : 0x1p-64; |
c3d466cb | 97 | z = log1 (x, &aa); /* x^y =e^(y log (X)) */ |
4da6db51 JM |
98 | t = y * CN; |
99 | y1 = t - (t - y); | |
100 | y2 = y - y1; | |
101 | t = z * CN; | |
102 | a1 = t - (t - z); | |
103 | a2 = (z - a1) + aa; | |
104 | a = y1 * a1; | |
105 | aa = y2 * a1 + y * a2; | |
106 | a1 = a + aa; | |
107 | a2 = (a - a1) + aa; | |
c3d466cb WD |
108 | |
109 | /* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits). | |
110 | Maximum relative error REexp of __exp1 is 8.8e-22 (69.9 bits). | |
111 | We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp). | |
112 | Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX), | |
113 | this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp). | |
114 | So the relative error is 710 * 1.0e-21 + 8.8e-22 = 7.1e-19 | |
115 | (60.2 bits). The worst-case ULP error is 0.5064. */ | |
116 | ||
117 | retval = __exp1 (a1, a2); | |
4da6db51 | 118 | } |
b7cd39e8 | 119 | |
d81f90cc | 120 | if (isinf (retval)) |
4da6db51 JM |
121 | retval = huge * huge; |
122 | else if (retval == 0) | |
123 | retval = tiny * tiny; | |
6ace3938 JM |
124 | else |
125 | math_check_force_underflow_nonneg (retval); | |
88576635 SP |
126 | return retval; |
127 | } | |
f7eac6eb | 128 | |
88576635 SP |
129 | if (x == 0) |
130 | { | |
131 | if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) | |
132 | || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */ | |
90ab295a | 133 | return y + y; |
0e9be4db | 134 | if (fabs (y) > 1.0e20) |
88576635 SP |
135 | return (y > 0) ? 0 : 1.0 / 0.0; |
136 | k = checkint (y); | |
137 | if (k == -1) | |
138 | return y < 0 ? 1.0 / x : x; | |
139 | else | |
140 | return y < 0 ? 1.0 / 0.0 : 0.0; /* return 0 */ | |
141 | } | |
8f3edfee | 142 | |
88576635 SP |
143 | qx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */ |
144 | qy = v.i[HIGH_HALF] & 0x7fffffff; /* no sign */ | |
8f3edfee | 145 | |
88576635 | 146 | if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) /* NaN */ |
90ab295a | 147 | return x + y; |
88576635 | 148 | if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0)) /* NaN */ |
90ab295a | 149 | return x == 1.0 && !issignaling (y) ? 1.0 : y + y; |
8f3edfee | 150 | |
e4d82761 | 151 | /* if x<0 */ |
88576635 SP |
152 | if (u.i[HIGH_HALF] < 0) |
153 | { | |
154 | k = checkint (y); | |
155 | if (k == 0) | |
156 | { | |
157 | if (qy == 0x7ff00000) | |
158 | { | |
159 | if (x == -1.0) | |
160 | return 1.0; | |
161 | else if (x > -1.0) | |
162 | return v.i[HIGH_HALF] < 0 ? INF.x : 0.0; | |
163 | else | |
164 | return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x; | |
165 | } | |
166 | else if (qx == 0x7ff00000) | |
167 | return y < 0 ? 0.0 : INF.x; | |
168 | return (x - x) / (x - x); /* y not integer and x<0 */ | |
169 | } | |
8f3edfee | 170 | else if (qx == 0x7ff00000) |
88576635 SP |
171 | { |
172 | if (k < 0) | |
173 | return y < 0 ? nZERO.x : nINF.x; | |
174 | else | |
175 | return y < 0 ? 0.0 : INF.x; | |
176 | } | |
177 | /* if y even or odd */ | |
4da6db51 JM |
178 | if (k == 1) |
179 | return __ieee754_pow (-x, y); | |
180 | else | |
181 | { | |
182 | double retval; | |
183 | { | |
184 | SET_RESTORE_ROUND (FE_TONEAREST); | |
185 | retval = -__ieee754_pow (-x, y); | |
186 | } | |
d81f90cc | 187 | if (isinf (retval)) |
4da6db51 JM |
188 | retval = -huge * huge; |
189 | else if (retval == 0) | |
190 | retval = -tiny * tiny; | |
191 | return retval; | |
192 | } | |
ca58f1db | 193 | } |
e4d82761 | 194 | /* x>0 */ |
f7eac6eb | 195 | |
88576635 | 196 | if (qx == 0x7ff00000) /* x= 2^-0x3ff */ |
d91da4ce | 197 | return y > 0 ? x : 0; |
f14bd805 | 198 | |
88576635 SP |
199 | if (qy > 0x45f00000 && qy < 0x7ff00000) |
200 | { | |
201 | if (x == 1.0) | |
202 | return 1.0; | |
203 | if (y > 0) | |
204 | return (x > 1.0) ? huge * huge : tiny * tiny; | |
205 | if (y < 0) | |
206 | return (x < 1.0) ? huge * huge : tiny * tiny; | |
207 | } | |
f7eac6eb | 208 | |
88576635 SP |
209 | if (x == 1.0) |
210 | return 1.0; | |
211 | if (y > 0) | |
212 | return (x > 1.0) ? INF.x : 0; | |
213 | if (y < 0) | |
214 | return (x < 1.0) ? INF.x : 0; | |
215 | return 0; /* unreachable, to make the compiler happy */ | |
e4d82761 | 216 | } |
88576635 | 217 | |
af968f62 | 218 | #ifndef __ieee754_pow |
0ac5ae23 | 219 | strong_alias (__ieee754_pow, __pow_finite) |
af968f62 | 220 | #endif |
f7eac6eb | 221 | |
88576635 | 222 | /* Compute log(x) (x is left argument). The result is the returned double + the |
c3d466cb | 223 | parameter DELTA. */ |
31d3cc00 UD |
224 | static double |
225 | SECTION | |
c3d466cb | 226 | log1 (double x, double *delta) |
88576635 | 227 | { |
0a1f1e78 AB |
228 | unsigned int i, j; |
229 | int m; | |
88576635 SP |
230 | double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0; |
231 | mynumber u, v; | |
27692f89 | 232 | #ifdef BIG_ENDI |
88576635 | 233 | mynumber /**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */ |
27692f89 | 234 | #else |
88576635 SP |
235 | # ifdef LITTLE_ENDI |
236 | mynumber /**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */ | |
237 | # endif | |
27692f89 | 238 | #endif |
ba1ffaa1 | 239 | |
e4d82761 UD |
240 | u.x = x; |
241 | m = u.i[HIGH_HALF]; | |
c3d466cb | 242 | if (m < 0x00100000) /* Handle denormal x. */ |
88576635 SP |
243 | { |
244 | x = x * t52.x; | |
245 | add = -52.0; | |
246 | u.x = x; | |
247 | m = u.i[HIGH_HALF]; | |
248 | } | |
f7eac6eb | 249 | |
88576635 SP |
250 | if ((m & 0x000fffff) < 0x0006a09e) |
251 | { | |
252 | u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000; | |
253 | two52.i[LOW_HALF] = (m >> 20); | |
254 | } | |
e4d82761 | 255 | else |
88576635 SP |
256 | { |
257 | u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000; | |
258 | two52.i[LOW_HALF] = (m >> 20) + 1; | |
259 | } | |
f7eac6eb | 260 | |
e4d82761 UD |
261 | v.x = u.x + bigu.x; |
262 | uu = v.x - bigu.x; | |
88576635 | 263 | i = (v.i[LOW_HALF] & 0x000003ff) << 2; |
c3d466cb | 264 | if (two52.i[LOW_HALF] == 1023) /* Exponent of x is 0. */ |
88576635 SP |
265 | { |
266 | if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */ | |
267 | { | |
e4d82761 | 268 | t = x - 1.0; |
88576635 SP |
269 | t1 = (t + 5.0e6) - 5.0e6; |
270 | t2 = t - t1; | |
271 | e1 = t - 0.5 * t1 * t1; | |
272 | e2 = (t * t * t * (r3 + t * (r4 + t * (r5 + t * (r6 + t | |
273 | * (r7 + t * r8))))) | |
274 | - 0.5 * t2 * (t + t1)); | |
275 | res = e1 + e2; | |
88576635 | 276 | *delta = (e1 - res) + e2; |
c3d466cb | 277 | /* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */ |
e4d82761 | 278 | return res; |
88576635 | 279 | } /* |x-1| < 1.5*2**-10 */ |
e4d82761 | 280 | else |
88576635 SP |
281 | { |
282 | v.x = u.x * (ui.x[i] + ui.x[i + 1]) + bigv.x; | |
283 | vv = v.x - bigv.x; | |
284 | j = v.i[LOW_HALF] & 0x0007ffff; | |
285 | j = j + j + j; | |
286 | eps = u.x - uu * vv; | |
287 | e1 = eps * ui.x[i]; | |
288 | e2 = eps * (ui.x[i + 1] + vj.x[j] * (ui.x[i] + ui.x[i + 1])); | |
289 | e = e1 + e2; | |
290 | e2 = ((e1 - e) + e2); | |
291 | t = ui.x[i + 2] + vj.x[j + 1]; | |
292 | t1 = t + e; | |
293 | t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e | |
294 | * (p2 + e * (p3 + e * p4))); | |
295 | res = t1 + t2; | |
88576635 | 296 | *delta = (t1 - res) + t2; |
c3d466cb | 297 | /* Max relative error is 1.0e-24, so accurate to 79.7 bits. */ |
e4d82761 | 298 | return res; |
88576635 | 299 | } |
c3d466cb WD |
300 | } |
301 | else /* Exponent of x != 0. */ | |
88576635 | 302 | { |
e4d82761 | 303 | eps = u.x - uu; |
88576635 SP |
304 | nx = (two52.x - two52e.x) + add; |
305 | e1 = eps * ui.x[i]; | |
306 | e2 = eps * ui.x[i + 1]; | |
307 | e = e1 + e2; | |
308 | e2 = (e1 - e) + e2; | |
309 | t = nx * ln2a.x + ui.x[i + 2]; | |
310 | t1 = t + e; | |
311 | t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e | |
312 | * (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6))))); | |
313 | res = t1 + t2; | |
88576635 | 314 | *delta = (t1 - res) + t2; |
c3d466cb | 315 | /* Max relative error is 1.0e-21, so accurate to 69.7 bits. */ |
e4d82761 | 316 | return res; |
88576635 | 317 | } |
e4d82761 | 318 | } |
f7eac6eb | 319 | |
c3d466cb | 320 | |
88576635 SP |
321 | /* This function receives a double x and checks if it is an integer. If not, |
322 | it returns 0, else it returns 1 if even or -1 if odd. */ | |
31d3cc00 UD |
323 | static int |
324 | SECTION | |
88576635 SP |
325 | checkint (double x) |
326 | { | |
327 | union | |
328 | { | |
329 | int4 i[2]; | |
330 | double x; | |
331 | } u; | |
648615e1 BE |
332 | int k; |
333 | unsigned int m, n; | |
e4d82761 | 334 | u.x = x; |
88576635 SP |
335 | m = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */ |
336 | if (m >= 0x7ff00000) | |
337 | return 0; /* x is +/-inf or NaN */ | |
338 | if (m >= 0x43400000) | |
339 | return 1; /* |x| >= 2**53 */ | |
340 | if (m < 0x40000000) | |
341 | return 0; /* |x| < 2, can not be 0 or 1 */ | |
e4d82761 | 342 | n = u.i[LOW_HALF]; |
88576635 SP |
343 | k = (m >> 20) - 1023; /* 1 <= k <= 52 */ |
344 | if (k == 52) | |
345 | return (n & 1) ? -1 : 1; /* odd or even */ | |
346 | if (k > 20) | |
347 | { | |
e223d1fe | 348 | if (n << (k - 20) != 0) |
88576635 | 349 | return 0; /* if not integer */ |
e223d1fe | 350 | return (n << (k - 21) != 0) ? -1 : 1; |
88576635 SP |
351 | } |
352 | if (n) | |
353 | return 0; /*if not integer */ | |
354 | if (k == 20) | |
355 | return (m & 1) ? -1 : 1; | |
e223d1fe | 356 | if (m << (k + 12) != 0) |
88576635 | 357 | return 0; |
e223d1fe | 358 | return (m << (k + 11) != 0) ? -1 : 1; |
f7eac6eb | 359 | } |