2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2014 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 */
27 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */
28 /* halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c */
30 /* root.tbl uexp.tbl upow.tbl */
31 /* An ultimate power routine. Given two IEEE double machine numbers y,x */
32 /* it computes the correctly rounded (to nearest) value of x^y. */
33 /* Assumption: Machine arithmetic operations are performed in */
34 /* round to nearest mode of IEEE 754 standard. */
36 /***************************************************************************/
43 #include <math_private.h>
50 static const double huge
= 1.0e300
, tiny
= 1.0e-300;
52 double __exp1 (double x
, double xx
, double error
);
53 static double log1 (double x
, double *delta
, double *error
);
54 static double my_log2 (double x
, double *delta
, double *error
);
55 double __slowpow (double x
, double y
, double z
);
56 static double power1 (double x
, double y
);
57 static int checkint (double x
);
59 /* An ultimate power routine. Given two IEEE double machine numbers y, x it
60 computes the correctly rounded (to nearest) value of X^y. */
63 __ieee754_pow (double x
, double y
)
65 double z
, a
, aa
, error
, t
, a1
, a2
, y1
, y2
;
71 if (v
.i
[LOW_HALF
] == 0)
73 qx
= u
.i
[HIGH_HALF
] & 0x7fffffff;
75 if (((qx
== 0x7ff00000) && (u
.i
[LOW_HALF
] != 0)) || (qx
> 0x7ff00000))
87 if (((u
.i
[HIGH_HALF
] > 0 && u
.i
[HIGH_HALF
] < 0x7ff00000) || /* x>0 and not x->0 */
88 (u
.i
[HIGH_HALF
] == 0 && u
.i
[LOW_HALF
] != 0)) &&
89 /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
90 (v
.i
[HIGH_HALF
] & 0x7fffffff) < 0x4ff00000)
91 { /* if y<-1 or y>1 */
94 SET_RESTORE_ROUND (FE_TONEAREST
);
96 /* Avoid internal underflow for tiny y. The exact value of y does
97 not matter if |y| <= 2**-64. */
98 if (ABS (y
) < 0x1p
-64)
99 y
= y
< 0 ? -0x1p
-64 : 0x1p
-64;
100 z
= log1 (x
, &aa
, &error
); /* x^y =e^(y log (X)) */
108 aa
= y2
* a1
+ y
* a2
;
111 error
= error
* ABS (y
);
112 t
= __exp1 (a1
, a2
, 1.9e16
* error
); /* return -10 or 0 if wasn't computed exactly */
113 retval
= (t
> 0) ? t
: power1 (x
, y
);
120 if (((v
.i
[HIGH_HALF
] & 0x7fffffff) == 0x7ff00000 && v
.i
[LOW_HALF
] != 0)
121 || (v
.i
[HIGH_HALF
] & 0x7fffffff) > 0x7ff00000) /* NaN */
123 if (ABS (y
) > 1.0e20
)
124 return (y
> 0) ? 0 : 1.0 / 0.0;
127 return y
< 0 ? 1.0 / x
: x
;
129 return y
< 0 ? 1.0 / 0.0 : 0.0; /* return 0 */
132 qx
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
133 qy
= v
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
135 if (qx
>= 0x7ff00000 && (qx
> 0x7ff00000 || u
.i
[LOW_HALF
] != 0)) /* NaN */
137 if (qy
>= 0x7ff00000 && (qy
> 0x7ff00000 || v
.i
[LOW_HALF
] != 0)) /* NaN */
138 return x
== 1.0 ? 1.0 : y
;
141 if (u
.i
[HIGH_HALF
] < 0)
146 if (qy
== 0x7ff00000)
151 return v
.i
[HIGH_HALF
] < 0 ? INF
.x
: 0.0;
153 return v
.i
[HIGH_HALF
] < 0 ? 0.0 : INF
.x
;
155 else if (qx
== 0x7ff00000)
156 return y
< 0 ? 0.0 : INF
.x
;
157 return (x
- x
) / (x
- x
); /* y not integer and x<0 */
159 else if (qx
== 0x7ff00000)
162 return y
< 0 ? nZERO
.x
: nINF
.x
;
164 return y
< 0 ? 0.0 : INF
.x
;
166 /* if y even or odd */
167 return (k
== 1) ? __ieee754_pow (-x
, y
) : -__ieee754_pow (-x
, y
);
171 if (qx
== 0x7ff00000) /* x= 2^-0x3ff */
172 return y
> 0 ? x
: 0;
174 if (qy
> 0x45f00000 && qy
< 0x7ff00000)
179 return (x
> 1.0) ? huge
* huge
: tiny
* tiny
;
181 return (x
< 1.0) ? huge
* huge
: tiny
* tiny
;
187 return (x
> 1.0) ? INF
.x
: 0;
189 return (x
< 1.0) ? INF
.x
: 0;
190 return 0; /* unreachable, to make the compiler happy */
193 #ifndef __ieee754_pow
194 strong_alias (__ieee754_pow
, __pow_finite
)
197 /* Compute x^y using more accurate but more slow log routine. */
200 power1 (double x
, double y
)
202 double z
, a
, aa
, error
, t
, a1
, a2
, y1
, y2
;
203 z
= my_log2 (x
, &aa
, &error
);
211 aa
= ((y1
* a1
- a
) + y1
* a2
+ y2
* a1
) + y2
* a2
+ aa
* y
;
214 error
= error
* ABS (y
);
215 t
= __exp1 (a1
, a2
, 1.9e16
* error
);
216 return (t
>= 0) ? t
: __slowpow (x
, y
, z
);
219 /* Compute log(x) (x is left argument). The result is the returned double + the
220 parameter DELTA. The result is bounded by ERROR. */
223 log1 (double x
, double *delta
, double *error
)
226 double uu
, vv
, eps
, nx
, e
, e1
, e2
, t
, t1
, t2
, res
, add
= 0;
229 mynumber
/**/ two52
= {{0x43300000, 0x00000000}}; /* 2**52 */
232 mynumber
/**/ two52
= {{0x00000000, 0x43300000}}; /* 2**52 */
240 if (m
< 0x00100000) /* 1<x<2^-1007 */
248 if ((m
& 0x000fffff) < 0x0006a09e)
250 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3ff00000;
251 two52
.i
[LOW_HALF
] = (m
>> 20);
255 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3fe00000;
256 two52
.i
[LOW_HALF
] = (m
>> 20) + 1;
261 i
= (v
.i
[LOW_HALF
] & 0x000003ff) << 2;
262 if (two52
.i
[LOW_HALF
] == 1023) /* nx = 0 */
264 if (i
> 1192 && i
< 1208) /* |x-1| < 1.5*2**-10 */
267 t1
= (t
+ 5.0e6
) - 5.0e6
;
269 e1
= t
- 0.5 * t1
* t1
;
270 e2
= (t
* t
* t
* (r3
+ t
* (r4
+ t
* (r5
+ t
* (r6
+ t
272 - 0.5 * t2
* (t
+ t1
));
274 *error
= 1.0e-21 * ABS (t
);
275 *delta
= (e1
- res
) + e2
;
277 } /* |x-1| < 1.5*2**-10 */
280 v
.x
= u
.x
* (ui
.x
[i
] + ui
.x
[i
+ 1]) + bigv
.x
;
282 j
= v
.i
[LOW_HALF
] & 0x0007ffff;
286 e2
= eps
* (ui
.x
[i
+ 1] + vj
.x
[j
] * (ui
.x
[i
] + ui
.x
[i
+ 1]));
288 e2
= ((e1
- e
) + e2
);
289 t
= ui
.x
[i
+ 2] + vj
.x
[j
+ 1];
291 t2
= ((((t
- t1
) + e
) + (ui
.x
[i
+ 3] + vj
.x
[j
+ 2])) + e2
+ e
* e
292 * (p2
+ e
* (p3
+ e
* p4
)));
295 *delta
= (t1
- res
) + t2
;
302 nx
= (two52
.x
- two52e
.x
) + add
;
304 e2
= eps
* ui
.x
[i
+ 1];
307 t
= nx
* ln2a
.x
+ ui
.x
[i
+ 2];
309 t2
= ((((t
- t1
) + e
) + nx
* ln2b
.x
+ ui
.x
[i
+ 3] + e2
) + e
* e
310 * (q2
+ e
* (q3
+ e
* (q4
+ e
* (q5
+ e
* q6
)))));
313 *delta
= (t1
- res
) + t2
;
318 /* Slower but more accurate routine of log. The returned result is double +
319 DELTA. The result is bounded by ERROR. */
322 my_log2 (double x
, double *delta
, double *error
)
325 double uu
, vv
, eps
, nx
, e
, e1
, e2
, t
, t1
, t2
, res
, add
= 0;
326 double ou1
, ou2
, lu1
, lu2
, ov
, lv1
, lv2
, a
, a1
, a2
;
327 double y
, yy
, z
, zz
, j1
, j2
, j7
, j8
;
329 double j3
, j4
, j5
, j6
;
333 mynumber
/**/ two52
= {{0x43300000, 0x00000000}}; /* 2**52 */
336 mynumber
/**/ two52
= {{0x00000000, 0x43300000}}; /* 2**52 */
353 if ((m
& 0x000fffff) < 0x0006a09e)
355 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3ff00000;
356 two52
.i
[LOW_HALF
] = (m
>> 20);
360 u
.i
[HIGH_HALF
] = (m
& 0x000fffff) | 0x3fe00000;
361 two52
.i
[LOW_HALF
] = (m
>> 20) + 1;
366 i
= (v
.i
[LOW_HALF
] & 0x000003ff) << 2;
367 /*------------------------------------- |x-1| < 2**-11------------------------------- */
368 if ((two52
.i
[LOW_HALF
] == 1023) && (i
== 1200))
371 EMULV (t
, s3
, y
, yy
, j1
, j2
, j3
, j4
, j5
);
372 ADD2 (-0.5, 0, y
, yy
, z
, zz
, j1
, j2
);
373 MUL2 (t
, 0, z
, zz
, y
, yy
, j1
, j2
, j3
, j4
, j5
, j6
, j7
, j8
);
374 MUL2 (t
, 0, y
, yy
, z
, zz
, j1
, j2
, j3
, j4
, j5
, j6
, j7
, j8
);
377 e2
= ((((t
- e1
) + z
) + zz
) + t
* t
* t
378 * (ss3
+ t
* (s4
+ t
* (s5
+ t
* (s6
+ t
* (s7
+ t
* s8
))))));
380 *error
= 1.0e-25 * ABS (t
);
381 *delta
= (e1
- res
) + e2
;
384 /*----------------------------- |x-1| > 2**-11 -------------------------- */
386 { /*Computing log(x) according to log table */
387 nx
= (two52
.x
- two52e
.x
) + add
;
392 v
.x
= u
.x
* (ou1
+ ou2
) + bigv
.x
;
394 j
= v
.i
[LOW_HALF
] & 0x0007ffff;
400 a
= (ou1
+ ou2
) * (1.0 + ov
);
401 a1
= (a
+ 1.0e10
) - 1.0e10
;
402 a2
= a
* (1.0 - a1
* uu
* vv
);
407 t
= nx
* ln2a
.x
+ lu1
+ lv1
;
409 t2
= ((((t
- t1
) + e
) + (lu2
+ lv2
+ nx
* ln2b
.x
+ e2
)) + e
* e
410 * (p2
+ e
* (p3
+ e
* p4
)));
413 *delta
= (t1
- res
) + t2
;
418 /* This function receives a double x and checks if it is an integer. If not,
419 it returns 0, else it returns 1 if even or -1 if odd. */
431 m
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign */
433 return 0; /* x is +/-inf or NaN */
435 return 1; /* |x| >= 2**53 */
437 return 0; /* |x| < 2, can not be 0 or 1 */
439 k
= (m
>> 20) - 1023; /* 1 <= k <= 52 */
441 return (n
& 1) ? -1 : 1; /* odd or even */
445 return 0; /* if not integer */
446 return (n
<< (k
- 21)) ? -1 : 1;
449 return 0; /*if not integer */
451 return (m
& 1) ? -1 : 1;
454 return (m
<< (k
+ 11)) ? -1 : 1;