1 /*
2 * Double-precision x^y function.
3 *
4 * Copyright (c) 2018-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "math_config.h"
9
10 /* Scalar version of pow used for fallbacks in vector implementations. */
11
12 /* Data is defined in v_pow_log_data.c. */
13 #define N_LOG (1 << V_POW_LOG_TABLE_BITS)
14 #define Off 0x3fe6955500000000
15 #define As __v_pow_log_data.poly
16
17 /* Data is defined in v_pow_exp_data.c. */
18 #define N_EXP (1 << V_POW_EXP_TABLE_BITS)
19 #define SignBias (0x800 << V_POW_EXP_TABLE_BITS)
20 #define SmallExp 0x3c9 /* top12(0x1p-54). */
21 #define BigExp 0x408 /* top12(512.0). */
22 #define ThresExp 0x03f /* BigExp - SmallExp. */
23 #define InvLn2N __v_pow_exp_data.n_over_ln2
24 #define Ln2HiN __v_pow_exp_data.ln2_over_n_hi
25 #define Ln2LoN __v_pow_exp_data.ln2_over_n_lo
26 #define SBits __v_pow_exp_data.sbits
27 #define Cs __v_pow_exp_data.poly
28
29 /* Constants associated with pow. */
30 #define SmallPowX 0x001 /* top12(0x1p-126). */
31 #define BigPowX 0x7ff /* top12(INFINITY). */
32 #define ThresPowX 0x7fe /* BigPowX - SmallPowX. */
33 #define SmallPowY 0x3be /* top12(0x1.e7b6p-65). */
34 #define BigPowY 0x43e /* top12(0x1.749p62). */
35 #define ThresPowY 0x080 /* BigPowY - SmallPowY. */
36
37 /* Top 12 bits of a double (sign and exponent bits). */
38 static inline uint32_t
top12(double x)39 top12 (double x)
40 {
41 return asuint64 (x) >> 52;
42 }
43
44 /* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
45 additional 15 bits precision. IX is the bit representation of x, but
46 normalized in the subnormal range using the sign bit for the exponent. */
47 static inline double
log_inline(uint64_t ix,double * tail)48 log_inline (uint64_t ix, double *tail)
49 {
50 /* x = 2^k z; where z is in range [Off,2*Off) and exact.
51 The range is split into N subintervals.
52 The ith subinterval contains z and c is near its center. */
53 uint64_t tmp = ix - Off;
54 int i = (tmp >> (52 - V_POW_LOG_TABLE_BITS)) & (N_LOG - 1);
55 int k = (int64_t) tmp >> 52; /* arithmetic shift. */
56 uint64_t iz = ix - (tmp & 0xfffULL << 52);
57 double z = asdouble (iz);
58 double kd = (double) k;
59
60 /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */
61 double invc = __v_pow_log_data.invc[i];
62 double logc = __v_pow_log_data.logc[i];
63 double logctail = __v_pow_log_data.logctail[i];
64
65 /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
66 |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
67 double r = fma (z, invc, -1.0);
68
69 /* k*Ln2 + log(c) + r. */
70 double t1 = kd * __v_pow_log_data.ln2_hi + logc;
71 double t2 = t1 + r;
72 double lo1 = kd * __v_pow_log_data.ln2_lo + logctail;
73 double lo2 = t1 - t2 + r;
74
75 /* Evaluation is optimized assuming superscalar pipelined execution. */
76 double ar = As[0] * r;
77 double ar2 = r * ar;
78 double ar3 = r * ar2;
79 /* k*Ln2 + log(c) + r + A[0]*r*r. */
80 double hi = t2 + ar2;
81 double lo3 = fma (ar, r, -ar2);
82 double lo4 = t2 - hi + ar2;
83 /* p = log1p(r) - r - A[0]*r*r. */
84 double p = (ar3
85 * (As[1] + r * As[2]
86 + ar2 * (As[3] + r * As[4] + ar2 * (As[5] + r * As[6]))));
87 double lo = lo1 + lo2 + lo3 + lo4 + p;
88 double y = hi + lo;
89 *tail = hi - y + lo;
90 return y;
91 }
92
93 /* Handle cases that may overflow or underflow when computing the result that
94 is scale*(1+TMP) without intermediate rounding. The bit representation of
95 scale is in SBITS, however it has a computed exponent that may have
96 overflown into the sign bit so that needs to be adjusted before using it as
97 a double. (int32_t)KI is the k used in the argument reduction and exponent
98 adjustment of scale, positive k here means the result may overflow and
99 negative k means the result may underflow. */
100 static inline double
special_case(double tmp,uint64_t sbits,uint64_t ki)101 special_case (double tmp, uint64_t sbits, uint64_t ki)
102 {
103 double scale, y;
104
105 if ((ki & 0x80000000) == 0)
106 {
107 /* k > 0, the exponent of scale might have overflowed by <= 460. */
108 sbits -= 1009ull << 52;
109 scale = asdouble (sbits);
110 y = 0x1p1009 * (scale + scale * tmp);
111 return check_oflow (eval_as_double (y));
112 }
113 /* k < 0, need special care in the subnormal range. */
114 sbits += 1022ull << 52;
115 /* Note: sbits is signed scale. */
116 scale = asdouble (sbits);
117 y = scale + scale * tmp;
118 #if WANT_SIMD_EXCEPT
119 if (fabs (y) < 1.0)
120 {
121 /* Round y to the right precision before scaling it into the subnormal
122 range to avoid double rounding that can cause 0.5+E/2 ulp error where
123 E is the worst-case ulp error outside the subnormal range. So this
124 is only useful if the goal is better than 1 ulp worst-case error. */
125 double hi, lo, one = 1.0;
126 if (y < 0.0)
127 one = -1.0;
128 lo = scale - y + scale * tmp;
129 hi = one + y;
130 lo = one - hi + y + lo;
131 y = eval_as_double (hi + lo) - one;
132 /* Fix the sign of 0. */
133 if (y == 0.0)
134 y = asdouble (sbits & 0x8000000000000000);
135 /* The underflow exception needs to be signaled explicitly. */
136 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
137 }
138 #endif
139 y = 0x1p-1022 * y;
140 return check_uflow (eval_as_double (y));
141 }
142
143 /* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
144 The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1. */
145 static inline double
exp_inline(double x,double xtail,uint32_t sign_bias)146 exp_inline (double x, double xtail, uint32_t sign_bias)
147 {
148 uint32_t abstop = top12 (x) & 0x7ff;
149 if (unlikely (abstop - SmallExp >= ThresExp))
150 {
151 if (abstop - SmallExp >= 0x80000000)
152 {
153 /* Avoid spurious underflow for tiny x. */
154 /* Note: 0 is common input. */
155 return sign_bias ? -1.0 : 1.0;
156 }
157 if (abstop >= top12 (1024.0))
158 {
159 /* Note: inf and nan are already handled. */
160 /* Skip errno handling. */
161 #if WANT_SIMD_EXCEPT
162 return asuint64 (x) >> 63 ? __math_uflow (sign_bias)
163 : __math_oflow (sign_bias);
164 #else
165 double res_uoflow = asuint64 (x) >> 63 ? 0.0 : INFINITY;
166 return sign_bias ? -res_uoflow : res_uoflow;
167 #endif
168 }
169 /* Large x is special cased below. */
170 abstop = 0;
171 }
172
173 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
174 /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */
175 double z = InvLn2N * x;
176 double kd = round (z);
177 uint64_t ki = lround (z);
178 double r = x - kd * Ln2HiN - kd * Ln2LoN;
179 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
180 r += xtail;
181 /* 2^(k/N) ~= scale. */
182 uint64_t idx = ki & (N_EXP - 1);
183 uint64_t top = (ki + sign_bias) << (52 - V_POW_EXP_TABLE_BITS);
184 /* This is only a valid scale when -1023*N < k < 1024*N. */
185 uint64_t sbits = SBits[idx] + top;
186 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */
187 /* Evaluation is optimized assuming superscalar pipelined execution. */
188 double r2 = r * r;
189 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
190 if (unlikely (abstop == 0))
191 return special_case (tmp, sbits, ki);
192 double scale = asdouble (sbits);
193 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
194 is no spurious underflow here even without fma. */
195 return eval_as_double (scale + scale * tmp);
196 }
197
198 /* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
199 A version of exp_inline that is not inlined and for which sign_bias is
200 equal to 0. */
201 static double NOINLINE
exp_nosignbias(double x,double xtail)202 exp_nosignbias (double x, double xtail)
203 {
204 uint32_t abstop = top12 (x) & 0x7ff;
205 if (unlikely (abstop - SmallExp >= ThresExp))
206 {
207 /* Avoid spurious underflow for tiny x. */
208 if (abstop - SmallExp >= 0x80000000)
209 return 1.0;
210 /* Note: inf and nan are already handled. */
211 if (abstop >= top12 (1024.0))
212 #if WANT_SIMD_EXCEPT
213 return asuint64 (x) >> 63 ? __math_uflow (0) : __math_oflow (0);
214 #else
215 return asuint64 (x) >> 63 ? 0.0 : INFINITY;
216 #endif
217 /* Large x is special cased below. */
218 abstop = 0;
219 }
220
221 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
222 /* x = ln2/N*k + r, with k integer and r in [-ln2/2N, ln2/2N]. */
223 double z = InvLn2N * x;
224 double kd = round (z);
225 uint64_t ki = lround (z);
226 double r = x - kd * Ln2HiN - kd * Ln2LoN;
227 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
228 r += xtail;
229 /* 2^(k/N) ~= scale. */
230 uint64_t idx = ki & (N_EXP - 1);
231 uint64_t top = ki << (52 - V_POW_EXP_TABLE_BITS);
232 /* This is only a valid scale when -1023*N < k < 1024*N. */
233 uint64_t sbits = SBits[idx] + top;
234 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */
235 double r2 = r * r;
236 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
237 if (unlikely (abstop == 0))
238 return special_case (tmp, sbits, ki);
239 double scale = asdouble (sbits);
240 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
241 is no spurious underflow here even without fma. */
242 return eval_as_double (scale + scale * tmp);
243 }
244
245 /* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
246 the bit representation of a non-zero finite floating-point value. */
247 static inline int
checkint(uint64_t iy)248 checkint (uint64_t iy)
249 {
250 int e = iy >> 52 & 0x7ff;
251 if (e < 0x3ff)
252 return 0;
253 if (e > 0x3ff + 52)
254 return 2;
255 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
256 return 0;
257 if (iy & (1ULL << (0x3ff + 52 - e)))
258 return 1;
259 return 2;
260 }
261
262 /* Returns 1 if input is the bit representation of 0, infinity or nan. */
263 static inline int
zeroinfnan(uint64_t i)264 zeroinfnan (uint64_t i)
265 {
266 return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;
267 }
268
269 static double NOINLINE
__pl_finite_pow(double x,double y)270 __pl_finite_pow (double x, double y)
271 {
272 uint32_t sign_bias = 0;
273 uint64_t ix, iy;
274 uint32_t topx, topy;
275
276 ix = asuint64 (x);
277 iy = asuint64 (y);
278 topx = top12 (x);
279 topy = top12 (y);
280 if (unlikely (topx - SmallPowX >= ThresPowX
281 || (topy & 0x7ff) - SmallPowY >= ThresPowY))
282 {
283 /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
284 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */
285 /* Special cases: (x < 0x1p-126 or inf or nan) or
286 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */
287 if (unlikely (zeroinfnan (iy)))
288 {
289 if (2 * iy == 0)
290 return issignaling_inline (x) ? x + y : 1.0;
291 if (ix == asuint64 (1.0))
292 return issignaling_inline (y) ? x + y : 1.0;
293 if (2 * ix > 2 * asuint64 (INFINITY)
294 || 2 * iy > 2 * asuint64 (INFINITY))
295 return x + y;
296 if (2 * ix == 2 * asuint64 (1.0))
297 return 1.0;
298 if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
299 return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
300 return y * y;
301 }
302 if (unlikely (zeroinfnan (ix)))
303 {
304 double x2 = x * x;
305 if (ix >> 63 && checkint (iy) == 1)
306 {
307 x2 = -x2;
308 sign_bias = 1;
309 }
310 #if WANT_SIMD_EXCEPT
311 if (2 * ix == 0 && iy >> 63)
312 return __math_divzero (sign_bias);
313 #endif
314 /* Without the barrier some versions of clang hoist the 1/x2 and
315 thus division by zero exception can be signaled spuriously. */
316 return iy >> 63 ? opt_barrier_double (1 / x2) : x2;
317 }
318 /* Here x and y are non-zero finite. */
319 if (ix >> 63)
320 {
321 /* Finite x < 0. */
322 int yint = checkint (iy);
323 if (yint == 0)
324 #if WANT_SIMD_EXCEPT
325 return __math_invalid (x);
326 #else
327 return __builtin_nan ("");
328 #endif
329 if (yint == 1)
330 sign_bias = SignBias;
331 ix &= 0x7fffffffffffffff;
332 topx &= 0x7ff;
333 }
334 if ((topy & 0x7ff) - SmallPowY >= ThresPowY)
335 {
336 /* Note: sign_bias == 0 here because y is not odd. */
337 if (ix == asuint64 (1.0))
338 return 1.0;
339 /* |y| < 2^-65, x^y ~= 1 + y*log(x). */
340 if ((topy & 0x7ff) < SmallPowY)
341 return 1.0;
342 #if WANT_SIMD_EXCEPT
343 return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0)
344 : __math_uflow (0);
345 #else
346 return (ix > asuint64 (1.0)) == (topy < 0x800) ? INFINITY : 0;
347 #endif
348 }
349 if (topx == 0)
350 {
351 /* Normalize subnormal x so exponent becomes negative. */
352 /* Without the barrier some versions of clang evalutate the mul
353 unconditionally causing spurious overflow exceptions. */
354 ix = asuint64 (opt_barrier_double (x) * 0x1p52);
355 ix &= 0x7fffffffffffffff;
356 ix -= 52ULL << 52;
357 }
358 }
359
360 double lo;
361 double hi = log_inline (ix, &lo);
362 double ehi = y * hi;
363 double elo = y * lo + fma (y, hi, -ehi);
364 return exp_inline (ehi, elo, sign_bias);
365 }
366