AArch64: Improve codegen for SVE pow

Move constants to struct.  Improve memory access with indexed/unpredicated
instructions.  Eliminate register spills.  Speedup on Neoverse V1: 24%.

Reviewed-by: Wilco Dijkstra  <Wilco.Dijkstra@arm.com>
(cherry picked from commit 0b195651db3ae793187c7dd6d78b5a7a8da9d5e6)
This commit is contained in:
Yat Long Poon 2025-02-13 18:02:01 +00:00 committed by Wilco Dijkstra
parent f3d9c116cb
commit 2025e27a81

View File

@ -44,19 +44,18 @@
/* Data is defined in v_pow_log_data.c. */
#define N_LOG (1 << V_POW_LOG_TABLE_BITS)
#define A __v_pow_log_data.poly
#define Off 0x3fe6955500000000
/* 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 C __v_pow_exp_data.poly
#define SmallExp 0x3c9 /* top12(0x1p-54). */
#define BigExp 0x408 /* top12(512.). */
#define ThresExp 0x03f /* BigExp - SmallExp. */
#define HugeExp 0x409 /* top12(1024.). */
/* Constants associated with pow. */
#define SmallBoundX 0x1p-126
#define SmallPowX 0x001 /* top12(0x1p-126). */
#define BigPowX 0x7ff /* top12(INFINITY). */
#define ThresPowX 0x7fe /* BigPowX - SmallPowX. */
@ -64,6 +63,31 @@
#define BigPowY 0x43e /* top12(0x1.749p62). */
#define ThresPowY 0x080 /* BigPowY - SmallPowY. */
static const struct data
{
double log_c0, log_c2, log_c4, log_c6, ln2_hi, ln2_lo;
double log_c1, log_c3, log_c5, off;
double n_over_ln2, exp_c2, ln2_over_n_hi, ln2_over_n_lo;
double exp_c0, exp_c1;
} data = {
.log_c0 = -0x1p-1,
.log_c1 = -0x1.555555555556p-1,
.log_c2 = 0x1.0000000000006p-1,
.log_c3 = 0x1.999999959554ep-1,
.log_c4 = -0x1.555555529a47ap-1,
.log_c5 = -0x1.2495b9b4845e9p0,
.log_c6 = 0x1.0002b8b263fc3p0,
.off = Off,
.exp_c0 = 0x1.fffffffffffd4p-2,
.exp_c1 = 0x1.5555571d6ef9p-3,
.exp_c2 = 0x1.5555576a5adcep-5,
.ln2_hi = 0x1.62e42fefa3800p-1,
.ln2_lo = 0x1.ef35793c76730p-45,
.n_over_ln2 = 0x1.71547652b82fep0 * N_EXP,
.ln2_over_n_hi = 0x1.62e42fefc0000p-9,
.ln2_over_n_lo = -0x1.c610ca86c3899p-45,
};
/* Check if x is an integer. */
static inline svbool_t
sv_isint (svbool_t pg, svfloat64_t x)
@ -82,7 +106,7 @@ sv_isnotint (svbool_t pg, svfloat64_t x)
static inline svbool_t
sv_isodd (svbool_t pg, svfloat64_t x)
{
svfloat64_t y = svmul_x (pg, x, 0.5);
svfloat64_t y = svmul_x (svptrue_b64 (), x, 0.5);
return sv_isnotint (pg, y);
}
@ -121,7 +145,7 @@ zeroinfnan (uint64_t i)
static inline svbool_t
sv_zeroinfnan (svbool_t pg, svuint64_t i)
{
return svcmpge (pg, svsub_x (pg, svmul_x (pg, i, 2), 1),
return svcmpge (pg, svsub_x (pg, svadd_x (pg, i, i), 1),
2 * asuint64 (INFINITY) - 1);
}
@ -174,16 +198,17 @@ sv_call_specialcase (svfloat64_t x1, svuint64_t u1, svuint64_t u2,
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 svfloat64_t
sv_log_inline (svbool_t pg, svuint64_t ix, svfloat64_t *tail)
sv_log_inline (svbool_t pg, svuint64_t ix, svfloat64_t *tail,
const struct data *d)
{
/* 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. */
svuint64_t tmp = svsub_x (pg, ix, Off);
svuint64_t tmp = svsub_x (pg, ix, d->off);
svuint64_t i = svand_x (pg, svlsr_x (pg, tmp, 52 - V_POW_LOG_TABLE_BITS),
sv_u64 (N_LOG - 1));
svint64_t k = svasr_x (pg, svreinterpret_s64 (tmp), 52);
svuint64_t iz = svsub_x (pg, ix, svand_x (pg, tmp, sv_u64 (0xfffULL << 52)));
svuint64_t iz = svsub_x (pg, ix, svlsl_x (pg, svreinterpret_u64 (k), 52));
svfloat64_t z = svreinterpret_f64 (iz);
svfloat64_t kd = svcvt_f64_x (pg, k);
@ -199,40 +224,85 @@ sv_log_inline (svbool_t pg, svuint64_t ix, svfloat64_t *tail)
|z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
svfloat64_t r = svmad_x (pg, z, invc, -1.0);
/* k*Ln2 + log(c) + r. */
svfloat64_t t1 = svmla_x (pg, logc, kd, __v_pow_log_data.ln2_hi);
svfloat64_t ln2_hilo = svld1rq_f64 (svptrue_b64 (), &d->ln2_hi);
svfloat64_t t1 = svmla_lane_f64 (logc, kd, ln2_hilo, 0);
svfloat64_t t2 = svadd_x (pg, t1, r);
svfloat64_t lo1 = svmla_x (pg, logctail, kd, __v_pow_log_data.ln2_lo);
svfloat64_t lo1 = svmla_lane_f64 (logctail, kd, ln2_hilo, 1);
svfloat64_t lo2 = svadd_x (pg, svsub_x (pg, t1, t2), r);
/* Evaluation is optimized assuming superscalar pipelined execution. */
svfloat64_t ar = svmul_x (pg, r, -0.5); /* A[0] = -0.5. */
svfloat64_t ar2 = svmul_x (pg, r, ar);
svfloat64_t ar3 = svmul_x (pg, r, ar2);
svfloat64_t log_c02 = svld1rq_f64 (svptrue_b64 (), &d->log_c0);
svfloat64_t ar = svmul_lane_f64 (r, log_c02, 0);
svfloat64_t ar2 = svmul_x (svptrue_b64 (), r, ar);
svfloat64_t ar3 = svmul_x (svptrue_b64 (), r, ar2);
/* k*Ln2 + log(c) + r + A[0]*r*r. */
svfloat64_t hi = svadd_x (pg, t2, ar2);
svfloat64_t lo3 = svmla_x (pg, svneg_x (pg, ar2), ar, r);
svfloat64_t lo3 = svmls_x (pg, ar2, ar, r);
svfloat64_t lo4 = svadd_x (pg, svsub_x (pg, t2, hi), ar2);
/* p = log1p(r) - r - A[0]*r*r. */
/* p = (ar3 * (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r *
A[6])))). */
svfloat64_t a56 = svmla_x (pg, sv_f64 (A[5]), r, A[6]);
svfloat64_t a34 = svmla_x (pg, sv_f64 (A[3]), r, A[4]);
svfloat64_t a12 = svmla_x (pg, sv_f64 (A[1]), r, A[2]);
svfloat64_t log_c46 = svld1rq_f64 (svptrue_b64 (), &d->log_c4);
svfloat64_t a56 = svmla_lane_f64 (sv_f64 (d->log_c5), r, log_c46, 1);
svfloat64_t a34 = svmla_lane_f64 (sv_f64 (d->log_c3), r, log_c46, 0);
svfloat64_t a12 = svmla_lane_f64 (sv_f64 (d->log_c1), r, log_c02, 1);
svfloat64_t p = svmla_x (pg, a34, ar2, a56);
p = svmla_x (pg, a12, ar2, p);
p = svmul_x (pg, ar3, p);
p = svmul_x (svptrue_b64 (), ar3, p);
svfloat64_t lo = svadd_x (
pg, svadd_x (pg, svadd_x (pg, svadd_x (pg, lo1, lo2), lo3), lo4), p);
pg, svadd_x (pg, svsub_x (pg, svadd_x (pg, lo1, lo2), lo3), lo4), p);
svfloat64_t y = svadd_x (pg, hi, lo);
*tail = svadd_x (pg, svsub_x (pg, hi, y), lo);
return y;
}
static inline svfloat64_t
sv_exp_core (svbool_t pg, svfloat64_t x, svfloat64_t xtail,
svuint64_t sign_bias, svfloat64_t *tmp, svuint64_t *sbits,
svuint64_t *ki, const struct data *d)
{
/* 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]. */
svfloat64_t n_over_ln2_and_c2 = svld1rq_f64 (svptrue_b64 (), &d->n_over_ln2);
svfloat64_t z = svmul_lane_f64 (x, n_over_ln2_and_c2, 0);
/* z - kd is in [-1, 1] in non-nearest rounding modes. */
svfloat64_t kd = svrinta_x (pg, z);
*ki = svreinterpret_u64 (svcvt_s64_x (pg, kd));
svfloat64_t ln2_over_n_hilo
= svld1rq_f64 (svptrue_b64 (), &d->ln2_over_n_hi);
svfloat64_t r = x;
r = svmls_lane_f64 (r, kd, ln2_over_n_hilo, 0);
r = svmls_lane_f64 (r, kd, ln2_over_n_hilo, 1);
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
r = svadd_x (pg, r, xtail);
/* 2^(k/N) ~= scale. */
svuint64_t idx = svand_x (pg, *ki, N_EXP - 1);
svuint64_t top
= svlsl_x (pg, svadd_x (pg, *ki, sign_bias), 52 - V_POW_EXP_TABLE_BITS);
/* This is only a valid scale when -1023*N < k < 1024*N. */
*sbits = svld1_gather_index (pg, __v_pow_exp_data.sbits, idx);
*sbits = svadd_x (pg, *sbits, top);
/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */
svfloat64_t r2 = svmul_x (svptrue_b64 (), r, r);
*tmp = svmla_lane_f64 (sv_f64 (d->exp_c1), r, n_over_ln2_and_c2, 1);
*tmp = svmla_x (pg, sv_f64 (d->exp_c0), r, *tmp);
*tmp = svmla_x (pg, r, r2, *tmp);
svfloat64_t scale = svreinterpret_f64 (*sbits);
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
is no spurious underflow here even without fma. */
z = svmla_x (pg, scale, scale, *tmp);
return z;
}
/* 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 svfloat64_t
sv_exp_inline (svbool_t pg, svfloat64_t x, svfloat64_t xtail,
svuint64_t sign_bias)
svuint64_t sign_bias, const struct data *d)
{
/* 3 types of special cases: tiny (uflow and spurious uflow), huge (oflow)
and other cases of large values of x (scale * (1 + TMP) oflow). */
@ -240,73 +310,46 @@ sv_exp_inline (svbool_t pg, svfloat64_t x, svfloat64_t xtail,
/* |x| is large (|x| >= 512) or tiny (|x| <= 0x1p-54). */
svbool_t uoflow = svcmpge (pg, svsub_x (pg, abstop, SmallExp), ThresExp);
/* Conditions special, uflow and oflow are all expressed as uoflow &&
something, hence do not bother computing anything if no lane in uoflow is
true. */
svbool_t special = svpfalse_b ();
svbool_t uflow = svpfalse_b ();
svbool_t oflow = svpfalse_b ();
svfloat64_t tmp;
svuint64_t sbits, ki;
if (__glibc_unlikely (svptest_any (pg, uoflow)))
{
svfloat64_t z
= sv_exp_core (pg, x, xtail, sign_bias, &tmp, &sbits, &ki, d);
/* |x| is tiny (|x| <= 0x1p-54). */
uflow = svcmpge (pg, svsub_x (pg, abstop, SmallExp), 0x80000000);
svbool_t uflow
= svcmpge (pg, svsub_x (pg, abstop, SmallExp), 0x80000000);
uflow = svand_z (pg, uoflow, uflow);
/* |x| is huge (|x| >= 1024). */
oflow = svcmpge (pg, abstop, HugeExp);
svbool_t oflow = svcmpge (pg, abstop, HugeExp);
oflow = svand_z (pg, uoflow, svbic_z (pg, oflow, uflow));
/* For large |x| values (512 < |x| < 1024) scale * (1 + TMP) can overflow
or underflow. */
special = svbic_z (pg, uoflow, svorr_z (pg, uflow, oflow));
or underflow. */
svbool_t special = svbic_z (pg, uoflow, svorr_z (pg, uflow, oflow));
/* Update result with special and large cases. */
z = sv_call_specialcase (tmp, sbits, ki, z, special);
/* Handle underflow and overflow. */
svbool_t x_is_neg = svcmplt (pg, x, 0);
svuint64_t sign_mask
= svlsl_x (pg, sign_bias, 52 - V_POW_EXP_TABLE_BITS);
svfloat64_t res_uoflow
= svsel (x_is_neg, sv_f64 (0.0), sv_f64 (INFINITY));
res_uoflow = svreinterpret_f64 (
svorr_x (pg, svreinterpret_u64 (res_uoflow), sign_mask));
/* Avoid spurious underflow for tiny x. */
svfloat64_t res_spurious_uflow
= svreinterpret_f64 (svorr_x (pg, sign_mask, 0x3ff0000000000000));
z = svsel (oflow, res_uoflow, z);
z = svsel (uflow, res_spurious_uflow, z);
return z;
}
/* 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]. */
svfloat64_t z = svmul_x (pg, x, __v_pow_exp_data.n_over_ln2);
/* z - kd is in [-1, 1] in non-nearest rounding modes. */
svfloat64_t shift = sv_f64 (__v_pow_exp_data.shift);
svfloat64_t kd = svadd_x (pg, z, shift);
svuint64_t ki = svreinterpret_u64 (kd);
kd = svsub_x (pg, kd, shift);
svfloat64_t r = x;
r = svmls_x (pg, r, kd, __v_pow_exp_data.ln2_over_n_hi);
r = svmls_x (pg, r, kd, __v_pow_exp_data.ln2_over_n_lo);
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
r = svadd_x (pg, r, xtail);
/* 2^(k/N) ~= scale. */
svuint64_t idx = svand_x (pg, ki, N_EXP - 1);
svuint64_t top
= svlsl_x (pg, svadd_x (pg, ki, sign_bias), 52 - V_POW_EXP_TABLE_BITS);
/* This is only a valid scale when -1023*N < k < 1024*N. */
svuint64_t sbits = svld1_gather_index (pg, __v_pow_exp_data.sbits, idx);
sbits = svadd_x (pg, sbits, top);
/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */
svfloat64_t r2 = svmul_x (pg, r, r);
svfloat64_t tmp = svmla_x (pg, sv_f64 (C[1]), r, C[2]);
tmp = svmla_x (pg, sv_f64 (C[0]), r, tmp);
tmp = svmla_x (pg, r, r2, tmp);
svfloat64_t scale = svreinterpret_f64 (sbits);
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
is no spurious underflow here even without fma. */
z = svmla_x (pg, scale, scale, tmp);
/* Update result with special and large cases. */
if (__glibc_unlikely (svptest_any (pg, special)))
z = sv_call_specialcase (tmp, sbits, ki, z, special);
/* Handle underflow and overflow. */
svuint64_t sign_bit = svlsr_x (pg, svreinterpret_u64 (x), 63);
svbool_t x_is_neg = svcmpne (pg, sign_bit, 0);
svuint64_t sign_mask = svlsl_x (pg, sign_bias, 52 - V_POW_EXP_TABLE_BITS);
svfloat64_t res_uoflow = svsel (x_is_neg, sv_f64 (0.0), sv_f64 (INFINITY));
res_uoflow = svreinterpret_f64 (
svorr_x (pg, svreinterpret_u64 (res_uoflow), sign_mask));
z = svsel (oflow, res_uoflow, z);
/* Avoid spurious underflow for tiny x. */
svfloat64_t res_spurious_uflow
= svreinterpret_f64 (svorr_x (pg, sign_mask, 0x3ff0000000000000));
z = svsel (uflow, res_spurious_uflow, z);
return z;
return sv_exp_core (pg, x, xtail, sign_bias, &tmp, &sbits, &ki, d);
}
static inline double
@ -341,47 +384,39 @@ pow_sc (double x, double y)
svfloat64_t SV_NAME_D2 (pow) (svfloat64_t x, svfloat64_t y, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
/* This preamble handles special case conditions used in the final scalar
fallbacks. It also updates ix and sign_bias, that are used in the core
computation too, i.e., exp( y * log (x) ). */
svuint64_t vix0 = svreinterpret_u64 (x);
svuint64_t viy0 = svreinterpret_u64 (y);
svuint64_t vtopx0 = svlsr_x (svptrue_b64 (), vix0, 52);
/* Negative x cases. */
svuint64_t sign_bit = svlsr_m (pg, vix0, 63);
svbool_t xisneg = svcmpeq (pg, sign_bit, 1);
svbool_t xisneg = svcmplt (pg, x, 0);
/* Set sign_bias and ix depending on sign of x and nature of y. */
svbool_t yisnotint_xisneg = svpfalse_b ();
svbool_t yint_or_xpos = pg;
svuint64_t sign_bias = sv_u64 (0);
svuint64_t vix = vix0;
svuint64_t vtopx1 = vtopx0;
if (__glibc_unlikely (svptest_any (pg, xisneg)))
{
/* Determine nature of y. */
yisnotint_xisneg = sv_isnotint (xisneg, y);
svbool_t yisint_xisneg = sv_isint (xisneg, y);
yint_or_xpos = sv_isint (xisneg, y);
svbool_t yisodd_xisneg = sv_isodd (xisneg, y);
/* ix set to abs(ix) if y is integer. */
vix = svand_m (yisint_xisneg, vix0, 0x7fffffffffffffff);
vtopx1 = svand_m (yisint_xisneg, vtopx0, 0x7ff);
vix = svand_m (yint_or_xpos, vix0, 0x7fffffffffffffff);
/* Set to SignBias if x is negative and y is odd. */
sign_bias = svsel (yisodd_xisneg, sv_u64 (SignBias), sv_u64 (0));
}
/* Special cases of x or y: zero, inf and nan. */
svbool_t xspecial = sv_zeroinfnan (pg, vix0);
svbool_t yspecial = sv_zeroinfnan (pg, viy0);
svbool_t special = svorr_z (pg, xspecial, yspecial);
/* Small cases of x: |x| < 0x1p-126. */
svuint64_t vabstopx0 = svand_x (pg, vtopx0, 0x7ff);
svbool_t xsmall = svcmplt (pg, vabstopx0, SmallPowX);
if (__glibc_unlikely (svptest_any (pg, xsmall)))
svbool_t xsmall = svaclt (yint_or_xpos, x, SmallBoundX);
if (__glibc_unlikely (svptest_any (yint_or_xpos, xsmall)))
{
/* Normalize subnormal x so exponent becomes negative. */
svbool_t topx_is_null = svcmpeq (xsmall, vtopx1, 0);
svuint64_t vtopx = svlsr_x (svptrue_b64 (), vix, 52);
svbool_t topx_is_null = svcmpeq (xsmall, vtopx, 0);
svuint64_t vix_norm = svreinterpret_u64 (svmul_m (xsmall, x, 0x1p52));
vix_norm = svand_m (xsmall, vix_norm, 0x7fffffffffffffff);
@ -391,20 +426,24 @@ svfloat64_t SV_NAME_D2 (pow) (svfloat64_t x, svfloat64_t y, const svbool_t pg)
/* y_hi = log(ix, &y_lo). */
svfloat64_t vlo;
svfloat64_t vhi = sv_log_inline (pg, vix, &vlo);
svfloat64_t vhi = sv_log_inline (yint_or_xpos, vix, &vlo, d);
/* z = exp(y_hi, y_lo, sign_bias). */
svfloat64_t vehi = svmul_x (pg, y, vhi);
svfloat64_t velo = svmul_x (pg, y, vlo);
svfloat64_t vemi = svmls_x (pg, vehi, y, vhi);
velo = svsub_x (pg, velo, vemi);
svfloat64_t vz = sv_exp_inline (pg, vehi, velo, sign_bias);
svfloat64_t vehi = svmul_x (svptrue_b64 (), y, vhi);
svfloat64_t vemi = svmls_x (yint_or_xpos, vehi, y, vhi);
svfloat64_t velo = svnmls_x (yint_or_xpos, vemi, y, vlo);
svfloat64_t vz = sv_exp_inline (yint_or_xpos, vehi, velo, sign_bias, d);
/* Cases of finite y and finite negative x. */
vz = svsel (yisnotint_xisneg, sv_f64 (__builtin_nan ("")), vz);
vz = svsel (yint_or_xpos, vz, sv_f64 (__builtin_nan ("")));
/* Special cases of x or y: zero, inf and nan. */
svbool_t xspecial = sv_zeroinfnan (svptrue_b64 (), vix0);
svbool_t yspecial = sv_zeroinfnan (svptrue_b64 (), viy0);
svbool_t special = svorr_z (svptrue_b64 (), xspecial, yspecial);
/* Cases of zero/inf/nan x or y. */
if (__glibc_unlikely (svptest_any (pg, special)))
if (__glibc_unlikely (svptest_any (svptrue_b64 (), special)))
vz = sv_call2_f64 (pow_sc, x, y, vz, special);
return vz;