1 /* ef_pow.c -- float version of e_pow.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
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
12 * is preserved.
13 * ====================================================
14 */
15
16 #include "fdlibm.h"
17 #include "math_config.h"
18
19 #if __OBSOLETE_MATH_FLOAT
20
21 #ifdef __v810__
22 #define const
23 #endif
24
25 static const float
26 bp[] = {1.0, 1.5,},
27 dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
28 dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
29 zero = 0.0,
30 one = 1.0,
31 two = 2.0,
32 two24 = 16777216.0, /* 0x4b800000 */
33 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
34 L1 = 6.0000002384e-01, /* 0x3f19999a */
35 L2 = 4.2857143283e-01, /* 0x3edb6db7 */
36 L3 = 3.3333334327e-01, /* 0x3eaaaaab */
37 L4 = 2.7272811532e-01, /* 0x3e8ba305 */
38 L5 = 2.3066075146e-01, /* 0x3e6c3255 */
39 L6 = 2.0697501302e-01, /* 0x3e53f142 */
40 P1 = 1.6666667163e-01, /* 0x3e2aaaab */
41 P2 = -2.7777778450e-03, /* 0xbb360b61 */
42 P3 = 6.6137559770e-05, /* 0x388ab355 */
43 P4 = -1.6533901999e-06, /* 0xb5ddea0e */
44 P5 = 4.1381369442e-08, /* 0x3331bb4c */
45 lg2 = 6.9314718246e-01, /* 0x3f317218 */
46 lg2_h = 6.93145752e-01, /* 0x3f317200 */
47 lg2_l = 1.42860654e-06, /* 0x35bfbe8c */
48 ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
49 cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
50 cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */
51 cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
52 ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
53 ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
54 ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
55
56 float
powf(float x,float y)57 powf(float x, float y)
58 {
59 float z, ax, z_h, z_l, p_h, p_l;
60 float y1, t1, t2, r, s, t, u, v, w;
61 __int32_t i, j, k, yisint, n;
62 __int32_t hx, hy, ix, iy, is;
63
64 GET_FLOAT_WORD(hx, x);
65 GET_FLOAT_WORD(hy, y);
66 ix = hx & 0x7fffffff;
67 iy = hy & 0x7fffffff;
68
69 /* y==zero: x**0 = 1 */
70 if (FLT_UWORD_IS_ZERO(iy)) {
71 if (issignalingf_inline(x))
72 return x + y;
73 return one;
74 }
75
76 /* x|y==NaN return NaN unless x==1 then return 1 */
77 if (FLT_UWORD_IS_NAN(ix) || FLT_UWORD_IS_NAN(iy)) {
78 if (hx == 0x3f800000 && !issignalingf_inline(y))
79 return one;
80 else
81 return x + y;
82 }
83
84 /* determine if y is an odd int when x < 0
85 * yisint = 0 ... y is not an integer
86 * yisint = 1 ... y is an odd int
87 * yisint = 2 ... y is an even int
88 */
89 yisint = 0;
90 if (hx < 0) {
91 if (iy >= 0x4b800000)
92 yisint = 2; /* even integer y */
93 else if (iy >= 0x3f800000) {
94 k = (iy >> 23) - 0x7f; /* exponent */
95 j = iy >> (23 - k);
96 if ((j << (23 - k)) == iy)
97 yisint = 2 - (j & 1);
98 }
99 }
100
101 /* special value of y */
102 if (FLT_UWORD_IS_INFINITE(iy)) { /* y is +-inf */
103 if (ix == 0x3f800000)
104 return one; /* +-1**+-inf = 1 */
105 else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
106 return (hy >= 0) ? y : zero;
107 else /* (|x|<1)**-,+inf = inf,0 */
108 return (hy < 0) ? -y : zero;
109 }
110 if (iy == 0x3f800000) { /* y is +-1 */
111 if (hy < 0) {
112 if (x == 0)
113 return __math_divzerof(hx < 0);
114 return one / x;
115 } else
116 return x;
117 }
118 if (hy == 0x40000000 && ix < 0x5f800000 && ix > 0x1a000000)
119 return x * x; /* y is 2 */
120 if (hy == 0x3f000000) { /* y is 0.5 */
121 if (hx >= 0) /* x >= +0 */
122 return sqrtf(x);
123 }
124
125 ax = fabsf(x);
126 /* special value of x */
127 if (FLT_UWORD_IS_INFINITE(ix) || ix == 0x3f800000) {
128 z = ax; /*x is +-inf,+-1*/
129 if (hy < 0)
130 z = one / z; /* z = (1/|x|) */
131 if (hx < 0) {
132 if (((ix - 0x3f800000) | yisint) == 0) {
133 return __math_invalidf(x); /* (-1)**non-int is NaN */
134 } else if (yisint == 1)
135 z = -z; /* (x<0)**odd = -(|x|**odd) */
136 }
137 return z;
138 }
139
140 if (FLT_UWORD_IS_ZERO(ix)) {
141 if (hy < 0)
142 return __math_divzerof(hx < 0 && yisint == 1);
143 if (yisint != 1)
144 x = ax;
145 return x;
146 }
147
148 /* (x<0)**(non-int) is NaN */
149 if (((((__uint32_t)hx >> 31) - 1) | yisint) == 0)
150 return __math_invalidf(x);
151
152 /* |y| is huge */
153 if (iy > 0x4d000000) { /* if |y| > 2**27 */
154 /* over/underflow if x is not close to one */
155 if (ix < 0x3f7ffff4)
156 return (hy < 0) ? __math_oflowf(0) : __math_uflowf(0);
157 if (ix > 0x3f800007)
158 return (hy > 0) ? __math_oflowf(0) : __math_uflowf(0);
159 /* now |1-x| is tiny <= 2**-20, suffice to compute
160 log(x) by x-x^2/2+x^3/3-x^4/4 */
161 t = ax - 1; /* t has 20 trailing zeros */
162 w = (t * t) *
163 ((float)0.5 - t * ((float)0.333333333333 - t * (float)0.25));
164 u = ivln2_h * t; /* ivln2_h has 16 sig. bits */
165 v = t * ivln2_l - w * ivln2;
166 t1 = u + v;
167 GET_FLOAT_WORD(is, t1);
168 SET_FLOAT_WORD(t1, is & 0xfffff000);
169 t2 = v - (t1 - u);
170 } else {
171 float s2, s_h, s_l, t_h, t_l;
172 n = 0;
173 /* take care subnormal number */
174 if (FLT_UWORD_IS_SUBNORMAL(ix)) {
175 ax *= two24;
176 n -= 24;
177 GET_FLOAT_WORD(ix, ax);
178 }
179 n += ((ix) >> 23) - 0x7f;
180 j = ix & 0x007fffff;
181 /* determine interval */
182 ix = j | 0x3f800000; /* normalize ix */
183 if (j <= 0x1cc471)
184 k = 0; /* |x|<sqrt(3/2) */
185 else if (j < 0x5db3d7)
186 k = 1; /* |x|<sqrt(3) */
187 else {
188 k = 0;
189 n += 1;
190 ix -= 0x00800000;
191 }
192 SET_FLOAT_WORD(ax, ix);
193
194 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
195 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
196 v = one / (ax + bp[k]);
197 s = u * v;
198 s_h = s;
199 GET_FLOAT_WORD(is, s_h);
200 SET_FLOAT_WORD(s_h, is & 0xfffff000);
201 /* t_h=ax+bp[k] High */
202 SET_FLOAT_WORD(t_h, ((ix >> 1) | 0x20000000) + 0x0040000 + (k << 21));
203 t_l = ax - (t_h - bp[k]);
204 s_l = v * ((u - s_h * t_h) - s_h * t_l);
205 /* compute log(ax) */
206 s2 = s * s;
207 r = s2 * s2 *
208 (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
209 r += s_l * (s_h + s);
210 s2 = s_h * s_h;
211 t_h = (float)3.0 + s2 + r;
212 GET_FLOAT_WORD(is, t_h);
213 SET_FLOAT_WORD(t_h, is & 0xfffff000);
214 t_l = r - ((t_h - (float)3.0) - s2);
215 /* u+v = s*(1+...) */
216 u = s_h * t_h;
217 v = s_l * t_h + t_l * s;
218 /* 2/(3log2)*(s+...) */
219 p_h = u + v;
220 GET_FLOAT_WORD(is, p_h);
221 SET_FLOAT_WORD(p_h, is & 0xfffff000);
222 p_l = v - (p_h - u);
223 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
224 z_l = cp_l * p_h + p_l * cp + dp_l[k];
225 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
226 t = (float)n;
227 t1 = (((z_h + z_l) + dp_h[k]) + t);
228 GET_FLOAT_WORD(is, t1);
229 SET_FLOAT_WORD(t1, is & 0xfffff000);
230 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
231 }
232
233 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
234 if (((((__uint32_t)hx >> 31) - 1) | (yisint - 1)) == 0)
235 s = -one; /* (-ve)**(odd int) */
236
237 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
238 GET_FLOAT_WORD(is, y);
239 SET_FLOAT_WORD(y1, is & 0xfffff000);
240 p_l = (y - y1) * t1 + y * t2;
241 p_h = y1 * t1;
242 z = p_l + p_h;
243 GET_FLOAT_WORD(j, z);
244 i = j & 0x7fffffff;
245 if (j > 0) {
246 if (i > FLT_UWORD_EXP_MAX)
247 return __math_oflowf(s < 0); /* overflow */
248 else if (i == FLT_UWORD_EXP_MAX)
249 if (p_l + ovt > z - p_h)
250 return __math_oflowf(s < 0); /* overflow */
251 } else {
252 if (i > FLT_UWORD_EXP_MIN)
253 return __math_uflowf(s < 0); /* underflow */
254 else if (i == FLT_UWORD_EXP_MIN)
255 if (p_l <= z - p_h)
256 return __math_uflowf(s < 0); /* underflow */
257 }
258 /*
259 * compute 2**(p_h+p_l)
260 */
261 k = (i >> 23) - 0x7f;
262 n = 0;
263 if (i > 0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
264 n = j + (0x00800000 >> (k + 1));
265 k = ((n & 0x7fffffff) >> 23) - 0x7f; /* new k for n */
266 SET_FLOAT_WORD(t, n & ~(0x007fffff >> k));
267 n = ((n & 0x007fffff) | 0x00800000) >> (23 - k);
268 if (j < 0)
269 n = -n;
270 p_h -= t;
271 }
272 t = p_l + p_h;
273 GET_FLOAT_WORD(is, t);
274 SET_FLOAT_WORD(t, is & 0xfffff000);
275 u = t * lg2_h;
276 v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
277 z = u + v;
278 w = v - (z - u);
279 t = z * z;
280 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
281 r = (z * t1) / (t1 - two) - (w + z * w);
282 z = one - (r - z);
283 GET_FLOAT_WORD(j, z);
284 j += (n << 23);
285 if ((j >> 23) <= 0)
286 z = scalbnf(z, (int)n); /* subnormal output */
287 else
288 SET_FLOAT_WORD(z, j);
289 return s * z;
290 }
291
292 #if defined(_HAVE_ALIAS_ATTRIBUTE)
293 #ifndef __clang__
294 #pragma GCC diagnostic ignored "-Wmissing-attributes"
295 #endif
296 __strong_reference(powf, _powf);
297 #endif
298
299 _MATH_ALIAS_f_ff(pow)
300
301 #else
302 #include "../common/sf_pow.c"
303 #endif /* __OBSOLETE_MATH_FLOAT */
304
