1 /* Double-precision x^y function.
2 Copyright (c) 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_DOUBLE
31
32 #include <math.h>
33 #include <stdint.h>
34 #include "math_config.h"
35
36 /*
37 Worst-case error: 0.54 ULP (~= ulperr_exp + 1024*Ln2*relerr_log*2^53)
38 relerr_log: 1.3 * 2^-68 (Relative error of log, 1.5 * 2^-68 without fma)
39 ulperr_exp: 0.509 ULP (ULP error of exp, 0.511 ULP without fma)
40 */
41
42 #define T __pow_log_data.tab
43 #define A __pow_log_data.poly
44 #define Ln2hi __pow_log_data.ln2hi
45 #define Ln2lo __pow_log_data.ln2lo
46 #define N (1 << POW_LOG_TABLE_BITS)
47 #define OFF 0x3fe6955500000000
48
49 /* Top 12 bits of a double (sign and exponent bits). */
50 static inline uint32_t
top12(double x)51 top12 (double x)
52 {
53 return asuint64 (x) >> 52;
54 }
55
56 /* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
57 additional 15 bits precision. IX is the bit representation of x, but
58 normalized in the subnormal range using the sign bit for the exponent. */
59 static inline double_t
log_inline(uint64_t ix,double_t * tail)60 log_inline (uint64_t ix, double_t *tail)
61 {
62 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
63 double_t z, r, y, invc, logc, logctail, kd, hi, t1, t2, lo, lo1, lo2, p;
64 uint64_t iz, tmp;
65 int k, i;
66
67 /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
68 The range is split into N subintervals.
69 The ith subinterval contains z and c is near its center. */
70 tmp = ix - OFF;
71 i = (tmp >> (52 - POW_LOG_TABLE_BITS)) % N;
72 k = (int64_t) tmp >> 52; /* arithmetic shift */
73 iz = ix - (tmp & 0xfffULL << 52);
74 z = asfloat64 (iz);
75 kd = (double_t) k;
76
77 /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */
78 invc = T[i].invc;
79 logc = T[i].logc;
80 logctail = T[i].logctail;
81
82 /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
83 |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
84 #if _HAVE_FAST_FMA
85 r = fma (z, invc, -1.0);
86 #else
87 /* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|. */
88 double_t zhi = asfloat64 ((iz + (1ULL << 31)) & (-1ULL << 32));
89 double_t zlo = z - zhi;
90 double_t rhi = zhi * invc - 1.0;
91 double_t rlo = zlo * invc;
92 r = rhi + rlo;
93 #endif
94
95 /* k*Ln2 + log(c) + r. */
96 t1 = kd * Ln2hi + logc;
97 t2 = t1 + r;
98 lo1 = kd * Ln2lo + logctail;
99 lo2 = t1 - t2 + r;
100
101 /* Evaluation is optimized assuming superscalar pipelined execution. */
102 double_t ar, ar2, ar3, lo3, lo4;
103 ar = A[0] * r; /* A[0] = -0.5. */
104 ar2 = r * ar;
105 ar3 = r * ar2;
106 /* k*Ln2 + log(c) + r + A[0]*r*r. */
107 #if _HAVE_FAST_FMA
108 hi = t2 + ar2;
109 lo3 = fma (ar, r, -ar2);
110 lo4 = t2 - hi + ar2;
111 #else
112 double_t arhi = A[0] * rhi;
113 double_t arhi2 = rhi * arhi;
114 hi = t2 + arhi2;
115 lo3 = rlo * (ar + arhi);
116 lo4 = t2 - hi + arhi2;
117 #endif
118 /* p = log1p(r) - r - A[0]*r*r. */
119 #if POW_LOG_POLY_ORDER == 8
120 p = (ar3
121 * (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6]))));
122 #endif
123 lo = lo1 + lo2 + lo3 + lo4 + p;
124 y = hi + lo;
125 *tail = hi - y + lo;
126 return y;
127 }
128
129 #undef N
130 #undef T
131 #define N (1 << EXP_TABLE_BITS)
132 #define InvLn2N __exp_data.invln2N
133 #define NegLn2hiN __exp_data.negln2hiN
134 #define NegLn2loN __exp_data.negln2loN
135 #define Shift __exp_data.shift
136 #define T __exp_data.tab
137 #define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
138 #define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
139 #define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
140 #define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
141 #define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
142
143 /* Handle cases that may overflow or underflow when computing the result that
144 is scale*(1+TMP) without intermediate rounding. The bit representation of
145 scale is in SBITS, however it has a computed exponent that may have
146 overflown into the sign bit so that needs to be adjusted before using it as
147 a double. (int32_t)KI is the k used in the argument reduction and exponent
148 adjustment of scale, positive k here means the result may overflow and
149 negative k means the result may underflow. */
150 static inline double
specialcase(double_t tmp,uint64_t sbits,uint64_t ki)151 specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
152 {
153 double_t scale, y;
154
155 if ((ki & 0x80000000) == 0)
156 {
157 /* k > 0, the exponent of scale might have overflowed by <= 460. */
158 sbits -= 1009ull << 52;
159 scale = asfloat64 (sbits);
160 y = 0x1p1009 * (scale + scale * tmp);
161 return check_oflow (y);
162 }
163 /* k < 0, need special care in the subnormal range. */
164 sbits += 1022ull << 52;
165 /* Note: sbits is signed scale. */
166 scale = asfloat64 (sbits);
167 y = scale + scale * tmp;
168 #if FLT_EVAL_METHOD == 2
169 #define fabs(x) fabsl(x)
170 #endif
171 if (fabs (y) < 1.0)
172 #undef fabs
173 {
174 /* Round y to the right precision before scaling it into the subnormal
175 range to avoid double rounding that can cause 0.5+E/2 ulp error where
176 E is the worst-case ulp error outside the subnormal range. So this
177 is only useful if the goal is better than 1 ulp worst-case error. */
178 double_t hi, lo, one = 1.0;
179 if (y < 0.0)
180 one = -1.0;
181 lo = scale - y + scale * tmp;
182 hi = one + y;
183 lo = one - hi + y + lo;
184 y = eval_as_double (hi + lo) - one;
185 /* Fix the sign of 0. */
186 if (y == 0.0)
187 y = asfloat64 (sbits & 0x8000000000000000);
188 /* The underflow exception needs to be signaled explicitly. */
189 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
190 }
191 y = 0x1p-1022 * y;
192 return check_uflow (y);
193 }
194
195 #define SIGN_BIAS ((int32_t) 0x800 << EXP_TABLE_BITS)
196
197 /* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
198 The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */
199 static inline double
exp_inline(double x,double xtail,uint32_t sign_bias)200 exp_inline (double x, double xtail, uint32_t sign_bias)
201 {
202 uint32_t abstop;
203 uint64_t ki, idx, top, sbits;
204 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
205 double_t kd, z, r, r2, scale, tail, tmp;
206
207 abstop = top12 (x) & 0x7ff;
208 if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
209 {
210 if (abstop - top12 (0x1p-54) >= 0x80000000)
211 {
212 /* Avoid spurious underflow for tiny x. */
213 /* Note: 0 is common input. */
214 double_t one = WANT_ROUNDING ? 1.0 + x : 1.0;
215 return sign_bias ? -one : one;
216 }
217 if (abstop >= top12 (1024.0))
218 {
219 /* Note: inf and nan are already handled. */
220 if (asuint64 (x) >> 63)
221 return __math_uflow (sign_bias);
222 else
223 return __math_oflow (sign_bias);
224 }
225 /* Large x is special cased below. */
226 abstop = 0;
227 }
228
229 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
230 /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */
231 z = InvLn2N * x;
232 #if TOINT_INTRINSICS
233 kd = roundtoint (z);
234 ki = converttoint (z);
235 #elif EXP_USE_TOINT_NARROW
236 /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes. */
237 kd = eval_as_double (z + Shift);
238 ki = asuint64 (kd) >> 16;
239 kd = (double_t) (int32_t) ki;
240 #else
241 /* z - kd is in [-1, 1] in non-nearest rounding modes. */
242 kd = eval_as_double (z + Shift);
243 ki = asuint64 (kd);
244 kd -= Shift;
245 #endif
246 r = x + kd * NegLn2hiN + kd * NegLn2loN;
247 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
248 r += xtail;
249 /* 2^(k/N) ~= scale * (1 + tail). */
250 idx = 2 * (ki % N);
251 top = (ki + sign_bias) << (52 - EXP_TABLE_BITS);
252 tail = asfloat64 (T[idx]);
253 /* This is only a valid scale when -1023*N < k < 1024*N. */
254 sbits = T[idx + 1] + top;
255 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */
256 /* Evaluation is optimized assuming superscalar pipelined execution. */
257 r2 = r * r;
258 /* Without fma the worst case error is 0.25/N ulp larger. */
259 /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */
260 #if EXP_POLY_ORDER == 4
261 tmp = tail + r + r2 * C2 + r * r2 * (C3 + r * C4);
262 #elif EXP_POLY_ORDER == 5
263 tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
264 #elif EXP_POLY_ORDER == 6
265 tmp = tail + r + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6);
266 #endif
267 if (unlikely (abstop == 0))
268 return specialcase (tmp, sbits, ki);
269 scale = asfloat64 (sbits);
270 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
271 is no spurious underflow here even without fma. */
272 return scale + scale * tmp;
273 }
274
275 /* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
276 the bit representation of a non-zero finite floating-point value. */
277 static inline int
checkint(uint64_t iy)278 checkint (uint64_t iy)
279 {
280 int e = iy >> 52 & 0x7ff;
281 if (e < 0x3ff)
282 return 0;
283 if (e > 0x3ff + 52)
284 return 2;
285 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
286 return 0;
287 if (iy & (1ULL << (0x3ff + 52 - e)))
288 return 1;
289 return 2;
290 }
291
292 /* Returns 1 if input is the bit representation of 0, infinity or nan. */
293 static inline int
zeroinfnan(uint64_t i)294 zeroinfnan (uint64_t i)
295 {
296 return 2 * i - 1 >= 2 * asuint64 ((double) INFINITY) - 1;
297 }
298
299 double
pow(double x,double y)300 pow (double x, double y)
301 {
302 uint32_t sign_bias = 0;
303 uint64_t ix, iy;
304 uint32_t topx, topy;
305
306 ix = asuint64 (x);
307 iy = asuint64 (y);
308 topx = top12 (x);
309 topy = top12 (y);
310 if (unlikely (topx - 0x001 >= 0x7ff - 0x001
311 || (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be))
312 {
313 /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
314 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */
315 /* Special cases: (x < 0x1p-126 or inf or nan) or
316 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */
317 if (unlikely (zeroinfnan (iy)))
318 {
319 if (2 * iy == 0)
320 return issignaling64_inline (x) ? x + y : 1.0;
321 if (ix == asuint64 (1.0))
322 return issignaling64_inline (y) ? x + y : 1.0;
323 if (2 * ix > 2 * asuint64 ((double) INFINITY)
324 || 2 * iy > 2 * asuint64 ((double) INFINITY))
325 return x + y;
326 if (2 * ix == 2 * asuint64 (1.0))
327 return 1.0;
328 if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
329 return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
330 return y * y;
331 }
332 if (unlikely (zeroinfnan (ix)))
333 {
334 double_t x2 = x * x;
335 if (ix >> 63 && checkint (iy) == 1)
336 {
337 x2 = -x2;
338 sign_bias = 1;
339 }
340 if (WANT_ERRNO && 2 * ix == 0 && iy >> 63)
341 return __math_divzero (sign_bias);
342 /* Without the barrier some versions of clang hoist the 1/x2 and
343 thus division by zero exception can be signaled spuriously. */
344 return iy >> 63 ? opt_barrier_double (1 / x2) : x2;
345 }
346 /* Here x and y are non-zero finite. */
347 if (ix >> 63)
348 {
349 /* Finite x < 0. */
350 int yint = checkint (iy);
351 if (yint == 0)
352 return __math_invalid (x);
353 if (yint == 1)
354 sign_bias = SIGN_BIAS;
355 ix &= 0x7fffffffffffffff;
356 topx &= 0x7ff;
357 }
358 if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be)
359 {
360 /* Note: sign_bias == 0 here because y is not odd. */
361 if (ix == asuint64 (1.0))
362 return 1.0;
363 if ((topy & 0x7ff) < 0x3be)
364 {
365 /* |y| < 2^-65, x^y ~= 1 + y*log(x). */
366 if (WANT_ROUNDING)
367 return ix > asuint64 (1.0) ? 1.0 + y : 1.0 - y;
368 else
369 return 1.0;
370 }
371 return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0)
372 : __math_uflow (0);
373 }
374 if (topx == 0)
375 {
376 /* Normalize subnormal x so exponent becomes negative. */
377 ix = asuint64 (x * 0x1p52);
378 ix &= 0x7fffffffffffffff;
379 ix -= 52ULL << 52;
380 }
381 }
382
383 double_t lo;
384 double_t hi = log_inline (ix, &lo);
385 double_t ehi, elo;
386 #if _HAVE_FAST_FMA
387 ehi = y * hi;
388 elo = y * lo + fma (y, hi, -ehi);
389 #else
390 double_t yhi = asfloat64 (iy & -1ULL << 27);
391 double_t ylo = y - yhi;
392 double_t lhi = asfloat64 (asuint64 (hi) & -1ULL << 27);
393 double_t llo = hi - lhi + lo;
394 ehi = yhi * lhi;
395 elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25. */
396 #endif
397 return exp_inline (ehi, elo, sign_bias);
398 }
399
400 #if defined(_HAVE_ALIAS_ATTRIBUTE)
401 #ifndef __clang__
402 #pragma GCC diagnostic ignored "-Wmissing-attributes"
403 #endif
404 __strong_reference(pow, _pow);
405 #endif
406
407 _MATH_ALIAS_d_dd(pow)
408
409 #endif
410
