1 /* Single-precision pow function.
2 Copyright (c) 2017-2018 Arm Ltd. All rights reserved.
3
4 SPDX-License-Identifier: BSD-3-Clause
5
6 Redistribution and use in source and binary forms, with or without
7 modification, are permitted provided that the following conditions
8 are met:
9 1. Redistributions of source code must retain the above copyright
10 notice, this list of conditions and the following disclaimer.
11 2. Redistributions in binary form must reproduce the above copyright
12 notice, this list of conditions and the following disclaimer in the
13 documentation and/or other materials provided with the distribution.
14 3. The name of the company may not be used to endorse or promote
15 products derived from this software without specific prior written
16 permission.
17
18 THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS'' AND ANY EXPRESS OR IMPLIED
19 WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
20 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
22 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
23 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
28
29 #include "fdlibm.h"
30 #if !__OBSOLETE_MATH_FLOAT
31
32 #include <math.h>
33 #include <stdint.h>
34 #include "math_config.h"
35
36 /*
37 POWF_LOG2_POLY_ORDER = 5
38 EXP2F_TABLE_BITS = 5
39
40 ULP error: 0.82 (~ 0.5 + relerr*2^24)
41 relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2)
42 relerr_log2: 1.83 * 2^-33 (Relative error of logx.)
43 relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).)
44 */
45
46 #define N (1 << POWF_LOG2_TABLE_BITS)
47 #define T __powf_log2_data.tab
48 #define A __powf_log2_data.poly
49 #define OFF 0x3f330000
50
51 /* Subnormal input is normalized so ix has negative biased exponent.
52 Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set. */
53 static inline double_t
log2_inline(uint32_t ix)54 log2_inline (uint32_t ix)
55 {
56 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
57 double_t z, r, r2, r4, p, q, y, y0, invc, logc;
58 uint32_t iz, top, tmp;
59 int k, i;
60
61 /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
62 The range is split into N subintervals.
63 The ith subinterval contains z and c is near its center. */
64 tmp = ix - OFF;
65 i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N;
66 top = tmp & 0xff800000;
67 iz = ix - top;
68 k = (int32_t) top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */
69 invc = T[i].invc;
70 logc = T[i].logc;
71 z = (double_t) asfloat (iz);
72
73 /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */
74 r = z * invc - 1;
75 y0 = logc + (double_t) k;
76
77 /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */
78 r2 = r * r;
79 y = A[0] * r + A[1];
80 p = A[2] * r + A[3];
81 r4 = r2 * r2;
82 q = A[4] * r + y0;
83 q = p * r2 + q;
84 y = y * r4 + q;
85 return y;
86 }
87
88 #undef N
89 #undef T
90 #define N (1 << EXP2F_TABLE_BITS)
91 #define T __exp2f_data.tab
92 #define SIGN_BIAS ((uint32_t) 1 << (EXP2F_TABLE_BITS + 11))
93
94 /* The output of log2 and thus the input of exp2 is either scaled by N
95 (in case of fast toint intrinsics) or not. The unscaled xd must be
96 in [-1021,1023], sign_bias sets the sign of the result. */
97 static inline double_t
exp2_inline(double_t xd,uint32_t sign_bias)98 exp2_inline (double_t xd, uint32_t sign_bias)
99 {
100 uint64_t ki, ski, t;
101 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
102 double_t kd, z, r, r2, y, s;
103
104 #if TOINT_INTRINSICS
105 # define C __exp2f_data.poly_scaled
106 /* N*x = k + r with r in [-1/2, 1/2] */
107 kd = roundtoint (xd); /* k */
108 ki = converttoint (xd);
109 #else
110 # define C __exp2f_data.poly
111 # define SHIFT __exp2f_data.shift_scaled
112 /* x = k/N + r with r in [-1/(2N), 1/(2N)] */
113 kd = (double) (xd + SHIFT); /* Rounding to double precision is required. */
114 ki = asuint64 (kd);
115 kd -= SHIFT; /* k/N */
116 #endif
117 r = xd - kd;
118
119 /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
120 t = T[ki % N];
121 ski = ki + sign_bias;
122 t += ski << (52 - EXP2F_TABLE_BITS);
123 s = asfloat64 (t);
124 z = C[0] * r + C[1];
125 r2 = r * r;
126 y = C[2] * r + 1;
127 y = z * r2 + y;
128 y = y * s;
129 return y;
130 }
131
132 /* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
133 the bit representation of a non-zero finite floating-point value. */
134 static inline int
checkint(uint32_t iy)135 checkint (uint32_t iy)
136 {
137 int e = iy >> 23 & 0xff;
138 if (e < 0x7f)
139 return 0;
140 if (e > 0x7f + 23)
141 return 2;
142 if (iy & (((uint32_t) 1 << (0x7f + 23 - e)) - 1))
143 return 0;
144 if (iy & ((uint32_t) 1 << (0x7f + 23 - e)))
145 return 1;
146 return 2;
147 }
148
149 static inline int
zeroinfnan(uint32_t ix)150 zeroinfnan (uint32_t ix)
151 {
152 return 2 * ix - 1 >= 2u * (uint32_t) 0x7f800000 - 1;
153 }
154
155 float
powf(float x,float y)156 powf (float x, float y)
157 {
158 uint32_t sign_bias = 0;
159 uint32_t ix, iy;
160
161 ix = asuint (x);
162 iy = asuint (y);
163 if (__builtin_expect (ix - 0x00800000 >= 0x7f800000 - 0x00800000
164 || zeroinfnan (iy),
165 0))
166 {
167 /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */
168 if (__builtin_expect (zeroinfnan (iy), 0))
169 {
170 if (2 * iy == 0)
171 return issignalingf_inline (x) ? x + y : 1.0f;
172 if (ix == 0x3f800000)
173 return issignalingf_inline (y) ? x + y : 1.0f;
174 if (2 * ix > 2u * (uint32_t) 0x7f800000 || 2 * iy > 2u * (uint32_t) 0x7f800000)
175 return x + y;
176 if (2 * ix == 2 * (uint32_t) 0x3f800000)
177 return 1.0f;
178 if ((2 * ix < 2 * (uint32_t) 0x3f800000) == !(iy & (uint32_t) 0x80000000))
179 return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
180 return y * y;
181 }
182 if (__builtin_expect (zeroinfnan (ix), 0))
183 {
184 float_t x2 = x * x;
185 if (ix & 0x80000000 && checkint (iy) == 1)
186 {
187 x2 = -x2;
188 sign_bias = 1;
189 }
190 if (!(iy & 0x80000000))
191 return opt_barrier_float(x2);
192 #if WANT_ERRNO
193 if (2 * ix == 0)
194 return __math_divzerof (sign_bias);
195 #endif
196 return 1 / x2;
197 }
198 /* x and y are non-zero finite. */
199 if (ix & 0x80000000)
200 {
201 /* Finite x < 0. */
202 int yint = checkint (iy);
203 if (yint == 0)
204 return __math_invalidf (x);
205 if (yint == 1)
206 sign_bias = SIGN_BIAS;
207 ix &= 0x7fffffff;
208 }
209 if (ix < 0x00800000)
210 {
211 /* Normalize subnormal x so exponent becomes negative. */
212 ix = asuint (x * 0x1p23f);
213 ix &= 0x7fffffff;
214 ix -= (uint32_t) 23 << 23;
215 }
216 }
217 double_t logx = log2_inline (ix);
218 double_t ylogx = (double) y * logx; /* Note: cannot overflow, y is single prec. */
219 if (__builtin_expect ((asuint64 (ylogx) >> 47 & 0xffff)
220 >= asuint64 (126.0 * POWF_SCALE) >> 47,
221 0))
222 {
223 /* |y*log(x)| >= 126. */
224 if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE)
225 /* |x^y| > 0x1.ffffffp127. */
226 return __math_oflowf (sign_bias);
227 if (WANT_ROUNDING && WANT_ERRNO
228 && ylogx > 0x1.fffffffa3aae2p+6 * POWF_SCALE)
229 /* |x^y| > 0x1.fffffep127, check if we round away from 0. */
230 if ((!sign_bias
231 && eval_as_float (1.0f + opt_barrier_float (0x1p-25f)) != 1.0f)
232 || (sign_bias
233 && eval_as_float (-1.0f - opt_barrier_float (0x1p-25f))
234 != -1.0f))
235 return __math_oflowf (sign_bias);
236 if (ylogx <= -150.0 * POWF_SCALE)
237 return __math_uflowf (sign_bias);
238 #if WANT_ERRNO_UFLOW
239 if (ylogx < -149.0 * POWF_SCALE)
240 return __math_may_uflowf (sign_bias);
241 #endif
242 }
243 return (float) exp2_inline (ylogx, sign_bias);
244 }
245
246 #if defined(_HAVE_ALIAS_ATTRIBUTE)
247 #ifndef __clang__
248 #pragma GCC diagnostic ignored "-Wmissing-attributes"
249 #endif
250 __strong_reference(powf, _powf);
251 #endif
252
253 #endif /* !__OBSOLETE_MATH_FLOAT */
254
