]>
Commit | Line | Data |
---|---|---|
4ea49f4c | 1 | /* Single-precision pow function. |
04277e02 | 2 | Copyright (C) 2017-2019 Free Software Foundation, Inc. |
c340290d PC |
3 | This file is part of the GNU C Library. |
4 | ||
5 | The GNU C Library is free software; you can redistribute it and/or | |
6 | modify it under the terms of the GNU Lesser General Public | |
7 | License as published by the Free Software Foundation; either | |
8 | version 2.1 of the License, or (at your option) any later version. | |
9 | ||
10 | The GNU C Library is distributed in the hope that it will be useful, | |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 | Lesser General Public License for more details. | |
14 | ||
15 | You should have received a copy of the GNU Lesser General Public | |
16 | License along with the GNU C Library; if not, see | |
5a82c748 | 17 | <https://www.gnu.org/licenses/>. */ |
f7eac6eb | 18 | |
1ed0291c | 19 | #include <math.h> |
505b5b29 SN |
20 | #include <math-barriers.h> |
21 | #include <math-narrow-eval.h> | |
4ea49f4c | 22 | #include <stdint.h> |
bd4430c2 | 23 | #include <shlib-compat.h> |
24b6515d | 24 | #include <libm-alias-float.h> |
4ea49f4c | 25 | #include "math_config.h" |
f7eac6eb | 26 | |
4ea49f4c SN |
27 | /* |
28 | POWF_LOG2_POLY_ORDER = 5 | |
29 | EXP2F_TABLE_BITS = 5 | |
30 | ||
31 | ULP error: 0.82 (~ 0.5 + relerr*2^24) | |
32 | relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2) | |
33 | relerr_log2: 1.83 * 2^-33 (Relative error of logx.) | |
34 | relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).) | |
35 | */ | |
36 | ||
37 | #define N (1 << POWF_LOG2_TABLE_BITS) | |
38 | #define T __powf_log2_data.tab | |
39 | #define A __powf_log2_data.poly | |
40 | #define OFF 0x3f330000 | |
41 | ||
42 | /* Subnormal input is normalized so ix has negative biased exponent. | |
43 | Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set. */ | |
44 | static inline double_t | |
45 | log2_inline (uint32_t ix) | |
f7eac6eb | 46 | { |
4ea49f4c SN |
47 | /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ |
48 | double_t z, r, r2, r4, p, q, y, y0, invc, logc; | |
49 | uint32_t iz, top, tmp; | |
50 | int k, i; | |
51 | ||
52 | /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. | |
53 | The range is split into N subintervals. | |
54 | The ith subinterval contains z and c is near its center. */ | |
55 | tmp = ix - OFF; | |
56 | i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N; | |
57 | top = tmp & 0xff800000; | |
58 | iz = ix - top; | |
59 | k = (int32_t) top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */ | |
60 | invc = T[i].invc; | |
61 | logc = T[i].logc; | |
62 | z = (double_t) asfloat (iz); | |
63 | ||
64 | /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ | |
65 | r = z * invc - 1; | |
66 | y0 = logc + (double_t) k; | |
67 | ||
68 | /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ | |
69 | r2 = r * r; | |
70 | y = A[0] * r + A[1]; | |
71 | p = A[2] * r + A[3]; | |
72 | r4 = r2 * r2; | |
73 | q = A[4] * r + y0; | |
74 | q = p * r2 + q; | |
75 | y = y * r4 + q; | |
76 | return y; | |
77 | } | |
f7eac6eb | 78 | |
4ea49f4c SN |
79 | #undef N |
80 | #undef T | |
81 | #define N (1 << EXP2F_TABLE_BITS) | |
82 | #define T __exp2f_data.tab | |
83 | #define SIGN_BIAS (1 << (EXP2F_TABLE_BITS + 11)) | |
84 | ||
85 | /* The output of log2 and thus the input of exp2 is either scaled by N | |
86 | (in case of fast toint intrinsics) or not. The unscaled xd must be | |
87 | in [-1021,1023], sign_bias sets the sign of the result. */ | |
88 | static inline double_t | |
2b445206 | 89 | exp2_inline (double_t xd, uint32_t sign_bias) |
4ea49f4c SN |
90 | { |
91 | uint64_t ki, ski, t; | |
92 | /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ | |
93 | double_t kd, z, r, r2, y, s; | |
94 | ||
95 | #if TOINT_INTRINSICS | |
96 | # define C __exp2f_data.poly_scaled | |
97 | /* N*x = k + r with r in [-1/2, 1/2] */ | |
98 | kd = roundtoint (xd); /* k */ | |
99 | ki = converttoint (xd); | |
100 | #else | |
101 | # define C __exp2f_data.poly | |
102 | # define SHIFT __exp2f_data.shift_scaled | |
103 | /* x = k/N + r with r in [-1/(2N), 1/(2N)] */ | |
104 | kd = (double) (xd + SHIFT); /* Rounding to double precision is required. */ | |
105 | ki = asuint64 (kd); | |
106 | kd -= SHIFT; /* k/N */ | |
107 | #endif | |
108 | r = xd - kd; | |
109 | ||
110 | /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ | |
111 | t = T[ki % N]; | |
112 | ski = ki + sign_bias; | |
113 | t += ski << (52 - EXP2F_TABLE_BITS); | |
114 | s = asdouble (t); | |
115 | z = C[0] * r + C[1]; | |
116 | r2 = r * r; | |
117 | y = C[2] * r + 1; | |
118 | y = z * r2 + y; | |
119 | y = y * s; | |
120 | return y; | |
121 | } | |
c340290d | 122 | |
d7347278 SN |
123 | /* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is |
124 | the bit representation of a non-zero finite floating-point value. */ | |
4ea49f4c SN |
125 | static inline int |
126 | checkint (uint32_t iy) | |
127 | { | |
128 | int e = iy >> 23 & 0xff; | |
129 | if (e < 0x7f) | |
130 | return 0; | |
131 | if (e > 0x7f + 23) | |
132 | return 2; | |
133 | if (iy & ((1 << (0x7f + 23 - e)) - 1)) | |
134 | return 0; | |
135 | if (iy & (1 << (0x7f + 23 - e))) | |
136 | return 1; | |
137 | return 2; | |
138 | } | |
ba1ffaa1 | 139 | |
4ea49f4c SN |
140 | static inline int |
141 | zeroinfnan (uint32_t ix) | |
142 | { | |
143 | return 2 * ix - 1 >= 2u * 0x7f800000 - 1; | |
144 | } | |
f7eac6eb | 145 | |
4ea49f4c | 146 | float |
bd4430c2 | 147 | __powf (float x, float y) |
4ea49f4c | 148 | { |
2b445206 | 149 | uint32_t sign_bias = 0; |
4ea49f4c SN |
150 | uint32_t ix, iy; |
151 | ||
152 | ix = asuint (x); | |
153 | iy = asuint (y); | |
154 | if (__glibc_unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000 | |
155 | || zeroinfnan (iy))) | |
156 | { | |
157 | /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */ | |
158 | if (__glibc_unlikely (zeroinfnan (iy))) | |
159 | { | |
160 | if (2 * iy == 0) | |
161 | return issignalingf_inline (x) ? x + y : 1.0f; | |
162 | if (ix == 0x3f800000) | |
163 | return issignalingf_inline (y) ? x + y : 1.0f; | |
164 | if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) | |
165 | return x + y; | |
166 | if (2 * ix == 2 * 0x3f800000) | |
167 | return 1.0f; | |
168 | if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) | |
169 | return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ | |
170 | return y * y; | |
171 | } | |
172 | if (__glibc_unlikely (zeroinfnan (ix))) | |
173 | { | |
174 | float_t x2 = x * x; | |
175 | if (ix & 0x80000000 && checkint (iy) == 1) | |
176 | { | |
177 | x2 = -x2; | |
178 | sign_bias = 1; | |
179 | } | |
180 | #if WANT_ERRNO | |
181 | if (2 * ix == 0 && iy & 0x80000000) | |
182 | return __math_divzerof (sign_bias); | |
183 | #endif | |
184 | return iy & 0x80000000 ? 1 / x2 : x2; | |
f7eac6eb | 185 | } |
4ea49f4c SN |
186 | /* x and y are non-zero finite. */ |
187 | if (ix & 0x80000000) | |
188 | { | |
189 | /* Finite x < 0. */ | |
190 | int yint = checkint (iy); | |
191 | if (yint == 0) | |
192 | return __math_invalidf (x); | |
193 | if (yint == 1) | |
194 | sign_bias = SIGN_BIAS; | |
195 | ix &= 0x7fffffff; | |
f7eac6eb | 196 | } |
4ea49f4c SN |
197 | if (ix < 0x00800000) |
198 | { | |
199 | /* Normalize subnormal x so exponent becomes negative. */ | |
200 | ix = asuint (x * 0x1p23f); | |
201 | ix &= 0x7fffffff; | |
202 | ix -= 23 << 23; | |
ba1ffaa1 | 203 | } |
4ea49f4c SN |
204 | } |
205 | double_t logx = log2_inline (ix); | |
206 | double_t ylogx = y * logx; /* Note: cannot overflow, y is single prec. */ | |
207 | if (__glibc_unlikely ((asuint64 (ylogx) >> 47 & 0xffff) | |
208 | >= asuint64 (126.0 * POWF_SCALE) >> 47)) | |
209 | { | |
210 | /* |y*log(x)| >= 126. */ | |
211 | if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE) | |
505b5b29 | 212 | /* |x^y| > 0x1.ffffffp127. */ |
4ea49f4c | 213 | return __math_oflowf (sign_bias); |
505b5b29 SN |
214 | if (WANT_ROUNDING && WANT_ERRNO |
215 | && ylogx > 0x1.fffffffa3aae2p+6 * POWF_SCALE) | |
216 | /* |x^y| > 0x1.fffffep127, check if we round away from 0. */ | |
217 | if ((!sign_bias | |
218 | && math_narrow_eval (1.0f + math_opt_barrier (0x1p-25f)) != 1.0f) | |
219 | || (sign_bias | |
220 | && math_narrow_eval (-1.0f - math_opt_barrier (0x1p-25f)) | |
221 | != -1.0f)) | |
222 | return __math_oflowf (sign_bias); | |
4ea49f4c SN |
223 | if (ylogx <= -150.0 * POWF_SCALE) |
224 | return __math_uflowf (sign_bias); | |
225 | #if WANT_ERRNO_UFLOW | |
226 | if (ylogx < -149.0 * POWF_SCALE) | |
227 | return __math_may_uflowf (sign_bias); | |
228 | #endif | |
229 | } | |
230 | return (float) exp2_inline (ylogx, sign_bias); | |
f7eac6eb | 231 | } |
bd4430c2 SN |
232 | #ifndef __powf |
233 | strong_alias (__powf, __ieee754_powf) | |
234 | strong_alias (__powf, __powf_finite) | |
235 | versioned_symbol (libm, __powf, powf, GLIBC_2_27); | |
24b6515d | 236 | libm_alias_float_other (__pow, pow) |
bd4430c2 | 237 | #endif |