1/* e_powf.c -- float version of e_pow.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4/* Copyright (C) 2017 Free Software Foundation, Inc.
5 This file is part of the GNU C Library.
6
7 The GNU C Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public
9 License as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11
12 The GNU C Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Lesser General Public License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with the GNU C Library; if not, see
19 <http://www.gnu.org/licenses/>. */
20
21/*
22 * ====================================================
23 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
24 *
25 * Developed at SunPro, a Sun Microsystems, Inc. business.
26 * Permission to use, copy, modify, and distribute this
27 * software is freely granted, provided that this notice
28 * is preserved.
29 * ====================================================
30 */
31
32#include <math.h>
33#include <math_private.h>
34
35static const float huge = 1.0e+30, tiny = 1.0e-30;
36
37static const float
38bp[] = {1.0, 1.5,},
39zero = 0.0,
40one = 1.0,
41two = 2.0,
42two24 = 16777216.0, /* 0x4b800000 */
43 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
44L1 = 6.0000002384e-01, /* 0x3f19999a */
45L2 = 4.2857143283e-01, /* 0x3edb6db7 */
46L3 = 3.3333334327e-01, /* 0x3eaaaaab */
47L4 = 2.7272811532e-01, /* 0x3e8ba305 */
48L5 = 2.3066075146e-01, /* 0x3e6c3255 */
49L6 = 2.0697501302e-01, /* 0x3e53f142 */
50P1 = 1.6666667163e-01, /* 0x3e2aaaab */
51P2 = -2.7777778450e-03, /* 0xbb360b61 */
52P3 = 6.6137559770e-05, /* 0x388ab355 */
53P4 = -1.6533901999e-06, /* 0xb5ddea0e */
54P5 = 4.1381369442e-08, /* 0x3331bb4c */
55ovt = 4.2995665694e-08; /* -(128-log2(ovfl+.5ulp)) */
56
57static const double
58 dp[] = { 0.0, 0x1.2b803473f7ad1p-1, }, /* log2(1.5) */
59 lg2 = M_LN2,
60 cp = 2.0/3.0/M_LN2,
61 invln2 = 1.0/M_LN2;
62
63float
64__ieee754_powf(float x, float y)
65{
66 float z, ax, s;
67 double d1, d2;
68 int32_t i,j,k,yisint,n;
69 int32_t hx,hy,ix,iy;
70
71 GET_FLOAT_WORD(hy,y);
72 iy = hy&0x7fffffff;
73
74 /* y==zero: x**0 = 1 */
75 if(iy==0 && !issignaling (x)) return one;
76
77 /* x==+-1 */
78 if(x == 1.0 && !issignaling (y)) return one;
79 if(x == -1.0 && isinf(y)) return one;
80
81 GET_FLOAT_WORD(hx,x);
82 ix = hx&0x7fffffff;
83
84 /* +-NaN return x+y */
85 if(__builtin_expect(ix > 0x7f800000 ||
86 iy > 0x7f800000, 0))
87 return x+y;
88
89 /* special value of y */
90 if (__builtin_expect(iy==0x7f800000, 0)) { /* y is +-inf */
91 if (ix==0x3f800000)
92 return y - y; /* inf**+-1 is NaN */
93 else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
94 return (hy>=0)? y: zero;
95 else /* (|x|<1)**-,+inf = inf,0 */
96 return (hy<0)?-y: zero;
97 }
98 if(iy==0x3f800000) { /* y is +-1 */
99 if(hy<0) return one/x; else return x;
100 }
101 if(hy==0x40000000) return x*x; /* y is 2 */
102 if(hy==0x3f000000) { /* y is 0.5 */
103 if(__builtin_expect(hx>=0, 1)) /* x >= +0 */
104 return __ieee754_sqrtf(x);
105 }
106
107 /* determine if y is an odd int when x < 0
108 * yisint = 0 ... y is not an integer
109 * yisint = 1 ... y is an odd int
110 * yisint = 2 ... y is an even int
111 */
112 yisint = 0;
113 if(hx<0) {
114 if(iy>=0x4b800000) yisint = 2; /* even integer y */
115 else if(iy>=0x3f800000) {
116 k = (iy>>23)-0x7f; /* exponent */
117 j = iy>>(23-k);
118 if((j<<(23-k))==iy) yisint = 2-(j&1);
119 }
120 }
121
122 ax = fabsf(x);
123 /* special value of x */
124 if(__builtin_expect(ix==0x7f800000||ix==0||ix==0x3f800000, 0)){
125 z = ax; /*x is +-0,+-inf,+-1*/
126 if(hy<0) z = one/z; /* z = (1/|x|) */
127 if(hx<0) {
128 if(((ix-0x3f800000)|yisint)==0) {
129 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
130 } else if(yisint==1)
131 z = -z; /* (x<0)**odd = -(|x|**odd) */
132 }
133 return z;
134 }
135
136 /* (x<0)**(non-int) is NaN */
137 if(__builtin_expect(((((u_int32_t)hx>>31)-1)|yisint)==0, 0))
138 return (x-x)/(x-x);
139
140 /* |y| is huge */
141 if(__builtin_expect(iy>0x4d000000, 0)) { /* if |y| > 2**27 */
142 /* over/underflow if x is not close to one */
143 if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny;
144 if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
145 /* now |1-x| is tiny <= 2**-20, suffice to compute
146 log(x) by x-x^2/2+x^3/3-x^4/4 */
147 d2 = ax-1; /* d2 has 20 trailing zeros. */
148 d2 = d2 * invln2 -
149 (d2 * d2) * (0.5 - d2 * (0.333333333333 - d2 * 0.25)) * invln2;
150 } else {
151 /* Avoid internal underflow for tiny y. The exact value
152 of y does not matter if |y| <= 2**-32. */
153 if (iy < 0x2f800000)
154 SET_FLOAT_WORD (y, (hy & 0x80000000) | 0x2f800000);
155 n = 0;
156 /* take care subnormal number */
157 if(ix<0x00800000)
158 {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
159 n += ((ix)>>23)-0x7f;
160 j = ix&0x007fffff;
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);
167
168 /* compute d1 = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
169 d1 = (ax-(double)bp[k])/(ax+(double)bp[k]);
170 /* compute d2 = log(ax) */
171 d2 = d1 * d1;
172 d2 = 3.0 + d2 + d2*d2*(L1+d2*(L2+d2*(L3+d2*(L4+d2*(L5+d2*L6)))));
173 /* 2/(3log2)*(d2+...) */
174 d2 = d1*d2*cp;
175 /* log2(ax) = (d2+..)*2/(3*log2) */
176 d2 = d2+dp[k]+(double)n;
177 }
178
179 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
180 if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
181 s = -one; /* (-ve)**(odd int) */
182
183 /* compute y * d2 */
184 d1 = y * d2;
185 z = d1;
186 GET_FLOAT_WORD(j,z);
187 if (__builtin_expect(j>0x43000000, 0)) /* if z > 128 */
188 return s*huge*huge; /* overflow */
189 else if (__builtin_expect(j==0x43000000, 0)) { /* if z == 128 */
190 if(ovt>(z-d1)) return s*huge*huge; /* overflow */
191 }
192 else if (__builtin_expect((j&0x7fffffff)>0x43160000, 0))/* z <= -150 */
193 return s*tiny*tiny; /* underflow */
194 else if (__builtin_expect((u_int32_t) j==0xc3160000, 0)){/* z == -150*/
195 if(0.0<=(z-d1)) return s*tiny*tiny; /* underflow */
196 }
197 /*
198 * compute 2**d1
199 */
200 i = j&0x7fffffff;
201 k = (i>>23)-0x7f;
202 n = 0;
203 if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
204 n = j+(0x00800000>>(k+1));
205 k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */
206 SET_FLOAT_WORD(z,n&~(0x007fffff>>k));
207 n = ((n&0x007fffff)|0x00800000)>>(23-k);
208 if(j<0) n = -n;
209 d1 -= z;
210 }
211 d1 = d1 * lg2;
212 d2 = d1*d1;
213 d2 = d1 - d2*(P1+d2*(P2+d2*(P3+d2*(P4+d2*P5))));
214 d2 = (d1*d2)/(d2-two);
215 z = one - (d2-d1);
216 GET_FLOAT_WORD(j,z);
217 j += (n<<23);
218 if((j>>23)<=0) /* subnormal output */
219 {
220 z = __scalbnf (z, n);
221 float force_underflow = z * z;
222 math_force_eval (force_underflow);
223 }
224 else SET_FLOAT_WORD(z,j);
225 return s*z;
226}
227strong_alias (__ieee754_powf, __powf_finite)
228