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 * ====================================================
16 #if defined(LIBM_SCCS) && !defined(lint)
17 static char rcsid
[] = "$NetBSD: e_powf.c,v 1.7 1996/04/08 15:43:44 phil Exp $";
21 #include "math_private.h"
23 static const float huge
= 1.0e+30, tiny
= 1.0e-30;
31 dp_h
[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
32 dp_l
[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
36 two24
= 16777216.0, /* 0x4b800000 */
37 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
38 L1
= 6.0000002384e-01, /* 0x3f19999a */
39 L2
= 4.2857143283e-01, /* 0x3edb6db7 */
40 L3
= 3.3333334327e-01, /* 0x3eaaaaab */
41 L4
= 2.7272811532e-01, /* 0x3e8ba305 */
42 L5
= 2.3066075146e-01, /* 0x3e6c3255 */
43 L6
= 2.0697501302e-01, /* 0x3e53f142 */
44 P1
= 1.6666667163e-01, /* 0x3e2aaaab */
45 P2
= -2.7777778450e-03, /* 0xbb360b61 */
46 P3
= 6.6137559770e-05, /* 0x388ab355 */
47 P4
= -1.6533901999e-06, /* 0xb5ddea0e */
48 P5
= 4.1381369442e-08, /* 0x3331bb4c */
49 lg2
= 6.9314718246e-01, /* 0x3f317218 */
50 lg2_h
= 6.93145752e-01, /* 0x3f317200 */
51 lg2_l
= 1.42860654e-06, /* 0x35bfbe8c */
52 ovt
= 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
53 cp
= 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
54 cp_h
= 9.6179199219e-01, /* 0x3f763800 =head of cp */
55 cp_l
= 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
56 ivln2
= 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
57 ivln2_h
= 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
58 ivln2_l
= 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
61 float __ieee754_powf(float x
, float y
)
63 float __ieee754_powf(x
,y
)
67 float z
,ax
,z_h
,z_l
,p_h
,p_l
;
68 float y1
,t1
,t2
,r
,s
,t
,u
,v
,w
;
69 int32_t i
,j
,k
,yisint
,n
;
70 int32_t hx
,hy
,ix
,iy
,is
;
74 ix
= hx
&0x7fffffff; iy
= hy
&0x7fffffff;
76 /* y==zero: x**0 = 1 */
80 if(x
== 1.0) return one
;
81 if(x
== -1.0 && isinf(y
)) return one
;
83 /* +-NaN return x+y */
88 /* determine if y is an odd int when x < 0
89 * yisint = 0 ... y is not an integer
90 * yisint = 1 ... y is an odd int
91 * yisint = 2 ... y is an even int
95 if(iy
>=0x4b800000) yisint
= 2; /* even integer y */
96 else if(iy
>=0x3f800000) {
97 k
= (iy
>>23)-0x7f; /* exponent */
99 if((j
<<(23-k
))==iy
) yisint
= 2-(j
&1);
103 /* special value of y */
104 if (iy
==0x7f800000) { /* y is +-inf */
106 return y
- y
; /* inf**+-1 is NaN */
107 else if (ix
> 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
108 return (hy
>=0)? y
: zero
;
109 else /* (|x|<1)**-,+inf = inf,0 */
110 return (hy
<0)?-y
: zero
;
112 if(iy
==0x3f800000) { /* y is +-1 */
113 if(hy
<0) return one
/x
; else return x
;
115 if(hy
==0x40000000) return x
*x
; /* y is 2 */
116 if(hy
==0x3f000000) { /* y is 0.5 */
117 if(hx
>=0) /* x >= +0 */
118 return __ieee754_sqrtf(x
);
122 /* special value of x */
123 if(ix
==0x7f800000||ix
==0||ix
==0x3f800000){
124 z
= ax
; /*x is +-0,+-inf,+-1*/
125 if(hy
<0) z
= one
/z
; /* z = (1/|x|) */
127 if(((ix
-0x3f800000)|yisint
)==0) {
128 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
130 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
135 /* (x<0)**(non-int) is NaN */
136 if(((((u_int32_t
)hx
>>31)-1)|yisint
)==0) return (x
-x
)/(x
-x
);
139 if(iy
>0x4d000000) { /* if |y| > 2**27 */
140 /* over/underflow if x is not close to one */
141 if(ix
<0x3f7ffff8) return (hy
<0)? huge
*huge
:tiny
*tiny
;
142 if(ix
>0x3f800007) return (hy
>0)? huge
*huge
:tiny
*tiny
;
143 /* now |1-x| is tiny <= 2**-20, suffice to compute
144 log(x) by x-x^2/2+x^3/3-x^4/4 */
145 t
= x
-1; /* t has 20 trailing zeros */
146 w
= (t
*t
)*((float)0.5-t
*((float)0.333333333333-t
*(float)0.25));
147 u
= ivln2_h
*t
; /* ivln2_h has 16 sig. bits */
148 v
= t
*ivln2_l
-w
*ivln2
;
150 GET_FLOAT_WORD(is
,t1
);
151 SET_FLOAT_WORD(t1
,is
&0xfffff000);
154 float s2
,s_h
,s_l
,t_h
,t_l
;
156 /* take care subnormal number */
158 {ax
*= two24
; n
-= 24; GET_FLOAT_WORD(ix
,ax
); }
159 n
+= ((ix
)>>23)-0x7f;
161 /* determine interval */
162 ix
= j
|0x3f800000; /* normalize ix */
163 if(j
<=0x1cc471) k
=0; /* |x|<sqrt(3/2) */
164 else if(j
<0x5db3d7) k
=1; /* |x|<sqrt(3) */
165 else {k
=0;n
+=1;ix
-= 0x00800000;}
166 SET_FLOAT_WORD(ax
,ix
);
168 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
169 u
= ax
-bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
173 GET_FLOAT_WORD(is
,s_h
);
174 SET_FLOAT_WORD(s_h
,is
&0xfffff000);
175 /* t_h=ax+bp[k] High */
176 SET_FLOAT_WORD(t_h
,((ix
>>1)|0x20000000)+0x0040000+(k
<<21));
177 t_l
= ax
- (t_h
-bp
[k
]);
178 s_l
= v
*((u
-s_h
*t_h
)-s_h
*t_l
);
179 /* compute log(ax) */
181 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
184 t_h
= (float)3.0+s2
+r
;
185 GET_FLOAT_WORD(is
,t_h
);
186 SET_FLOAT_WORD(t_h
,is
&0xfffff000);
187 t_l
= r
-((t_h
-(float)3.0)-s2
);
188 /* u+v = s*(1+...) */
191 /* 2/(3log2)*(s+...) */
193 GET_FLOAT_WORD(is
,p_h
);
194 SET_FLOAT_WORD(p_h
,is
&0xfffff000);
196 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
197 z_l
= cp_l
*p_h
+p_l
*cp
+dp_l
[k
];
198 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
200 t1
= (((z_h
+z_l
)+dp_h
[k
])+t
);
201 GET_FLOAT_WORD(is
,t1
);
202 SET_FLOAT_WORD(t1
,is
&0xfffff000);
203 t2
= z_l
-(((t1
-t
)-dp_h
[k
])-z_h
);
206 s
= one
; /* s (sign of result -ve**odd) = -1 else = 1 */
207 if(((((u_int32_t
)hx
>>31)-1)|(yisint
-1))==0)
208 s
= -one
; /* (-ve)**(odd int) */
210 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
211 GET_FLOAT_WORD(is
,y
);
212 SET_FLOAT_WORD(y1
,is
&0xfffff000);
213 p_l
= (y
-y1
)*t1
+y
*t2
;
217 if (j
>0x43000000) /* if z > 128 */
218 return s
*huge
*huge
; /* overflow */
219 else if (j
==0x43000000) { /* if z == 128 */
220 if(p_l
+ovt
>z
-p_h
) return s
*huge
*huge
; /* overflow */
222 else if ((j
&0x7fffffff)>0x43160000) /* z <= -150 */
223 return s
*tiny
*tiny
; /* underflow */
224 else if ((u_int32_t
) j
==0xc3160000){ /* z == -150 */
225 if(p_l
<=z
-p_h
) return s
*tiny
*tiny
; /* underflow */
228 * compute 2**(p_h+p_l)
233 if(i
>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
234 n
= j
+(0x00800000>>(k
+1));
235 k
= ((n
&0x7fffffff)>>23)-0x7f; /* new k for n */
236 SET_FLOAT_WORD(t
,n
&~(0x007fffff>>k
));
237 n
= ((n
&0x007fffff)|0x00800000)>>(23-k
);
242 GET_FLOAT_WORD(is
,t
);
243 SET_FLOAT_WORD(t
,is
&0xfffff000);
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
-two
)-(w
+z
*w
);
254 if((j
>>23)<=0) z
= __scalbnf(z
,n
); /* subnormal output */
255 else SET_FLOAT_WORD(z
,j
);