mirror of git://sourceware.org/git/glibc.git
Aarch64: Improve codegen in SVE exp and users, and update expf_inline
Use unpredicted muls, and improve memory access. 7%, 3% and 1% improvement in throughput microbenchmark on Neoverse V1, for exp, exp2 and cosh respectively. Reviewed-by: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
This commit is contained in:
Luna Lamb
Wilco Dijkstra
parent
8f0e7fe61e
commit
c0ff447edf
|
|
@ -23,7 +23,7 @@ static const struct data
|
|||
{
|
||||
float64_t poly[3];
|
||||
float64_t inv_ln2, ln2_hi, ln2_lo, shift, thres;
|
||||
uint64_t index_mask, special_bound;
|
||||
uint64_t special_bound;
|
||||
} data = {
|
||||
.poly = { 0x1.fffffffffffd4p-2, 0x1.5555571d6b68cp-3,
|
||||
0x1.5555576a59599p-5, },
|
||||
|
|
@ -35,14 +35,16 @@ static const struct data
|
|||
.shift = 0x1.8p+52,
|
||||
.thres = 704.0,
|
||||
|
||||
.index_mask = 0xff,
|
||||
/* 0x1.6p9, above which exp overflows. */
|
||||
.special_bound = 0x4086000000000000,
|
||||
};
|
||||
|
||||
static svfloat64_t NOINLINE
|
||||
special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
|
||||
special_case (svfloat64_t x, svbool_t pg, svfloat64_t t, svbool_t special)
|
||||
{
|
||||
svfloat64_t half_t = svmul_x (svptrue_b64 (), t, 0.5);
|
||||
svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
|
||||
svfloat64_t y = svadd_x (pg, half_t, half_over_t);
|
||||
return sv_call_f64 (cosh, x, y, special);
|
||||
}
|
||||
|
||||
|
|
@ -60,12 +62,12 @@ exp_inline (svfloat64_t x, const svbool_t pg, const struct data *d)
|
|||
|
||||
svuint64_t u = svreinterpret_u64 (z);
|
||||
svuint64_t e = svlsl_x (pg, u, 52 - V_EXP_TAIL_TABLE_BITS);
|
||||
svuint64_t i = svand_x (pg, u, d->index_mask);
|
||||
svuint64_t i = svand_x (svptrue_b64 (), u, 0xff);
|
||||
|
||||
svfloat64_t y = svmla_x (pg, sv_f64 (d->poly[1]), r, d->poly[2]);
|
||||
y = svmla_x (pg, sv_f64 (d->poly[0]), r, y);
|
||||
y = svmla_x (pg, sv_f64 (1.0), r, y);
|
||||
y = svmul_x (pg, r, y);
|
||||
y = svmul_x (svptrue_b64 (), r, y);
|
||||
|
||||
/* s = 2^(n/N). */
|
||||
u = svld1_gather_index (pg, __v_exp_tail_data, i);
|
||||
|
|
@ -94,12 +96,12 @@ svfloat64_t SV_NAME_D1 (cosh) (svfloat64_t x, const svbool_t pg)
|
|||
/* Up to the point that exp overflows, we can use it to calculate cosh by
|
||||
exp(|x|) / 2 + 1 / (2 * exp(|x|)). */
|
||||
svfloat64_t t = exp_inline (ax, pg, d);
|
||||
svfloat64_t half_t = svmul_x (pg, t, 0.5);
|
||||
svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
|
||||
|
||||
/* Fall back to scalar for any special cases. */
|
||||
if (__glibc_unlikely (svptest_any (pg, special)))
|
||||
return special_case (x, svadd_x (pg, half_t, half_over_t), special);
|
||||
return special_case (x, pg, t, special);
|
||||
|
||||
svfloat64_t half_t = svmul_x (svptrue_b64 (), t, 0.5);
|
||||
svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
|
||||
return svadd_x (pg, half_t, half_over_t);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -18,21 +18,23 @@
|
|||
<https://www.gnu.org/licenses/>. */
|
||||
|
||||
#include "sv_math.h"
|
||||
#include "poly_sve_f64.h"
|
||||
|
||||
#define SpecialBound 307.0 /* floor (log10 (2^1023)). */
|
||||
|
||||
static const struct data
|
||||
{
|
||||
double poly[5];
|
||||
double c1, c3, c2, c4, c0;
|
||||
double shift, log10_2, log2_10_hi, log2_10_lo, scale_thres, special_bound;
|
||||
} data = {
|
||||
/* Coefficients generated using Remez algorithm.
|
||||
rel error: 0x1.9fcb9b3p-60
|
||||
abs error: 0x1.a20d9598p-60 in [ -log10(2)/128, log10(2)/128 ]
|
||||
max ulp err 0.52 +0.5. */
|
||||
.poly = { 0x1.26bb1bbb55516p1, 0x1.53524c73cd32ap1, 0x1.0470591daeafbp1,
|
||||
0x1.2bd77b1361ef6p0, 0x1.142b5d54e9621p-1 },
|
||||
.c0 = 0x1.26bb1bbb55516p1,
|
||||
.c1 = 0x1.53524c73cd32ap1,
|
||||
.c2 = 0x1.0470591daeafbp1,
|
||||
.c3 = 0x1.2bd77b1361ef6p0,
|
||||
.c4 = 0x1.142b5d54e9621p-1,
|
||||
/* 1.5*2^46+1023. This value is further explained below. */
|
||||
.shift = 0x1.800000000ffc0p+46,
|
||||
.log10_2 = 0x1.a934f0979a371p1, /* 1/log2(10). */
|
||||
|
|
@ -70,9 +72,9 @@ special_case (svbool_t pg, svfloat64_t s, svfloat64_t y, svfloat64_t n,
|
|||
/* |n| > 1280 => 2^(n) overflows. */
|
||||
svbool_t p_cmp = svacgt (pg, n, d->scale_thres);
|
||||
|
||||
svfloat64_t r1 = svmul_x (pg, s1, s1);
|
||||
svfloat64_t r1 = svmul_x (svptrue_b64 (), s1, s1);
|
||||
svfloat64_t r2 = svmla_x (pg, s2, s2, y);
|
||||
svfloat64_t r0 = svmul_x (pg, r2, s1);
|
||||
svfloat64_t r0 = svmul_x (svptrue_b64 (), r2, s1);
|
||||
|
||||
return svsel (p_cmp, r1, r0);
|
||||
}
|
||||
|
|
@ -103,11 +105,14 @@ svfloat64_t SV_NAME_D1 (exp10) (svfloat64_t x, svbool_t pg)
|
|||
comes at significant performance cost. */
|
||||
svuint64_t u = svreinterpret_u64 (z);
|
||||
svfloat64_t scale = svexpa (u);
|
||||
|
||||
svfloat64_t c24 = svld1rq (svptrue_b64 (), &d->c2);
|
||||
/* Approximate exp10(r) using polynomial. */
|
||||
svfloat64_t r2 = svmul_x (pg, r, r);
|
||||
svfloat64_t y = svmla_x (pg, svmul_x (pg, r, d->poly[0]), r2,
|
||||
sv_pairwise_poly_3_f64_x (pg, r, r2, d->poly + 1));
|
||||
svfloat64_t r2 = svmul_x (svptrue_b64 (), r, r);
|
||||
svfloat64_t p12 = svmla_lane (sv_f64 (d->c1), r, c24, 0);
|
||||
svfloat64_t p34 = svmla_lane (sv_f64 (d->c3), r, c24, 1);
|
||||
svfloat64_t p14 = svmla_x (pg, p12, p34, r2);
|
||||
|
||||
svfloat64_t y = svmla_x (pg, svmul_x (svptrue_b64 (), r, d->c0), r2, p14);
|
||||
|
||||
/* Assemble result as exp10(x) = 2^n * exp10(r). If |x| > SpecialBound
|
||||
multiplication may overflow, so use special case routine. */
|
||||
|
|
|
|||
|
|
@ -18,7 +18,6 @@
|
|||
<https://www.gnu.org/licenses/>. */
|
||||
|
||||
#include "sv_math.h"
|
||||
#include "poly_sve_f64.h"
|
||||
|
||||
#define N (1 << V_EXP_TABLE_BITS)
|
||||
|
||||
|
|
@ -27,15 +26,15 @@
|
|||
|
||||
static const struct data
|
||||
{
|
||||
double poly[4];
|
||||
double c0, c2;
|
||||
double c1, c3;
|
||||
double shift, big_bound, uoflow_bound;
|
||||
} data = {
|
||||
/* Coefficients are computed using Remez algorithm with
|
||||
minimisation of the absolute error. */
|
||||
.poly = { 0x1.62e42fefa3686p-1, 0x1.ebfbdff82c241p-3, 0x1.c6b09b16de99ap-5,
|
||||
0x1.3b2abf5571ad8p-7 },
|
||||
.shift = 0x1.8p52 / N,
|
||||
.uoflow_bound = UOFlowBound,
|
||||
.c0 = 0x1.62e42fefa3686p-1, .c1 = 0x1.ebfbdff82c241p-3,
|
||||
.c2 = 0x1.c6b09b16de99ap-5, .c3 = 0x1.3b2abf5571ad8p-7,
|
||||
.shift = 0x1.8p52 / N, .uoflow_bound = UOFlowBound,
|
||||
.big_bound = BigBound,
|
||||
};
|
||||
|
||||
|
|
@ -67,9 +66,9 @@ special_case (svbool_t pg, svfloat64_t s, svfloat64_t y, svfloat64_t n,
|
|||
/* |n| > 1280 => 2^(n) overflows. */
|
||||
svbool_t p_cmp = svacgt (pg, n, d->uoflow_bound);
|
||||
|
||||
svfloat64_t r1 = svmul_x (pg, s1, s1);
|
||||
svfloat64_t r1 = svmul_x (svptrue_b64 (), s1, s1);
|
||||
svfloat64_t r2 = svmla_x (pg, s2, s2, y);
|
||||
svfloat64_t r0 = svmul_x (pg, r2, s1);
|
||||
svfloat64_t r0 = svmul_x (svptrue_b64 (), r2, s1);
|
||||
|
||||
return svsel (p_cmp, r1, r0);
|
||||
}
|
||||
|
|
@ -99,11 +98,14 @@ svfloat64_t SV_NAME_D1 (exp2) (svfloat64_t x, svbool_t pg)
|
|||
svuint64_t top = svlsl_x (pg, ki, 52 - V_EXP_TABLE_BITS);
|
||||
svfloat64_t scale = svreinterpret_f64 (svadd_x (pg, sbits, top));
|
||||
|
||||
svfloat64_t c13 = svld1rq (svptrue_b64 (), &d->c1);
|
||||
/* Approximate exp2(r) using polynomial. */
|
||||
svfloat64_t r2 = svmul_x (pg, r, r);
|
||||
svfloat64_t p = sv_pairwise_poly_3_f64_x (pg, r, r2, d->poly);
|
||||
svfloat64_t y = svmul_x (pg, r, p);
|
||||
|
||||
/* y = exp2(r) - 1 ~= C0 r + C1 r^2 + C2 r^3 + C3 r^4. */
|
||||
svfloat64_t r2 = svmul_x (svptrue_b64 (), r, r);
|
||||
svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), r, c13, 0);
|
||||
svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), r, c13, 1);
|
||||
svfloat64_t p = svmla_x (pg, p01, p23, r2);
|
||||
svfloat64_t y = svmul_x (svptrue_b64 (), r, p);
|
||||
/* Assemble exp2(x) = exp2(r) * scale. */
|
||||
if (__glibc_unlikely (svptest_any (pg, special)))
|
||||
return special_case (pg, scale, y, kd, d);
|
||||
|
|
|
|||
|
|
@ -21,12 +21,15 @@
|
|||
|
||||
static const struct data
|
||||
{
|
||||
double poly[4];
|
||||
double c0, c2;
|
||||
double c1, c3;
|
||||
double ln2_hi, ln2_lo, inv_ln2, shift, thres;
|
||||
|
||||
} data = {
|
||||
.poly = { /* ulp error: 0.53. */
|
||||
0x1.fffffffffdbcdp-2, 0x1.555555555444cp-3, 0x1.555573c6a9f7dp-5,
|
||||
0x1.1111266d28935p-7 },
|
||||
.c0 = 0x1.fffffffffdbcdp-2,
|
||||
.c1 = 0x1.555555555444cp-3,
|
||||
.c2 = 0x1.555573c6a9f7dp-5,
|
||||
.c3 = 0x1.1111266d28935p-7,
|
||||
.ln2_hi = 0x1.62e42fefa3800p-1,
|
||||
.ln2_lo = 0x1.ef35793c76730p-45,
|
||||
/* 1/ln2. */
|
||||
|
|
@ -36,7 +39,6 @@ static const struct data
|
|||
.thres = 704.0,
|
||||
};
|
||||
|
||||
#define C(i) sv_f64 (d->poly[i])
|
||||
#define SpecialOffset 0x6000000000000000 /* 0x1p513. */
|
||||
/* SpecialBias1 + SpecialBias1 = asuint(1.0). */
|
||||
#define SpecialBias1 0x7000000000000000 /* 0x1p769. */
|
||||
|
|
@ -56,20 +58,20 @@ special_case (svbool_t pg, svfloat64_t s, svfloat64_t y, svfloat64_t n)
|
|||
svuint64_t b
|
||||
= svdup_u64_z (p_sign, SpecialOffset); /* Inactive lanes set to 0. */
|
||||
|
||||
/* Set s1 to generate overflow depending on sign of exponent n. */
|
||||
svfloat64_t s1 = svreinterpret_f64 (
|
||||
svsubr_x (pg, b, SpecialBias1)); /* 0x70...0 - b. */
|
||||
/* Offset s to avoid overflow in final result if n is below threshold. */
|
||||
/* Set s1 to generate overflow depending on sign of exponent n,
|
||||
ie. s1 = 0x70...0 - b. */
|
||||
svfloat64_t s1 = svreinterpret_f64 (svsubr_x (pg, b, SpecialBias1));
|
||||
/* Offset s to avoid overflow in final result if n is below threshold.
|
||||
ie. s2 = as_u64 (s) - 0x3010...0 + b. */
|
||||
svfloat64_t s2 = svreinterpret_f64 (
|
||||
svadd_x (pg, svsub_x (pg, svreinterpret_u64 (s), SpecialBias2),
|
||||
b)); /* as_u64 (s) - 0x3010...0 + b. */
|
||||
svadd_x (pg, svsub_x (pg, svreinterpret_u64 (s), SpecialBias2), b));
|
||||
|
||||
/* |n| > 1280 => 2^(n) overflows. */
|
||||
svbool_t p_cmp = svacgt (pg, n, 1280.0);
|
||||
|
||||
svfloat64_t r1 = svmul_x (pg, s1, s1);
|
||||
svfloat64_t r1 = svmul_x (svptrue_b64 (), s1, s1);
|
||||
svfloat64_t r2 = svmla_x (pg, s2, s2, y);
|
||||
svfloat64_t r0 = svmul_x (pg, r2, s1);
|
||||
svfloat64_t r0 = svmul_x (svptrue_b64 (), r2, s1);
|
||||
|
||||
return svsel (p_cmp, r1, r0);
|
||||
}
|
||||
|
|
@ -103,16 +105,16 @@ svfloat64_t SV_NAME_D1 (exp) (svfloat64_t x, const svbool_t pg)
|
|||
svfloat64_t z = svmla_x (pg, sv_f64 (d->shift), x, d->inv_ln2);
|
||||
svuint64_t u = svreinterpret_u64 (z);
|
||||
svfloat64_t n = svsub_x (pg, z, d->shift);
|
||||
|
||||
svfloat64_t c13 = svld1rq (svptrue_b64 (), &d->c1);
|
||||
/* r = x - n * ln2, r is in [-ln2/(2N), ln2/(2N)]. */
|
||||
svfloat64_t ln2 = svld1rq (svptrue_b64 (), &d->ln2_hi);
|
||||
svfloat64_t r = svmls_lane (x, n, ln2, 0);
|
||||
r = svmls_lane (r, n, ln2, 1);
|
||||
|
||||
/* y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5. */
|
||||
svfloat64_t r2 = svmul_x (pg, r, r);
|
||||
svfloat64_t p01 = svmla_x (pg, C (0), C (1), r);
|
||||
svfloat64_t p23 = svmla_x (pg, C (2), C (3), r);
|
||||
svfloat64_t r2 = svmul_x (svptrue_b64 (), r, r);
|
||||
svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), r, c13, 0);
|
||||
svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), r, c13, 1);
|
||||
svfloat64_t p04 = svmla_x (pg, p01, p23, r2);
|
||||
svfloat64_t y = svmla_x (pg, r, p04, r2);
|
||||
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ expf_inline (svfloat32_t x, const svbool_t pg, const struct sv_expf_data *d)
|
|||
/* scale = 2^(n/N). */
|
||||
svfloat32_t scale = svexpa (svreinterpret_u32 (z));
|
||||
|
||||
/* y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5 + C4 r^6. */
|
||||
/* poly(r) = exp(r) - 1 ~= C0 r + C1 r^2 + C2 r^3 + C3 r^4 + C4 r^5. */
|
||||
svfloat32_t p12 = svmla_lane (sv_f32 (d->c1), r, lane_consts, 2);
|
||||
svfloat32_t p34 = svmla_lane (sv_f32 (d->c3), r, lane_consts, 3);
|
||||
svfloat32_t r2 = svmul_x (svptrue_b32 (), r, r);
|
||||
|
|
@ -71,5 +71,4 @@ expf_inline (svfloat32_t x, const svbool_t pg, const struct sv_expf_data *d)
|
|||
|
||||
return svmla_x (pg, scale, scale, poly);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
|
|||
Loading…
Reference in New Issue