1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 Free Software Foundation, Inc.
5 *
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
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
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
18 */
19/***************************************************************************/
20/* MODULE_NAME: upow.c */
21/* */
22/* FUNCTIONS: upow */
23/* log1 */
24/* checkint */
25/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */
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/***************************************************************************/
33#include <math.h>
34#include "endian.h"
35#include "upow.h"
36#include <dla.h>
37#include "mydefs.h"
38#include "MathLib.h"
39#include "upow.tbl"
40#include <math_private.h>
41#include <math-underflow.h>
42#include <fenv.h>
43
44#ifndef SECTION
45# define SECTION
46#endif
47
48static const double huge = 1.0e300, tiny = 1.0e-300;
49
50double __exp1 (double x, double xx);
51static double log1 (double x, double *delta);
52static int checkint (double x);
53
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. */
56double
57SECTION
58__ieee754_pow (double x, double y)
59{
60 double z, a, aa, t, a1, a2, y1, y2;
61 mynumber u, v;
62 int k;
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? */
70 if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
71 && (y != 0 || issignaling (x)))
72 return x + x;
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 }
82 /* else */
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;
89
90 {
91 SET_RESTORE_ROUND (FE_TONEAREST);
92
93 /* Avoid internal underflow for tiny y. The exact value of y does
94 not matter if |y| <= 2**-64. */
95 if (fabs (y) < 0x1p-64)
96 y = y < 0 ? -0x1p-64 : 0x1p-64;
97 z = log1 (x, &aa); /* x^y =e^(y log (X)) */
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;
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);
118 }
119
120 if (isinf (retval))
121 retval = huge * huge;
122 else if (retval == 0)
123 retval = tiny * tiny;
124 else
125 math_check_force_underflow_nonneg (retval);
126 return retval;
127 }
128
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 */
133 return y + y;
134 if (fabs (y) > 1.0e20)
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 }
142
143 qx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */
144 qy = v.i[HIGH_HALF] & 0x7fffffff; /* no sign */
145
146 if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) /* NaN */
147 return x + y;
148 if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0)) /* NaN */
149 return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
150
151 /* if x<0 */
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 }
170 else if (qx == 0x7ff00000)
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 */
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 }
187 if (isinf (retval))
188 retval = -huge * huge;
189 else if (retval == 0)
190 retval = -tiny * tiny;
191 return retval;
192 }
193 }
194 /* x>0 */
195
196 if (qx == 0x7ff00000) /* x= 2^-0x3ff */
197 return y > 0 ? x : 0;
198
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 }
208
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 */
216}
217
218#ifndef __ieee754_pow
219strong_alias (__ieee754_pow, __pow_finite)
220#endif
221
222/* Compute log(x) (x is left argument). The result is the returned double + the
223 parameter DELTA. */
224static double
225SECTION
226log1 (double x, double *delta)
227{
228 unsigned int i, j;
229 int m;
230 double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
231 mynumber u, v;
232#ifdef BIG_ENDI
233 mynumber /**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */
234#else
235# ifdef LITTLE_ENDI
236 mynumber /**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */
237# endif
238#endif
239
240 u.x = x;
241 m = u.i[HIGH_HALF];
242 if (m < 0x00100000) /* Handle denormal x. */
243 {
244 x = x * t52.x;
245 add = -52.0;
246 u.x = x;
247 m = u.i[HIGH_HALF];
248 }
249
250 if ((m & 0x000fffff) < 0x0006a09e)
251 {
252 u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
253 two52.i[LOW_HALF] = (m >> 20);
254 }
255 else
256 {
257 u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
258 two52.i[LOW_HALF] = (m >> 20) + 1;
259 }
260
261 v.x = u.x + bigu.x;
262 uu = v.x - bigu.x;
263 i = (v.i[LOW_HALF] & 0x000003ff) << 2;
264 if (two52.i[LOW_HALF] == 1023) /* Exponent of x is 0. */
265 {
266 if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */
267 {
268 t = x - 1.0;
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;
276 *delta = (e1 - res) + e2;
277 /* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */
278 return res;
279 } /* |x-1| < 1.5*2**-10 */
280 else
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;
296 *delta = (t1 - res) + t2;
297 /* Max relative error is 1.0e-24, so accurate to 79.7 bits. */
298 return res;
299 }
300 }
301 else /* Exponent of x != 0. */
302 {
303 eps = u.x - uu;
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;
314 *delta = (t1 - res) + t2;
315 /* Max relative error is 1.0e-21, so accurate to 69.7 bits. */
316 return res;
317 }
318}
319
320
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. */
323static int
324SECTION
325checkint (double x)
326{
327 union
328 {
329 int4 i[2];
330 double x;
331 } u;
332 int k;
333 unsigned int m, n;
334 u.x = x;
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 */
342 n = u.i[LOW_HALF];
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 {
348 if (n << (k - 20) != 0)
349 return 0; /* if not integer */
350 return (n << (k - 21) != 0) ? -1 : 1;
351 }
352 if (n)
353 return 0; /*if not integer */
354 if (k == 20)
355 return (m & 1) ? -1 : 1;
356 if (m << (k + 12) != 0)
357 return 0;
358 return (m << (k + 11) != 0) ? -1 : 1;
359}
360