1 /* origin: FreeBSD /usr/src/lib/msun/src/e_powf.c */
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
13 * ====================================================
20 dp_h
[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
21 dp_l
[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
22 two24
= 16777216.0, /* 0x4b800000 */
25 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
26 L1
= 6.0000002384e-01, /* 0x3f19999a */
27 L2
= 4.2857143283e-01, /* 0x3edb6db7 */
28 L3
= 3.3333334327e-01, /* 0x3eaaaaab */
29 L4
= 2.7272811532e-01, /* 0x3e8ba305 */
30 L5
= 2.3066075146e-01, /* 0x3e6c3255 */
31 L6
= 2.0697501302e-01, /* 0x3e53f142 */
32 P1
= 1.6666667163e-01, /* 0x3e2aaaab */
33 P2
= -2.7777778450e-03, /* 0xbb360b61 */
34 P3
= 6.6137559770e-05, /* 0x388ab355 */
35 P4
= -1.6533901999e-06, /* 0xb5ddea0e */
36 P5
= 4.1381369442e-08, /* 0x3331bb4c */
37 lg2
= 6.9314718246e-01, /* 0x3f317218 */
38 lg2_h
= 6.93145752e-01, /* 0x3f317200 */
39 lg2_l
= 1.42860654e-06, /* 0x35bfbe8c */
40 ovt
= 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
41 cp
= 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
42 cp_h
= 9.6191406250e-01, /* 0x3f764000 =12b cp */
43 cp_l
= -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */
44 ivln2
= 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
45 ivln2_h
= 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
46 ivln2_l
= 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
48 float powf(float x
, float y
)
50 float z
,ax
,z_h
,z_l
,p_h
,p_l
;
51 float y1
,t1
,t2
,r
,s
,sn
,t
,u
,v
,w
;
52 int32_t i
,j
,k
,yisint
,n
;
53 int32_t hx
,hy
,ix
,iy
,is
;
55 GET_FLOAT_WORD(hx
, x
);
56 GET_FLOAT_WORD(hy
, y
);
60 /* x**0 = 1, even if x is NaN */
63 /* 1**y = 1, even if y is NaN */
66 /* NaN if either arg is NaN */
67 if (ix
> 0x7f800000 || iy
> 0x7f800000)
70 /* determine if y is an odd int when x < 0
71 * yisint = 0 ... y is not an integer
72 * yisint = 1 ... y is an odd int
73 * yisint = 2 ... y is an even int
78 yisint
= 2; /* even integer y */
79 else if (iy
>= 0x3f800000) {
80 k
= (iy
>>23) - 0x7f; /* exponent */
82 if ((j
<<(23-k
)) == iy
)
87 /* special value of y */
88 if (iy
== 0x7f800000) { /* y is +-inf */
89 if (ix
== 0x3f800000) /* (-1)**+-inf is 1 */
91 else if (ix
> 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
92 return hy
>= 0 ? y
: 0.0f
;
93 else if (ix
!= 0) /* (|x|<1)**+-inf = 0,inf if x!=0 */
94 return hy
>= 0 ? 0.0f
: -y
;
96 if (iy
== 0x3f800000) /* y is +-1 */
97 return hy
>= 0 ? x
: 1.0f
/x
;
98 if (hy
== 0x40000000) /* y is 2 */
100 if (hy
== 0x3f000000) { /* y is 0.5 */
101 if (hx
>= 0) /* x >= +0 */
106 /* special value of x */
107 if (ix
== 0x7f800000 || ix
== 0 || ix
== 0x3f800000) { /* x is +-0,+-inf,+-1 */
109 if (hy
< 0) /* z = (1/|x|) */
112 if (((ix
-0x3f800000)|yisint
) == 0) {
113 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
114 } else if (yisint
== 1)
115 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
120 sn
= 1.0f
; /* sign of result */
122 if (yisint
== 0) /* (x<0)**(non-int) is NaN */
124 if (yisint
== 1) /* (x<0)**(odd int) */
129 if (iy
> 0x4d000000) { /* if |y| > 2**27 */
130 /* over/underflow if x is not close to one */
132 return hy
< 0 ? sn
*huge
*huge
: sn
*tiny
*tiny
;
134 return hy
> 0 ? sn
*huge
*huge
: sn
*tiny
*tiny
;
135 /* now |1-x| is tiny <= 2**-20, suffice to compute
136 log(x) by x-x^2/2+x^3/3-x^4/4 */
137 t
= ax
- 1; /* t has 20 trailing zeros */
138 w
= (t
*t
)*(0.5f
- t
*(0.333333333333f
- t
*0.25f
));
139 u
= ivln2_h
*t
; /* ivln2_h has 16 sig. bits */
140 v
= t
*ivln2_l
- w
*ivln2
;
142 GET_FLOAT_WORD(is
, t1
);
143 SET_FLOAT_WORD(t1
, is
& 0xfffff000);
146 float s2
,s_h
,s_l
,t_h
,t_l
;
148 /* take care subnormal number */
149 if (ix
< 0x00800000) {
152 GET_FLOAT_WORD(ix
, ax
);
154 n
+= ((ix
)>>23) - 0x7f;
156 /* determine interval */
157 ix
= j
| 0x3f800000; /* normalize ix */
158 if (j
<= 0x1cc471) /* |x|<sqrt(3/2) */
160 else if (j
< 0x5db3d7) /* |x|<sqrt(3) */
167 SET_FLOAT_WORD(ax
, ix
);
169 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
170 u
= ax
- bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
174 GET_FLOAT_WORD(is
, s_h
);
175 SET_FLOAT_WORD(s_h
, is
& 0xfffff000);
176 /* t_h=ax+bp[k] High */
177 is
= ((ix
>>1) & 0xfffff000) | 0x20000000;
178 SET_FLOAT_WORD(t_h
, is
+ 0x00400000 + (k
<<21));
179 t_l
= ax
- (t_h
- bp
[k
]);
180 s_l
= v
*((u
- s_h
*t_h
) - s_h
*t_l
);
181 /* compute log(ax) */
183 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
187 GET_FLOAT_WORD(is
, t_h
);
188 SET_FLOAT_WORD(t_h
, is
& 0xfffff000);
189 t_l
= r
- ((t_h
- 3.0f
) - s2
);
190 /* u+v = s*(1+...) */
193 /* 2/(3log2)*(s+...) */
195 GET_FLOAT_WORD(is
, p_h
);
196 SET_FLOAT_WORD(p_h
, is
& 0xfffff000);
198 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
199 z_l
= cp_l
*p_h
+ p_l
*cp
+dp_l
[k
];
200 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
202 t1
= (((z_h
+ z_l
) + dp_h
[k
]) + t
);
203 GET_FLOAT_WORD(is
, t1
);
204 SET_FLOAT_WORD(t1
, is
& 0xfffff000);
205 t2
= z_l
- (((t1
- t
) - dp_h
[k
]) - z_h
);
208 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
209 GET_FLOAT_WORD(is
, y
);
210 SET_FLOAT_WORD(y1
, is
& 0xfffff000);
211 p_l
= (y
-y1
)*t1
+ y
*t2
;
214 GET_FLOAT_WORD(j
, z
);
215 if (j
> 0x43000000) /* if z > 128 */
216 return sn
*huge
*huge
; /* overflow */
217 else if (j
== 0x43000000) { /* if z == 128 */
218 if (p_l
+ ovt
> z
- p_h
)
219 return sn
*huge
*huge
; /* overflow */
220 } else if ((j
&0x7fffffff) > 0x43160000) /* z < -150 */ // FIXME: check should be (uint32_t)j > 0xc3160000
221 return sn
*tiny
*tiny
; /* underflow */
222 else if (j
== 0xc3160000) { /* z == -150 */
224 return sn
*tiny
*tiny
; /* underflow */
227 * compute 2**(p_h+p_l)
232 if (i
> 0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
233 n
= j
+ (0x00800000>>(k
+1));
234 k
= ((n
&0x7fffffff)>>23) - 0x7f; /* new k for n */
235 SET_FLOAT_WORD(t
, n
& ~(0x007fffff>>k
));
236 n
= ((n
&0x007fffff)|0x00800000)>>(23-k
);
242 GET_FLOAT_WORD(is
, t
);
243 SET_FLOAT_WORD(t
, is
& 0xffff8000);
245 v
= (p_l
-(t
-p_h
))*lg2
+ t
*lg2_l
;
249 t1
= z
- t
*(P1
+t
*(P2
+t
*(P3
+t
*(P4
+t
*P5
))));
250 r
= (z
*t1
)/(t1
-2.0f
) - (w
+z
*w
);
252 GET_FLOAT_WORD(j
, z
);
254 if ((j
>>23) <= 0) /* subnormal output */
257 SET_FLOAT_WORD(z
, j
);