]>
git.proxmox.com Git - wasi-libc.git/blob - libc-top-half/musl/src/math/pow.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_pow.c */
3 * ====================================================
4 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
11 /* pow(x,y) return x**y
14 * Method: Let x = 2 * (1+f)
15 * 1. Compute and return log2(x) in two pieces:
17 * where w1 has 53-24 = 29 bit trailing zeros.
18 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
19 * arithmetic, where |y'|<=0.5.
20 * 3. Return x**y = 2**n*exp(y'*log2)
23 * 1. (anything) ** 0 is 1
24 * 2. 1 ** (anything) is 1
25 * 3. (anything except 1) ** NAN is NAN
26 * 4. NAN ** (anything except 0) is NAN
27 * 5. +-(|x| > 1) ** +INF is +INF
28 * 6. +-(|x| > 1) ** -INF is +0
29 * 7. +-(|x| < 1) ** +INF is +0
30 * 8. +-(|x| < 1) ** -INF is +INF
32 * 10. +0 ** (+anything except 0, NAN) is +0
33 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
34 * 12. +0 ** (-anything except 0, NAN) is +INF, raise divbyzero
35 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF, raise divbyzero
36 * 14. -0 ** (+odd integer) is -0
37 * 15. -0 ** (-odd integer) is -INF, raise divbyzero
38 * 16. +INF ** (+anything except 0,NAN) is +INF
39 * 17. +INF ** (-anything except 0,NAN) is +0
40 * 18. -INF ** (+odd integer) is -INF
41 * 19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer)
42 * 20. (anything) ** 1 is (anything)
43 * 21. (anything) ** -1 is 1/(anything)
44 * 22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
45 * 23. (-anything except 0 and inf) ** (non-integer) is NAN
48 * pow(x,y) returns x**y nearly rounded. In particular
49 * pow(integer,integer)
50 * always returns the correct integer provided it is
54 * The hexadecimal values are the intended ones for the following
55 * constants. The decimal values may be used, provided that the
56 * compiler will convert from decimal to binary accurately enough
57 * to produce the hexadecimal values shown.
64 dp_h
[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
65 dp_l
[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
66 two53
= 9007199254740992.0, /* 0x43400000, 0x00000000 */
69 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
70 L1
= 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
71 L2
= 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
72 L3
= 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
73 L4
= 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
74 L5
= 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
75 L6
= 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
76 P1
= 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
77 P2
= -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
78 P3
= 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
79 P4
= -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
80 P5
= 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
81 lg2
= 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
82 lg2_h
= 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
83 lg2_l
= -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
84 ovt
= 8.0085662595372944372e-017, /* -(1024-log2(ovfl+.5ulp)) */
85 cp
= 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
86 cp_h
= 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
87 cp_l
= -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
88 ivln2
= 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
89 ivln2_h
= 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
90 ivln2_l
= 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
92 double pow(double x
, double y
)
94 double z
,ax
,z_h
,z_l
,p_h
,p_l
;
95 double y1
,t1
,t2
,r
,s
,t
,u
,v
,w
;
96 int32_t i
,j
,k
,yisint
,n
;
100 EXTRACT_WORDS(hx
, lx
, x
);
101 EXTRACT_WORDS(hy
, ly
, y
);
102 ix
= hx
& 0x7fffffff;
103 iy
= hy
& 0x7fffffff;
105 /* x**0 = 1, even if x is NaN */
108 /* 1**y = 1, even if y is NaN */
109 if (hx
== 0x3ff00000 && lx
== 0)
111 /* NaN if either arg is NaN */
112 if (ix
> 0x7ff00000 || (ix
== 0x7ff00000 && lx
!= 0) ||
113 iy
> 0x7ff00000 || (iy
== 0x7ff00000 && ly
!= 0))
116 /* determine if y is an odd int when x < 0
117 * yisint = 0 ... y is not an integer
118 * yisint = 1 ... y is an odd int
119 * yisint = 2 ... y is an even int
123 if (iy
>= 0x43400000)
124 yisint
= 2; /* even integer y */
125 else if (iy
>= 0x3ff00000) {
126 k
= (iy
>>20) - 0x3ff; /* exponent */
128 uint32_t j
= ly
>>(52-k
);
129 if ((j
<<(52-k
)) == ly
)
131 } else if (ly
== 0) {
132 uint32_t j
= iy
>>(20-k
);
133 if ((j
<<(20-k
)) == iy
)
139 /* special value of y */
141 if (iy
== 0x7ff00000) { /* y is +-inf */
142 if (((ix
-0x3ff00000)|lx
) == 0) /* (-1)**+-inf is 1 */
144 else if (ix
>= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
145 return hy
>= 0 ? y
: 0.0;
146 else /* (|x|<1)**+-inf = 0,inf */
147 return hy
>= 0 ? 0.0 : -y
;
149 if (iy
== 0x3ff00000) { /* y is +-1 */
153 #if FLT_EVAL_METHOD!=0
155 union {double f
; uint64_t i
;} u
= {y
};
156 uint64_t i
= u
.i
& -1ULL/2;
157 if (i
>>52 == 0 && (i
&(i
-1)))
158 FORCE_EVAL((float)y
);
163 if (hy
== 0x40000000) /* y is 2 */
165 if (hy
== 0x3fe00000) { /* y is 0.5 */
166 if (hx
>= 0) /* x >= +0 */
172 /* special value of x */
174 if (ix
== 0x7ff00000 || ix
== 0 || ix
== 0x3ff00000) { /* x is +-0,+-inf,+-1 */
176 if (hy
< 0) /* z = (1/|x|) */
179 if (((ix
-0x3ff00000)|yisint
) == 0) {
180 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
181 } else if (yisint
== 1)
182 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
188 s
= 1.0; /* sign of result */
190 if (yisint
== 0) /* (x<0)**(non-int) is NaN */
192 if (yisint
== 1) /* (x<0)**(odd int) */
197 if (iy
> 0x41e00000) { /* if |y| > 2**31 */
198 if (iy
> 0x43f00000) { /* if |y| > 2**64, must o/uflow */
199 if (ix
<= 0x3fefffff)
200 return hy
< 0 ? huge
*huge
: tiny
*tiny
;
201 if (ix
>= 0x3ff00000)
202 return hy
> 0 ? huge
*huge
: tiny
*tiny
;
204 /* over/underflow if x is not close to one */
206 return hy
< 0 ? s
*huge
*huge
: s
*tiny
*tiny
;
208 return hy
> 0 ? s
*huge
*huge
: s
*tiny
*tiny
;
209 /* now |1-x| is tiny <= 2**-20, suffice to compute
210 log(x) by x-x^2/2+x^3/3-x^4/4 */
211 t
= ax
- 1.0; /* t has 20 trailing zeros */
212 w
= (t
*t
)*(0.5 - t
*(0.3333333333333333333333-t
*0.25));
213 u
= ivln2_h
*t
; /* ivln2_h has 21 sig. bits */
214 v
= t
*ivln2_l
- w
*ivln2
;
219 double ss
,s2
,s_h
,s_l
,t_h
,t_l
;
221 /* take care subnormal number */
222 if (ix
< 0x00100000) {
225 GET_HIGH_WORD(ix
,ax
);
227 n
+= ((ix
)>>20) - 0x3ff;
229 /* determine interval */
230 ix
= j
| 0x3ff00000; /* normalize ix */
231 if (j
<= 0x3988E) /* |x|<sqrt(3/2) */
233 else if (j
< 0xBB67A) /* |x|<sqrt(3) */
240 SET_HIGH_WORD(ax
, ix
);
242 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
243 u
= ax
- bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
247 SET_LOW_WORD(s_h
, 0);
248 /* t_h=ax+bp[k] High */
250 SET_HIGH_WORD(t_h
, ((ix
>>1)|0x20000000) + 0x00080000 + (k
<<18));
251 t_l
= ax
- (t_h
-bp
[k
]);
252 s_l
= v
*((u
-s_h
*t_h
)-s_h
*t_l
);
253 /* compute log(ax) */
255 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
259 SET_LOW_WORD(t_h
, 0);
260 t_l
= r
- ((t_h
-3.0)-s2
);
261 /* u+v = ss*(1+...) */
263 v
= s_l
*t_h
+ t_l
*ss
;
264 /* 2/(3log2)*(ss+...) */
266 SET_LOW_WORD(p_h
, 0);
268 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
269 z_l
= cp_l
*p_h
+p_l
*cp
+ dp_l
[k
];
270 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
272 t1
= ((z_h
+ z_l
) + dp_h
[k
]) + t
;
274 t2
= z_l
- (((t1
- t
) - dp_h
[k
]) - z_h
);
277 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
280 p_l
= (y
-y1
)*t1
+ y
*t2
;
283 EXTRACT_WORDS(j
, i
, z
);
284 if (j
>= 0x40900000) { /* z >= 1024 */
285 if (((j
-0x40900000)|i
) != 0) /* if z > 1024 */
286 return s
*huge
*huge
; /* overflow */
287 if (p_l
+ ovt
> z
- p_h
)
288 return s
*huge
*huge
; /* overflow */
289 } else if ((j
&0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */ // FIXME: instead of abs(j) use unsigned j
290 if (((j
-0xc090cc00)|i
) != 0) /* z < -1075 */
291 return s
*tiny
*tiny
; /* underflow */
293 return s
*tiny
*tiny
; /* underflow */
296 * compute 2**(p_h+p_l)
301 if (i
> 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
302 n
= j
+ (0x00100000>>(k
+1));
303 k
= ((n
&0x7fffffff)>>20) - 0x3ff; /* new k for n */
305 SET_HIGH_WORD(t
, n
& ~(0x000fffff>>k
));
306 n
= ((n
&0x000fffff)|0x00100000)>>(20-k
);
314 v
= (p_l
-(t
-p_h
))*lg2
+ t
*lg2_l
;
318 t1
= z
- t
*(P1
+t
*(P2
+t
*(P3
+t
*(P4
+t
*P5
))));
319 r
= (z
*t1
)/(t1
-2.0) - (w
+ z
*w
);
323 if ((j
>>20) <= 0) /* subnormal output */