1 /* e_powf.c -- float version of e_pow.c.
2 * 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 * ====================================================
15 * $NetBSD: e_powf.c,v 1.12 2006/03/19 20:46:25 christos Exp $
16 * $DragonFly: src/lib/libm/src/e_powf.c,v 1.2 2007/06/17 00:57:06 pavalos Exp $
20 #include "math_private.h"
22 static const float huge
= 1.0e+30, tiny
= 1.0e-30;
26 dp_h
[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
27 dp_l
[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
31 two24
= 16777216.0, /* 0x4b800000 */
32 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
33 L1
= 6.0000002384e-01, /* 0x3f19999a */
34 L2
= 4.2857143283e-01, /* 0x3edb6db7 */
35 L3
= 3.3333334327e-01, /* 0x3eaaaaab */
36 L4
= 2.7272811532e-01, /* 0x3e8ba305 */
37 L5
= 2.3066075146e-01, /* 0x3e6c3255 */
38 L6
= 2.0697501302e-01, /* 0x3e53f142 */
39 P1
= 1.6666667163e-01, /* 0x3e2aaaab */
40 P2
= -2.7777778450e-03, /* 0xbb360b61 */
41 P3
= 6.6137559770e-05, /* 0x388ab355 */
42 P4
= -1.6533901999e-06, /* 0xb5ddea0e */
43 P5
= 4.1381369442e-08, /* 0x3331bb4c */
44 lg2
= 6.9314718246e-01, /* 0x3f317218 */
45 lg2_h
= 6.93145752e-01, /* 0x3f317200 */
46 lg2_l
= 1.42860654e-06, /* 0x35bfbe8c */
47 ovt
= 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
48 cp
= 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
49 cp_h
= 9.6179199219e-01, /* 0x3f763800 =head of cp */
50 cp_l
= 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
51 ivln2
= 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
52 ivln2_h
= 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
53 ivln2_l
= 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
56 powf(float x
, float y
)
58 float z
,ax
,z_h
,z_l
,p_h
,p_l
;
59 float y1_
,t1
,t2
,r
,s
,t
,u
,v
,w
;
60 int32_t i
,j
,k
,yisint
,n
;
61 int32_t hx
,hy
,ix
,iy
,is
;
65 ix
= hx
&0x7fffffff; iy
= hy
&0x7fffffff;
67 /* y==zero: x**0 = 1 */
70 /* +-NaN return x+y */
75 /* determine if y is an odd int when x < 0
76 * yisint = 0 ... y is not an integer
77 * yisint = 1 ... y is an odd int
78 * yisint = 2 ... y is an even int
82 if(iy
>=0x4b800000) yisint
= 2; /* even integer y */
83 else if(iy
>=0x3f800000) {
84 k
= (iy
>>23)-0x7f; /* exponent */
86 if((j
<<(23-k
))==iy
) yisint
= 2-(j
&1);
90 /* special value of y */
91 if (iy
==0x7f800000) { /* y is +-inf */
93 return y
- y
; /* inf**+-1 is NaN */
94 else if (ix
> 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
95 return (hy
>=0)? y
: zero
;
96 else /* (|x|<1)**-,+inf = inf,0 */
97 return (hy
<0)?-y
: zero
;
99 if(iy
==0x3f800000) { /* y is +-1 */
100 if(hy
<0) return one
/x
; else return x
;
102 if(hy
==0x40000000) return x
*x
; /* y is 2 */
103 if(hy
==0x3f000000) { /* y is 0.5 */
104 if(hx
>=0) /* x >= +0 */
109 /* special value of x */
110 if(ix
==0x7f800000||ix
==0||ix
==0x3f800000){
111 z
= ax
; /*x is +-0,+-inf,+-1*/
112 if(hy
<0) z
= one
/z
; /* z = (1/|x|) */
114 if(((ix
-0x3f800000)|yisint
)==0) {
115 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
117 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
122 /* (x<0)**(non-int) is NaN */
123 if(((((u_int32_t
)hx
>>31)-1)|yisint
)==0) return (x
-x
)/(x
-x
);
126 if(iy
>0x4d000000) { /* if |y| > 2**27 */
127 /* over/underflow if x is not close to one */
128 if(ix
<0x3f7ffff8) return (hy
<0)? huge
*huge
:tiny
*tiny
;
129 if(ix
>0x3f800007) return (hy
>0)? huge
*huge
:tiny
*tiny
;
130 /* now |1-x| is tiny <= 2**-20, suffice to compute
131 log(x) by x-x^2/2+x^3/3-x^4/4 */
132 t
= ax
-one
; /* t has 20 trailing zeros */
133 w
= (t
*t
)*((float)0.5-t
*((float)0.333333333333-t
*(float)0.25));
134 u
= ivln2_h
*t
; /* ivln2_h has 16 sig. bits */
135 v
= t
*ivln2_l
-w
*ivln2
;
137 GET_FLOAT_WORD(is
,t1
);
138 SET_FLOAT_WORD(t1
,is
&0xfffff000);
141 float s2
,s_h
,s_l
,t_h
,t_l
;
143 /* take care subnormal number */
145 {ax
*= two24
; n
-= 24; GET_FLOAT_WORD(ix
,ax
); }
146 n
+= ((ix
)>>23)-0x7f;
148 /* determine interval */
149 ix
= j
|0x3f800000; /* normalize ix */
150 if(j
<=0x1cc471) k
=0; /* |x|<sqrt(3/2) */
151 else if(j
<0x5db3d7) k
=1; /* |x|<sqrt(3) */
152 else {k
=0;n
+=1;ix
-= 0x00800000;}
153 SET_FLOAT_WORD(ax
,ix
);
155 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
156 u
= ax
-bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
160 GET_FLOAT_WORD(is
,s_h
);
161 SET_FLOAT_WORD(s_h
,is
&0xfffff000);
162 /* t_h=ax+bp[k] High */
163 SET_FLOAT_WORD(t_h
,((ix
>>1)|0x20000000)+0x0040000+(k
<<21));
164 t_l
= ax
- (t_h
-bp
[k
]);
165 s_l
= v
*((u
-s_h
*t_h
)-s_h
*t_l
);
166 /* compute log(ax) */
168 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
171 t_h
= (float)3.0+s2
+r
;
172 GET_FLOAT_WORD(is
,t_h
);
173 SET_FLOAT_WORD(t_h
,is
&0xfffff000);
174 t_l
= r
-((t_h
-(float)3.0)-s2
);
175 /* u+v = s*(1+...) */
178 /* 2/(3log2)*(s+...) */
180 GET_FLOAT_WORD(is
,p_h
);
181 SET_FLOAT_WORD(p_h
,is
&0xfffff000);
183 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
184 z_l
= cp_l
*p_h
+p_l
*cp
+dp_l
[k
];
185 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
187 t1
= (((z_h
+z_l
)+dp_h
[k
])+t
);
188 GET_FLOAT_WORD(is
,t1
);
189 SET_FLOAT_WORD(t1
,is
&0xfffff000);
190 t2
= z_l
-(((t1
-t
)-dp_h
[k
])-z_h
);
193 s
= one
; /* s (sign of result -ve**odd) = -1 else = 1 */
194 if(((((u_int32_t
)hx
>>31)-1)|(yisint
-1))==0)
195 s
= -one
; /* (-ve)**(odd int) */
197 /* split up y into y1_+y2 and compute (y1_+y2)*(t1+t2) */
198 GET_FLOAT_WORD(is
,y
);
199 SET_FLOAT_WORD(y1_
,is
&0xfffff000);
200 p_l
= (y
-y1_
)*t1
+y
*t2
;
204 if (j
>0x43000000) /* if z > 128 */
205 return s
*huge
*huge
; /* overflow */
206 else if (j
==0x43000000) { /* if z == 128 */
207 if(p_l
+ovt
>z
-p_h
) return s
*huge
*huge
; /* overflow */
209 else if (j
==0xc3160000){ /* z == -150 */
210 if(p_l
<=z
-p_h
) return s
*tiny
*tiny
; /* underflow */
212 else if ((j
&0x7fffffff)>0x43160000) /* z <= -150 */
213 return s
*tiny
*tiny
; /* underflow */
215 * compute 2**(p_h+p_l)
220 if(i
>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
221 n
= j
+(0x00800000>>(k
+1));
222 k
= ((n
&0x7fffffff)>>23)-0x7f; /* new k for n */
223 SET_FLOAT_WORD(t
,n
&~(0x007fffff>>k
));
224 n
= ((n
&0x007fffff)|0x00800000)>>(23-k
);
229 GET_FLOAT_WORD(is
,t
);
230 SET_FLOAT_WORD(t
,is
&0xfffff000);
232 v
= (p_l
-(t
-p_h
))*lg2
+t
*lg2_l
;
236 t1
= z
- t
*(P1
+t
*(P2
+t
*(P3
+t
*(P4
+t
*P5
))));
237 r
= (z
*t1
)/(t1
-two
)-(w
+z
*w
);
241 if((j
>>23)<=0) z
= scalbnf(z
,n
); /* subnormal output */
242 else SET_FLOAT_WORD(z
,j
);