2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 Free Software Foundation, Inc.
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.
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.
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/>.
19 /***************************************************************************/
20 /* MODULE_NAME: upow.c */
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. */
32 /***************************************************************************/
40 #include <math_private.h>
41 #include <fenv_private.h>
42 #include <math-underflow.h>
49 static const double huge
= 1.0e300
, tiny
= 1.0e-300;
51 double __exp1 (double x
, double xx
);
52 static double log1 (double x
, double *delta
);
53 static int checkint (double x
);
55 /* An ultimate power routine. Given two IEEE double machine numbers y, x it
56 computes the correctly rounded (to nearest) value of X^y. */
59 __ieee754_pow (double x
, double y
)
61 double z
, a
, aa
, t
, a1
, a2
, y1
, y2
;
67 if (v
.i
[LOW_HALF
] == 0)
69 qx
= u
.i
[HIGH_HALF
] & 0x7fffffff;
71 if ((((qx
== 0x7ff00000) && (u
.i
[LOW_HALF
] != 0)) || (qx
> 0x7ff00000))
72 && (y
!= 0 || issignaling (x
)))
84 if (((u
.i
[HIGH_HALF
] > 0 && u
.i
[HIGH_HALF
] < 0x7ff00000) || /* x>0 and not x->0 */
85 (u
.i
[HIGH_HALF
] == 0 && u
.i
[LOW_HALF
] != 0)) &&
86 /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
87 (v
.i
[HIGH_HALF
] & 0x7fffffff) < 0x4ff00000)
88 { /* if y<-1 or y>1 */
92 SET_RESTORE_ROUND (FE_TONEAREST
);
94 /* Avoid internal underflow for tiny y. The exact value of y does
95 not matter if |y| <= 2**-64. */
96 if (fabs (y
) < 0x1p
-64)
97 y
= y
< 0 ? -0x1p
-64 : 0x1p
-64;
98 z
= log1 (x
, &aa
); /* x^y =e^(y log (X)) */
106 aa
= y2
* a1
+ y
* a2
;
110 /* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits).
111 Maximum relative error REexp of __exp1 is 1.0e-18 (59.8 bits).
112 We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp).
113 Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX),
114 this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp).
115 So the relative error is 710 * 1.0e-21 + 1.0e-18 = 1.7e-18
116 (59 bits). The worst-case ULP error is 0.515. */
118 retval
= __exp1 (a1
, a2
);
122 retval
= huge
* huge
;
123 else if (retval
== 0)
124 retval
= tiny
* tiny
;
126 math_check_force_underflow_nonneg (retval
);
132 if (((v
.i
[HIGH_HALF
] & 0x7fffffff) == 0x7ff00000 && v
.i
[LOW_HALF
] != 0)
133 || (v
.i
[HIGH_HALF
] & 0x7fffffff) > 0x7ff00000) /* NaN */
135 if (fabs (y
) > 1.0e20
)
136 return (y
> 0) ? 0 : 1.0 / 0.0;
139 return y
< 0 ? 1.0 / x
: x
;
141 return y
< 0 ? 1.0 / 0.0 : 0.0; /* return 0 */
144 qx
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
145 qy
= v
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
147 if (qx
>= 0x7ff00000 && (qx
> 0x7ff00000 || u
.i
[LOW_HALF
] != 0)) /* NaN */
149 if (qy
>= 0x7ff00000 && (qy
> 0x7ff00000 || v
.i
[LOW_HALF
] != 0)) /* NaN */
150 return x
== 1.0 && !issignaling (y
) ? 1.0 : y
+ y
;
153 if (u
.i
[HIGH_HALF
] < 0)
158 if (qy
== 0x7ff00000)
163 return v
.i
[HIGH_HALF
] < 0 ? INF
.x
: 0.0;
165 return v
.i
[HIGH_HALF
] < 0 ? 0.0 : INF
.x
;
167 else if (qx
== 0x7ff00000)
168 return y
< 0 ? 0.0 : INF
.x
;
169 return (x
- x
) / (x
- x
); /* y not integer and x<0 */
171 else if (qx
== 0x7ff00000)
174 return y
< 0 ? nZERO
.x
: nINF
.x
;
176 return y
< 0 ? 0.0 : INF
.x
;
178 /* if y even or odd */
180 return __ieee754_pow (-x
, y
);
185 SET_RESTORE_ROUND (FE_TONEAREST
);
186 retval
= -__ieee754_pow (-x
, y
);
189 retval
= -huge
* huge
;
190 else if (retval
== 0)
191 retval
= -tiny
* tiny
;
197 if (qx
== 0x7ff00000) /* x= 2^-0x3ff */
198 return y
> 0 ? x
: 0;
200 if (qy
> 0x45f00000 && qy
< 0x7ff00000)
205 return (x
> 1.0) ? huge
* huge
: tiny
* tiny
;
207 return (x
< 1.0) ? huge
* huge
: tiny
* tiny
;
213 return (x
> 1.0) ? INF
.x
: 0;
215 return (x
< 1.0) ? INF
.x
: 0;
216 return 0; /* unreachable, to make the compiler happy */
219 #ifndef __ieee754_pow
220 strong_alias (__ieee754_pow
, __pow_finite
)
223 /* Compute log(x) (x is left argument). The result is the returned double + the
227 log1 (double x
, double *delta
)
231 double uu
, vv
, eps
, nx
, e
, e1
, e2
, t
, t1
, t2
, res
, add
= 0;
234 mynumber
/**/ two52
= {{0x43300000, 0x00000000}}; /* 2**52 */
237 mynumber
/**/ two52
= {{0x00000000, 0x43300000}}; /* 2**52 */
243 if (m
< 0x00100000) /* Handle denormal x. */
251 if ((m
& 0x000fffff) < 0x0006a09e)
253 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3ff00000;
254 two52
.i
[LOW_HALF
] = (m
>> 20);
258 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3fe00000;
259 two52
.i
[LOW_HALF
] = (m
>> 20) + 1;
264 i
= (v
.i
[LOW_HALF
] & 0x000003ff) << 2;
265 if (two52
.i
[LOW_HALF
] == 1023) /* Exponent of x is 0. */
267 if (i
> 1192 && i
< 1208) /* |x-1| < 1.5*2**-10 */
270 t1
= (t
+ 5.0e6
) - 5.0e6
;
272 e1
= t
- 0.5 * t1
* t1
;
273 e2
= (t
* t
* t
* (r3
+ t
* (r4
+ t
* (r5
+ t
* (r6
+ t
275 - 0.5 * t2
* (t
+ t1
));
277 *delta
= (e1
- res
) + e2
;
278 /* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */
280 } /* |x-1| < 1.5*2**-10 */
283 v
.x
= u
.x
* (ui
.x
[i
] + ui
.x
[i
+ 1]) + bigv
.x
;
285 j
= v
.i
[LOW_HALF
] & 0x0007ffff;
289 e2
= eps
* (ui
.x
[i
+ 1] + vj
.x
[j
] * (ui
.x
[i
] + ui
.x
[i
+ 1]));
291 e2
= ((e1
- e
) + e2
);
292 t
= ui
.x
[i
+ 2] + vj
.x
[j
+ 1];
294 t2
= ((((t
- t1
) + e
) + (ui
.x
[i
+ 3] + vj
.x
[j
+ 2])) + e2
+ e
* e
295 * (p2
+ e
* (p3
+ e
* p4
)));
297 *delta
= (t1
- res
) + t2
;
298 /* Max relative error is 1.0e-24, so accurate to 79.7 bits. */
302 else /* Exponent of x != 0. */
305 nx
= (two52
.x
- two52e
.x
) + add
;
307 e2
= eps
* ui
.x
[i
+ 1];
310 t
= nx
* ln2a
.x
+ ui
.x
[i
+ 2];
312 t2
= ((((t
- t1
) + e
) + nx
* ln2b
.x
+ ui
.x
[i
+ 3] + e2
) + e
* e
313 * (q2
+ e
* (q3
+ e
* (q4
+ e
* (q5
+ e
* q6
)))));
315 *delta
= (t1
- res
) + t2
;
316 /* Max relative error is 1.0e-21, so accurate to 69.7 bits. */
322 /* This function receives a double x and checks if it is an integer. If not,
323 it returns 0, else it returns 1 if even or -1 if odd. */
336 m
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
338 return 0; /* x is +/-inf or NaN */
340 return 1; /* |x| >= 2**53 */
342 return 0; /* |x| < 2, can not be 0 or 1 */
344 k
= (m
>> 20) - 1023; /* 1 <= k <= 52 */
346 return (n
& 1) ? -1 : 1; /* odd or even */
349 if (n
<< (k
- 20) != 0)
350 return 0; /* if not integer */
351 return (n
<< (k
- 21) != 0) ? -1 : 1;
354 return 0; /*if not integer */
356 return (m
& 1) ? -1 : 1;
357 if (m
<< (k
+ 12) != 0)
359 return (m
<< (k
+ 11) != 0) ? -1 : 1;