glibc/sysdeps/aarch64/fpu/finite_pow.h

374 lines
12 KiB
C
Raw Normal View History

/* Double-precision x^y function.
Copyright (C) 2024 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include "math_config.h"
/* Scalar version of pow used for fallbacks in vector implementations. */
/* Data is defined in v_pow_log_data.c. */
#define N_LOG (1 << V_POW_LOG_TABLE_BITS)
#define Off 0x3fe6955500000000
#define As __v_pow_log_data.poly
/* Data is defined in v_pow_exp_data.c. */
#define N_EXP (1 << V_POW_EXP_TABLE_BITS)
#define SignBias (0x800 << V_POW_EXP_TABLE_BITS)
#define SmallExp 0x3c9 /* top12(0x1p-54). */
#define BigExp 0x408 /* top12(512.0). */
#define ThresExp 0x03f /* BigExp - SmallExp. */
#define InvLn2N __v_pow_exp_data.n_over_ln2
#define Ln2HiN __v_pow_exp_data.ln2_over_n_hi
#define Ln2LoN __v_pow_exp_data.ln2_over_n_lo
#define SBits __v_pow_exp_data.sbits
#define Cs __v_pow_exp_data.poly
/* Constants associated with pow. */
#define SmallPowX 0x001 /* top12(0x1p-126). */
#define BigPowX 0x7ff /* top12(INFINITY). */
#define ThresPowX 0x7fe /* BigPowX - SmallPowX. */
#define SmallPowY 0x3be /* top12(0x1.e7b6p-65). */
#define BigPowY 0x43e /* top12(0x1.749p62). */
#define ThresPowY 0x080 /* BigPowY - SmallPowY. */
/* Top 12 bits of a double (sign and exponent bits). */
static inline uint32_t
top12 (double x)
{
return asuint64 (x) >> 52;
}
/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
additional 15 bits precision. IX is the bit representation of x, but
normalized in the subnormal range using the sign bit for the exponent. */
static inline double
log_inline (uint64_t ix, double *tail)
{
/* x = 2^k z; where z is in range [Off,2*Off) and exact.
The range is split into N subintervals.
The ith subinterval contains z and c is near its center. */
uint64_t tmp = ix - Off;
int i = (tmp >> (52 - V_POW_LOG_TABLE_BITS)) & (N_LOG - 1);
int k = (int64_t) tmp >> 52; /* arithmetic shift. */
uint64_t iz = ix - (tmp & 0xfffULL << 52);
double z = asdouble (iz);
double kd = (double) k;
/* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */
double invc = __v_pow_log_data.invc[i];
double logc = __v_pow_log_data.logc[i];
double logctail = __v_pow_log_data.logctail[i];
/* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
|z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
double r = fma (z, invc, -1.0);
/* k*Ln2 + log(c) + r. */
double t1 = kd * __v_pow_log_data.ln2_hi + logc;
double t2 = t1 + r;
double lo1 = kd * __v_pow_log_data.ln2_lo + logctail;
double lo2 = t1 - t2 + r;
/* Evaluation is optimized assuming superscalar pipelined execution. */
double ar = As[0] * r;
double ar2 = r * ar;
double ar3 = r * ar2;
/* k*Ln2 + log(c) + r + A[0]*r*r. */
double hi = t2 + ar2;
double lo3 = fma (ar, r, -ar2);
double lo4 = t2 - hi + ar2;
/* p = log1p(r) - r - A[0]*r*r. */
double p = (ar3
* (As[1] + r * As[2]
+ ar2 * (As[3] + r * As[4] + ar2 * (As[5] + r * As[6]))));
double lo = lo1 + lo2 + lo3 + lo4 + p;
double y = hi + lo;
*tail = hi - y + lo;
return y;
}
/* Handle cases that may overflow or underflow when computing the result that
is scale*(1+TMP) without intermediate rounding. The bit representation of
scale is in SBITS, however it has a computed exponent that may have
overflown into the sign bit so that needs to be adjusted before using it as
a double. (int32_t)KI is the k used in the argument reduction and exponent
adjustment of scale, positive k here means the result may overflow and
negative k means the result may underflow. */
static inline double
special_case (double tmp, uint64_t sbits, uint64_t ki)
{
double scale, y;
if ((ki & 0x80000000) == 0)
{
/* k > 0, the exponent of scale might have overflowed by <= 460. */
sbits -= 1009ull << 52;
scale = asdouble (sbits);
y = 0x1p1009 * (scale + scale * tmp);
return y;
}
/* k < 0, need special care in the subnormal range. */
sbits += 1022ull << 52;
/* Note: sbits is signed scale. */
scale = asdouble (sbits);
y = scale + scale * tmp;
#if WANT_SIMD_EXCEPT
if (fabs (y) < 1.0)
{
/* Round y to the right precision before scaling it into the subnormal
range to avoid double rounding that can cause 0.5+E/2 ulp error where
E is the worst-case ulp error outside the subnormal range. So this
is only useful if the goal is better than 1 ulp worst-case error. */
double hi, lo, one = 1.0;
if (y < 0.0)
one = -1.0;
lo = scale - y + scale * tmp;
hi = one + y;
lo = one - hi + y + lo;
y = (hi + lo) - one;
/* Fix the sign of 0. */
if (y == 0.0)
y = asdouble (sbits & 0x8000000000000000);
/* The underflow exception needs to be signaled explicitly. */
force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
}
#endif
y = 0x1p-1022 * y;
return y;
}
/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1. */
static inline double
exp_inline (double x, double xtail, uint32_t sign_bias)
{
uint32_t abstop = top12 (x) & 0x7ff;
if (__glibc_unlikely (abstop - SmallExp >= ThresExp))
{
if (abstop - SmallExp >= 0x80000000)
{
/* Avoid spurious underflow for tiny x. */
/* Note: 0 is common input. */
return sign_bias ? -1.0 : 1.0;
}
if (abstop >= top12 (1024.0))
{
/* Note: inf and nan are already handled. */
/* Skip errno handling. */
#if WANT_SIMD_EXCEPT
return asuint64 (x) >> 63 ? __math_uflow (sign_bias)
: __math_oflow (sign_bias);
#else
double res_uoflow = asuint64 (x) >> 63 ? 0.0 : INFINITY;
return sign_bias ? -res_uoflow : res_uoflow;
#endif
}
/* Large x is special cased below. */
abstop = 0;
}
/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
/* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */
double z = InvLn2N * x;
double kd = round (z);
uint64_t ki = lround (z);
double r = x - kd * Ln2HiN - kd * Ln2LoN;
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
r += xtail;
/* 2^(k/N) ~= scale. */
uint64_t idx = ki & (N_EXP - 1);
uint64_t top = (ki + sign_bias) << (52 - V_POW_EXP_TABLE_BITS);
/* This is only a valid scale when -1023*N < k < 1024*N. */
uint64_t sbits = SBits[idx] + top;
/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */
/* Evaluation is optimized assuming superscalar pipelined execution. */
double r2 = r * r;
double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
if (__glibc_unlikely (abstop == 0))
return special_case (tmp, sbits, ki);
double scale = asdouble (sbits);
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
is no spurious underflow here even without fma. */
return scale + scale * tmp;
}
/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
A version of exp_inline that is not inlined and for which sign_bias is
equal to 0. */
static double NOINLINE
exp_nosignbias (double x, double xtail)
{
uint32_t abstop = top12 (x) & 0x7ff;
if (__glibc_unlikely (abstop - SmallExp >= ThresExp))
{
/* Avoid spurious underflow for tiny x. */
if (abstop - SmallExp >= 0x80000000)
return 1.0;
/* Note: inf and nan are already handled. */
if (abstop >= top12 (1024.0))
#if WANT_SIMD_EXCEPT
return asuint64 (x) >> 63 ? __math_uflow (0) : __math_oflow (0);
#else
return asuint64 (x) >> 63 ? 0.0 : INFINITY;
#endif
/* Large x is special cased below. */
abstop = 0;
}
/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
/* x = ln2/N*k + r, with k integer and r in [-ln2/2N, ln2/2N]. */
double z = InvLn2N * x;
double kd = round (z);
uint64_t ki = lround (z);
double r = x - kd * Ln2HiN - kd * Ln2LoN;
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
r += xtail;
/* 2^(k/N) ~= scale. */
uint64_t idx = ki & (N_EXP - 1);
uint64_t top = ki << (52 - V_POW_EXP_TABLE_BITS);
/* This is only a valid scale when -1023*N < k < 1024*N. */
uint64_t sbits = SBits[idx] + top;
/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */
double r2 = r * r;
double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
if (__glibc_unlikely (abstop == 0))
return special_case (tmp, sbits, ki);
double scale = asdouble (sbits);
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
is no spurious underflow here even without fma. */
return scale + scale * tmp;
}
/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
the bit representation of a non-zero finite floating-point value. */
static inline int
checkint (uint64_t iy)
{
int e = iy >> 52 & 0x7ff;
if (e < 0x3ff)
return 0;
if (e > 0x3ff + 52)
return 2;
if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
return 0;
if (iy & (1ULL << (0x3ff + 52 - e)))
return 1;
return 2;
}
/* Returns 1 if input is the bit representation of 0, infinity or nan. */
static inline int
zeroinfnan (uint64_t i)
{
return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;
}
static double NOINLINE
pow_scalar_special_case (double x, double y)
{
uint32_t sign_bias = 0;
uint64_t ix, iy;
uint32_t topx, topy;
ix = asuint64 (x);
iy = asuint64 (y);
topx = top12 (x);
topy = top12 (y);
if (__glibc_unlikely (topx - SmallPowX >= ThresPowX
|| (topy & 0x7ff) - SmallPowY >= ThresPowY))
{
/* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */
/* Special cases: (x < 0x1p-126 or inf or nan) or
(|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */
if (__glibc_unlikely (zeroinfnan (iy)))
{
if (2 * iy == 0)
return issignaling_inline (x) ? x + y : 1.0;
if (ix == asuint64 (1.0))
return issignaling_inline (y) ? x + y : 1.0;
if (2 * ix > 2 * asuint64 (INFINITY)
|| 2 * iy > 2 * asuint64 (INFINITY))
return x + y;
if (2 * ix == 2 * asuint64 (1.0))
return 1.0;
if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
return y * y;
}
if (__glibc_unlikely (zeroinfnan (ix)))
{
double x2 = x * x;
if (ix >> 63 && checkint (iy) == 1)
{
x2 = -x2;
sign_bias = 1;
}
#if WANT_SIMD_EXCEPT
if (2 * ix == 0 && iy >> 63)
return __math_divzero (sign_bias);
#endif
return iy >> 63 ? 1 / x2 : x2;
}
/* Here x and y are non-zero finite. */
if (ix >> 63)
{
/* Finite x < 0. */
int yint = checkint (iy);
if (yint == 0)
#if WANT_SIMD_EXCEPT
return __math_invalid (x);
#else
return __builtin_nan ("");
#endif
if (yint == 1)
sign_bias = SignBias;
ix &= 0x7fffffffffffffff;
topx &= 0x7ff;
}
if ((topy & 0x7ff) - SmallPowY >= ThresPowY)
{
/* Note: sign_bias == 0 here because y is not odd. */
if (ix == asuint64 (1.0))
return 1.0;
/* |y| < 2^-65, x^y ~= 1 + y*log(x). */
if ((topy & 0x7ff) < SmallPowY)
return 1.0;
#if WANT_SIMD_EXCEPT
return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0)
: __math_uflow (0);
#else
return (ix > asuint64 (1.0)) == (topy < 0x800) ? INFINITY : 0;
#endif
}
if (topx == 0)
{
/* Normalize subnormal x so exponent becomes negative. */
ix = asuint64 (x * 0x1p52);
ix &= 0x7fffffffffffffff;
ix -= 52ULL << 52;
}
}
double lo;
double hi = log_inline (ix, &lo);
double ehi = y * hi;
double elo = y * lo + fma (y, hi, -ehi);
return exp_inline (ehi, elo, sign_bias);
}