2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
12 /* Expansions and modifications for 128-bit long double are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, see
31 <https://www.gnu.org/licenses/>. */
33 /* __ieee754_powl(x,y) return x**y
36 * Method: Let x = 2 * (1+f)
37 * 1. Compute and return log2(x) in two pieces:
39 * where w1 has 113-53 = 60 bit trailing zeros.
40 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
41 * arithmetic, where |y'|<=0.5.
42 * 3. Return x**y = 2**n*exp(y'*log2)
45 * 1. (anything) ** 0 is 1
46 * 2. (anything) ** 1 is itself
47 * 3. (anything) ** NAN is NAN
48 * 4. NAN ** (anything except 0) is NAN
49 * 5. +-(|x| > 1) ** +INF is +INF
50 * 6. +-(|x| > 1) ** -INF is +0
51 * 7. +-(|x| < 1) ** +INF is +0
52 * 8. +-(|x| < 1) ** -INF is +INF
53 * 9. +-1 ** +-INF is NAN
54 * 10. +0 ** (+anything except 0, NAN) is +0
55 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
56 * 12. +0 ** (-anything except 0, NAN) is +INF
57 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
58 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
59 * 15. +INF ** (+anything except 0,NAN) is +INF
60 * 16. +INF ** (-anything except 0,NAN) is +0
61 * 17. -INF ** (anything) = -0 ** (-anything)
62 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
63 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
68 #include <math-barriers.h>
69 #include <math_private.h>
71 static const _Float128 bp
[] = {
77 static const _Float128 dp_h
[] = {
79 L(5.8496250072115607565592654282227158546448E-1)
82 /* Low part of log_2(1.5) */
83 static const _Float128 dp_l
[] = {
85 L(1.0579781240112554492329533686862998106046E-16)
88 static const _Float128 zero
= 0,
91 two113
= L(1.0384593717069655257060992658440192E34
),
95 /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
98 Peak relative error 2.3e-37 */
99 static const _Float128 LN
[] =
101 L(-3.0779177200290054398792536829702930623200E1
),
102 L(6.5135778082209159921251824580292116201640E1
),
103 L(-4.6312921812152436921591152809994014413540E1
),
104 L(1.2510208195629420304615674658258363295208E1
),
105 L(-9.9266909031921425609179910128531667336670E-1)
107 static const _Float128 LD
[] =
109 L(-5.129862866715009066465422805058933131960E1
),
110 L(1.452015077564081884387441590064272782044E2
),
111 L(-1.524043275549860505277434040464085593165E2
),
112 L(7.236063513651544224319663428634139768808E1
),
113 L(-1.494198912340228235853027849917095580053E1
)
117 /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
119 Peak relative error 5.7e-38 */
120 static const _Float128 PN
[] =
122 L(5.081801691915377692446852383385968225675E8
),
123 L(9.360895299872484512023336636427675327355E6
),
124 L(4.213701282274196030811629773097579432957E4
),
125 L(5.201006511142748908655720086041570288182E1
),
126 L(9.088368420359444263703202925095675982530E-3),
128 static const _Float128 PD
[] =
130 L(3.049081015149226615468111430031590411682E9
),
131 L(1.069833887183886839966085436512368982758E8
),
132 L(8.259257717868875207333991924545445705394E5
),
133 L(1.872583833284143212651746812884298360922E3
),
137 static const _Float128
139 lg2
= L(6.9314718055994530941723212145817656807550E-1),
140 lg2_h
= L(6.9314718055994528622676398299518041312695E-1),
141 lg2_l
= L(2.3190468138462996154948554638754786504121E-17),
142 ovt
= L(8.0085662595372944372e-0017),
144 cp
= L(9.6179669392597560490661645400126142495110E-1),
145 cp_h
= L(9.6179669392597555432899980587535537779331E-1),
146 cp_l
= L(5.0577616648125906047157785230014751039424E-17);
149 __ieee754_powl (_Float128 x
, _Float128 y
)
151 _Float128 z
, ax
, z_h
, z_l
, p_h
, p_l
;
152 _Float128 y1
, t1
, t2
, r
, s
, sgn
, t
, u
, v
, w
;
153 _Float128 s2
, s_h
, s_l
, t_h
, t_l
, ay
;
154 int32_t i
, j
, k
, yisint
, n
;
157 ieee854_long_double_shape_type o
, p
, q
;
161 ix
= hx
& 0x7fffffff;
165 iy
= hy
& 0x7fffffff;
168 /* y==zero: x**0 = 1 */
169 if ((iy
| q
.parts32
.w1
| q
.parts32
.w2
| q
.parts32
.w3
) == 0
173 /* 1.0**y = 1; -1.0**+-Inf = 1 */
174 if (x
== one
&& !issignaling (y
))
176 if (x
== -1 && iy
== 0x7fff0000
177 && (q
.parts32
.w1
| q
.parts32
.w2
| q
.parts32
.w3
) == 0)
180 /* +-NaN return x+y */
181 if ((ix
> 0x7fff0000)
182 || ((ix
== 0x7fff0000)
183 && ((p
.parts32
.w1
| p
.parts32
.w2
| p
.parts32
.w3
) != 0))
185 || ((iy
== 0x7fff0000)
186 && ((q
.parts32
.w1
| q
.parts32
.w2
| q
.parts32
.w3
) != 0)))
189 /* determine if y is an odd int when x < 0
190 * yisint = 0 ... y is not an integer
191 * yisint = 1 ... y is an odd int
192 * yisint = 2 ... y is an even int
197 if (iy
>= 0x40700000) /* 2^113 */
198 yisint
= 2; /* even integer y */
199 else if (iy
>= 0x3fff0000) /* 1.0 */
212 /* special value of y */
213 if ((q
.parts32
.w1
| q
.parts32
.w2
| q
.parts32
.w3
) == 0)
215 if (iy
== 0x7fff0000) /* y is +-inf */
217 if (((ix
- 0x3fff0000) | p
.parts32
.w1
| p
.parts32
.w2
| p
.parts32
.w3
)
219 return y
- y
; /* +-1**inf is NaN */
220 else if (ix
>= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */
221 return (hy
>= 0) ? y
: zero
;
222 else /* (|x|<1)**-,+inf = inf,0 */
223 return (hy
< 0) ? -y
: zero
;
225 if (iy
== 0x3fff0000)
232 if (hy
== 0x40000000)
233 return x
* x
; /* y is 2 */
234 if (hy
== 0x3ffe0000)
236 if (hx
>= 0) /* x >= +0 */
242 /* special value of x */
243 if ((p
.parts32
.w1
| p
.parts32
.w2
| p
.parts32
.w3
) == 0)
245 if (ix
== 0x7fff0000 || ix
== 0 || ix
== 0x3fff0000)
247 z
= ax
; /*x is +-0,+-inf,+-1 */
249 z
= one
/ z
; /* z = (1/|x|) */
252 if (((ix
- 0x3fff0000) | yisint
) == 0)
254 z
= (z
- z
) / (z
- z
); /* (-1)**non-int is NaN */
256 else if (yisint
== 1)
257 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
263 /* (x<0)**(non-int) is NaN */
264 if (((((uint32_t) hx
>> 31) - 1) | yisint
) == 0)
265 return (x
- x
) / (x
- x
);
267 /* sgn (sign of result -ve**odd) = -1 else = 1 */
269 if (((((uint32_t) hx
>> 31) - 1) | (yisint
- 1)) == 0)
270 sgn
= -one
; /* (-ve)**(odd int) */
273 2^-16495 = 1/2 of smallest representable value.
274 If (1 - 1/131072)^y underflows, y > 1.4986e9 */
277 /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
280 if (ix
<= 0x3ffeffff)
281 return (hy
< 0) ? huge
* huge
: tiny
* tiny
;
282 if (ix
>= 0x3fff0000)
283 return (hy
> 0) ? huge
* huge
: tiny
* tiny
;
285 /* over/underflow if x is not close to one */
287 return (hy
< 0) ? sgn
* huge
* huge
: sgn
* tiny
* tiny
;
289 return (hy
> 0) ? sgn
* huge
* huge
: sgn
* tiny
* tiny
;
294 y
= y
< 0 ? -0x1p
-128 : 0x1p
-128;
297 /* take care subnormal number */
305 n
+= ((ix
) >> 16) - 0x3fff;
307 /* determine interval */
308 ix
= j
| 0x3fff0000; /* normalize ix */
310 k
= 0; /* |x|<sqrt(3/2) */
312 k
= 1; /* |x|<sqrt(3) */
324 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
325 u
= ax
- bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
326 v
= one
/ (ax
+ bp
[k
]);
332 o
.parts32
.w2
&= 0xf8000000;
334 /* t_h=ax+bp[k] High */
338 o
.parts32
.w2
&= 0xf8000000;
340 t_l
= ax
- (t_h
- bp
[k
]);
341 s_l
= v
* ((u
- s_h
* t_h
) - s_h
* t_l
);
342 /* compute log(ax) */
344 u
= LN
[0] + s2
* (LN
[1] + s2
* (LN
[2] + s2
* (LN
[3] + s2
* LN
[4])));
345 v
= LD
[0] + s2
* (LD
[1] + s2
* (LD
[2] + s2
* (LD
[3] + s2
* (LD
[4] + s2
))));
347 r
+= s_l
* (s_h
+ s
);
352 o
.parts32
.w2
&= 0xf8000000;
354 t_l
= r
- ((t_h
- 3.0) - s2
);
355 /* u+v = s*(1+...) */
357 v
= s_l
* t_h
+ t_l
* s
;
358 /* 2/(3log2)*(s+...) */
362 o
.parts32
.w2
&= 0xf8000000;
365 z_h
= cp_h
* p_h
; /* cp_h+cp_l = 2/(3*log2) */
366 z_l
= cp_l
* p_h
+ p_l
* cp
+ dp_l
[k
];
367 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
369 t1
= (((z_h
+ z_l
) + dp_h
[k
]) + t
);
372 o
.parts32
.w2
&= 0xf8000000;
374 t2
= z_l
- (((t1
- t
) - dp_h
[k
]) - z_h
);
376 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
380 o
.parts32
.w2
&= 0xf8000000;
382 p_l
= (y
- y1
) * t1
+ y
* t2
;
387 if (j
>= 0x400d0000) /* z >= 16384 */
390 if (((j
- 0x400d0000) | o
.parts32
.w1
| o
.parts32
.w2
| o
.parts32
.w3
) != 0)
391 return sgn
* huge
* huge
; /* overflow */
394 if (p_l
+ ovt
> z
- p_h
)
395 return sgn
* huge
* huge
; /* overflow */
398 else if ((j
& 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
401 if (((j
- 0xc00d01bc) | o
.parts32
.w1
| o
.parts32
.w2
| o
.parts32
.w3
)
403 return sgn
* tiny
* tiny
; /* underflow */
407 return sgn
* tiny
* tiny
; /* underflow */
410 /* compute 2**(p_h+p_l) */
412 k
= (i
>> 16) - 0x3fff;
415 { /* if |z| > 0.5, set n = [z+0.5] */
416 n
= floorl (z
+ L(0.5));
423 o
.parts32
.w2
&= 0xf8000000;
426 v
= (p_l
- (t
- p_h
)) * lg2
+ t
* lg2_l
;
431 u
= PN
[0] + t
* (PN
[1] + t
* (PN
[2] + t
* (PN
[3] + t
* PN
[4])));
432 v
= PD
[0] + t
* (PD
[1] + t
* (PD
[2] + t
* (PD
[3] + t
)));
434 r
= (z
* t1
) / (t1
- two
) - (w
+ z
* w
);
441 z
= __scalbnl (z
, n
); /* subnormal output */
442 _Float128 force_underflow
= z
* z
;
443 math_force_eval (force_underflow
);
452 strong_alias (__ieee754_powl
, __powl_finite
)