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.61914062e-01, /* 0x3f764000 =head of cp */
51 cp_l = -1.17368574e-04, /* 0xb8f623c6 =tail of cp */
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 is = ((ix >> 1) & 0xfffff000U) | 0x20000000;
203 SET_FLOAT_WORD(t_h, is + 0x00400000 + (k << 21));
204 t_l = ax - (t_h - bp[k]);
205 s_l = v * ((u - s_h * t_h) - s_h * t_l);
206 /* compute log(ax) */
207 s2 = s * s;
208 r = s2 * s2 *
209 (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
210 r += s_l * (s_h + s);
211 s2 = s_h * s_h;
212 t_h = (float)3.0 + s2 + r;
213 GET_FLOAT_WORD(is, t_h);
214 SET_FLOAT_WORD(t_h, is & 0xfffff000);
215 t_l = r - ((t_h - (float)3.0) - s2);
216 /* u+v = s*(1+...) */
217 u = s_h * t_h;
218 v = s_l * t_h + t_l * s;
219 /* 2/(3log2)*(s+...) */
220 p_h = u + v;
221 GET_FLOAT_WORD(is, p_h);
222 SET_FLOAT_WORD(p_h, is & 0xfffff000);
223 p_l = v - (p_h - u);
224 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
225 z_l = cp_l * p_h + p_l * cp + dp_l[k];
226 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
227 t = (float)n;
228 t1 = (((z_h + z_l) + dp_h[k]) + t);
229 GET_FLOAT_WORD(is, t1);
230 SET_FLOAT_WORD(t1, is & 0xfffff000);
231 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
232 }
233
234 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
235 if (((((__uint32_t)hx >> 31) - 1) | (yisint - 1)) == 0)
236 s = -one; /* (-ve)**(odd int) */
237
238 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
239 GET_FLOAT_WORD(is, y);
240 SET_FLOAT_WORD(y1, is & 0xfffff000);
241 p_l = (y - y1) * t1 + y * t2;
242 p_h = y1 * t1;
243 z = p_l + p_h;
244 GET_FLOAT_WORD(j, z);
245 i = j & 0x7fffffff;
246 if (j > 0) {
247 if (i > FLT_UWORD_EXP_MAX)
248 return __math_oflowf(s < 0); /* overflow */
249 else if (i == FLT_UWORD_EXP_MAX)
250 if (p_l + ovt > z - p_h)
251 return __math_oflowf(s < 0); /* overflow */
252 } else {
253 if (i > FLT_UWORD_EXP_MIN)
254 return __math_uflowf(s < 0); /* underflow */
255 else if (i == FLT_UWORD_EXP_MIN)
256 if (p_l <= z - p_h)
257 return __math_uflowf(s < 0); /* underflow */
258 }
259 /*
260 * compute 2**(p_h+p_l)
261 */
262 k = (i >> 23) - 0x7f;
263 n = 0;
264 if (i > 0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
265 n = j + (0x00800000 >> (k + 1));
266 k = ((n & 0x7fffffff) >> 23) - 0x7f; /* new k for n */
267 SET_FLOAT_WORD(t, n & ~(0x007fffff >> k));
268 n = ((n & 0x007fffff) | 0x00800000) >> (23 - k);
269 if (j < 0)
270 n = -n;
271 p_h -= t;
272 }
273 t = p_l + p_h;
274 GET_FLOAT_WORD(is, t);
275 SET_FLOAT_WORD(t, is & 0xfffff000);
276 u = t * lg2_h;
277 v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
278 z = u + v;
279 w = v - (z - u);
280 t = z * z;
281 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
282 r = (z * t1) / (t1 - two) - (w + z * w);
283 z = one - (r - z);
284 GET_FLOAT_WORD(j, z);
285 j += (n << 23);
286 if ((j >> 23) <= 0)
287 z = scalbnf(z, (int)n); /* subnormal output */
288 else
289 SET_FLOAT_WORD(z, j);
290 return s * z;
291 }
292
293 #if defined(_HAVE_ALIAS_ATTRIBUTE)
294 #ifndef __clang__
295 #pragma GCC diagnostic ignored "-Wmissing-attributes"
296 #endif
297 __strong_reference(powf, _powf);
298 #endif
299
300 _MATH_ALIAS_f_ff(pow)
301
302 #else
303 #include "../common/sf_pow.c"
304 #endif /* __OBSOLETE_MATH_FLOAT */
305
