From fe5f790b64cfe52abe58c265b0a00534a69b7e80 Mon Sep 17 00:00:00 2001 From: Dylan Fleming Date: Tue, 31 Oct 2023 10:43:46 +0000 Subject: pl/math: Use Neon indexed fmla/mul with constants stored in neon registers. Replaced various vfmaq instructions with vfmaq_laneq variants to reduce constant vector loads. --- pl/math/v_cosh_2u.c | 9 ++++----- pl/math/v_exp10f_2u4.c | 21 +++++++++++---------- pl/math/v_expf_inline.h | 13 ++++++------- pl/math/v_expm1_2u5.c | 9 ++++----- pl/math/v_expm1f_1u6.c | 15 ++++++++------- pl/math/v_expm1f_inline.h | 18 +++++++++--------- pl/math/v_sinh_3u.c | 9 ++++----- pl/math/v_tan_3u5.c | 20 ++++++++++---------- pl/math/v_tanf_3u5.c | 25 +++++++++++++------------ 9 files changed, 69 insertions(+), 70 deletions(-) diff --git a/pl/math/v_cosh_2u.c b/pl/math/v_cosh_2u.c index ccb67d7..060d7ec 100644 --- a/pl/math/v_cosh_2u.c +++ b/pl/math/v_cosh_2u.c @@ -12,7 +12,7 @@ static const struct data { float64x2_t poly[3]; - float64x2_t inv_ln2, ln2_hi, ln2_lo, shift, thres; + float64x2_t inv_ln2, ln2, shift, thres; uint64x2_t index_mask, special_bound; } data = { .poly = { V2 (0x1.fffffffffffd4p-2), V2 (0x1.5555571d6b68cp-3), @@ -20,8 +20,7 @@ static const struct data .inv_ln2 = V2 (0x1.71547652b82fep8), /* N/ln2. */ /* -ln2/N. */ - .ln2_hi = V2 (-0x1.62e42fefa39efp-9), - .ln2_lo = V2 (-0x1.abc9e3b39803f3p-64), + .ln2 = {-0x1.62e42fefa39efp-9, -0x1.abc9e3b39803f3p-64}, .shift = V2 (0x1.8p+52), .thres = V2 (704.0), @@ -49,8 +48,8 @@ exp_inline (float64x2_t x) float64x2_t n = vsubq_f64 (z, d->shift); /* r = x - n*ln2/N. */ - float64x2_t r = vfmaq_f64 (x, d->ln2_hi, n); - r = vfmaq_f64 (r, d->ln2_lo, n); + float64x2_t r = vfmaq_laneq_f64 (x, n, d->ln2, 0); + r = vfmaq_laneq_f64 (r, n, d->ln2, 1); uint64x2_t e = vshlq_n_u64 (u, 52 - V_EXP_TAIL_TABLE_BITS); uint64x2_t i = vandq_u64 (u, d->index_mask); diff --git a/pl/math/v_exp10f_2u4.c b/pl/math/v_exp10f_2u4.c index 8a80fbd..0d37dcb 100644 --- a/pl/math/v_exp10f_2u4.c +++ b/pl/math/v_exp10f_2u4.c @@ -16,7 +16,8 @@ static const struct data { float32x4_t poly[5]; - float32x4_t shift, log10_2, log2_10_hi, log2_10_lo; + float32x4_t log10_2_and_inv, shift; + #if !WANT_SIMD_EXCEPT float32x4_t scale_thresh; #endif @@ -29,9 +30,9 @@ static const struct data .poly = { V4 (0x1.26bb16p+1f), V4 (0x1.5350d2p+1f), V4 (0x1.04744ap+1f), V4 (0x1.2d8176p+0f), V4 (0x1.12b41ap-1f) }, .shift = V4 (0x1.8p23f), - .log10_2 = V4 (0x1.a934fp+1), - .log2_10_hi = V4 (0x1.344136p-2), - .log2_10_lo = V4 (-0x1.ec10cp-27), + + /* Stores constants 1/log10(2), log10(2)_high, log10(2)_low, 0. */ + .log10_2_and_inv = { 0x1.a934fp+1, 0x1.344136p-2, -0x1.ec10cp-27, 0 }, #if !WANT_SIMD_EXCEPT .scale_thresh = V4 (ScaleBound) #endif @@ -56,9 +57,9 @@ special_case (float32x4_t x, float32x4_t y, uint32x4_t cmp) #else -# define SpecialBound 126.0f /* rint (log2 (2^127 / (1 + sqrt (2)))). */ -# define SpecialOffset v_u32 (0x82000000) -# define SpecialBias v_u32 (0x7f000000) +# define SpecialBound 126.0f /* rint (log2 (2^127 / (1 + sqrt (2)))). */ +# define SpecialOffset v_u32 (0x82000000) +# define SpecialBias v_u32 (0x7f000000) static float32x4_t VPCS_ATTR NOINLINE special_case (float32x4_t poly, float32x4_t n, uint32x4_t e, uint32x4_t cmp1, @@ -103,10 +104,10 @@ float32x4_t VPCS_ATTR V_NAME_F1 (exp10) (float32x4_t x) /* exp10(x) = 2^n * 10^r = 2^n * (1 + poly (r)), with poly(r) in [1/sqrt(2), sqrt(2)] and x = r + n * log10 (2), with r in [-log10(2)/2, log10(2)/2]. */ - float32x4_t z = vfmaq_f32 (d->shift, x, d->log10_2); + float32x4_t z = vfmaq_laneq_f32 (d->shift, x, d->log10_2_and_inv, 0); float32x4_t n = vsubq_f32 (z, d->shift); - float32x4_t r = vfmsq_f32 (x, n, d->log2_10_hi); - r = vfmsq_f32 (r, n, d->log2_10_lo); + float32x4_t r = vfmsq_laneq_f32 (x, n, d->log10_2_and_inv, 1); + r = vfmsq_laneq_f32 (r, n, d->log10_2_and_inv, 2); uint32x4_t e = vshlq_n_u32 (vreinterpretq_u32_f32 (z), 23); float32x4_t scale = vreinterpretq_f32_u32 (vaddq_u32 (e, ExponentBias)); diff --git a/pl/math/v_expf_inline.h b/pl/math/v_expf_inline.h index 682cbc5..1666837 100644 --- a/pl/math/v_expf_inline.h +++ b/pl/math/v_expf_inline.h @@ -14,7 +14,7 @@ struct v_expf_data { float32x4_t poly[5]; - float32x4_t shift, invln2, ln2_hi, ln2_lo; + float32x4_t shift, invln2_and_ln2; }; /* maxerr: 1.45358 +0.5 ulp. */ @@ -22,9 +22,8 @@ struct v_expf_data { \ .poly = { V4 (0x1.0e4020p-7f), V4 (0x1.573e2ep-5f), V4 (0x1.555e66p-3f), \ V4 (0x1.fffdb6p-2f), V4 (0x1.ffffecp-1f) }, \ - \ - .shift = V4 (0x1.8p23f), .invln2 = V4 (0x1.715476p+0f), \ - .ln2_hi = V4 (0x1.62e4p-1f), .ln2_lo = V4 (0x1.7f7d1cp-20f), \ + .shift = V4 (0x1.8p23f), \ + .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 }, \ } #define ExponentBias v_u32 (0x3f800000) /* asuint(1.0f). */ @@ -40,10 +39,10 @@ v_expf_inline (float32x4_t x, const struct v_expf_data *d) /* exp(x) = 2^n (1 + poly(r)), with 1 + poly(r) in [1/sqrt(2),sqrt(2)] x = ln2*n + r, with r in [-ln2/2, ln2/2]. */ float32x4_t n, r, z; - z = vfmaq_f32 (d->shift, x, d->invln2); + z = vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0); n = vsubq_f32 (z, d->shift); - r = vfmsq_f32 (x, n, d->ln2_hi); - r = vfmsq_f32 (r, n, d->ln2_lo); + r = vfmsq_laneq_f32 (x, n, d->invln2_and_ln2, 1); + r = vfmsq_laneq_f32 (r, n, d->invln2_and_ln2, 2); uint32x4_t e = vshlq_n_u32 (vreinterpretq_u32_f32 (z), 23); float32x4_t scale = vreinterpretq_f32_u32 (vaddq_u32 (e, ExponentBias)); diff --git a/pl/math/v_expm1_2u5.c b/pl/math/v_expm1_2u5.c index 0c2aa07..f867361 100644 --- a/pl/math/v_expm1_2u5.c +++ b/pl/math/v_expm1_2u5.c @@ -13,7 +13,7 @@ static const struct data { float64x2_t poly[11]; - float64x2_t invln2, ln2_lo, ln2_hi, shift; + float64x2_t invln2, ln2, shift; int64x2_t exponent_bias; #if WANT_SIMD_EXCEPT uint64x2_t thresh, tiny_bound; @@ -28,8 +28,7 @@ static const struct data V2 (0x1.71ddf82db5bb4p-19), V2 (0x1.27e517fc0d54bp-22), V2 (0x1.af5eedae67435p-26), V2 (0x1.1f143d060a28ap-29) }, .invln2 = V2 (0x1.71547652b82fep0), - .ln2_hi = V2 (0x1.62e42fefa39efp-1), - .ln2_lo = V2 (0x1.abc9e3b39803fp-56), + .ln2 = { 0x1.62e42fefa39efp-1, 0x1.abc9e3b39803fp-56 }, .shift = V2 (0x1.8p52), .exponent_bias = V2 (0x3ff0000000000000), #if WANT_SIMD_EXCEPT @@ -83,8 +82,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (expm1) (float64x2_t x) where 2^i is exact because i is an integer. */ float64x2_t n = vsubq_f64 (vfmaq_f64 (d->shift, d->invln2, x), d->shift); int64x2_t i = vcvtq_s64_f64 (n); - float64x2_t f = vfmsq_f64 (x, n, d->ln2_hi); - f = vfmsq_f64 (f, n, d->ln2_lo); + float64x2_t f = vfmsq_laneq_f64 (x, n, d->ln2, 0); + f = vfmsq_laneq_f64 (f, n, d->ln2, 1); /* Approximate expm1(f) using polynomial. Taylor expansion for expm1(x) has the form: diff --git a/pl/math/v_expm1f_1u6.c b/pl/math/v_expm1f_1u6.c index b755b30..ea6f255 100644 --- a/pl/math/v_expm1f_1u6.c +++ b/pl/math/v_expm1f_1u6.c @@ -13,7 +13,8 @@ static const struct data { float32x4_t poly[5]; - float32x4_t invln2, ln2_lo, ln2_hi, shift; + float32x4_t invln2_and_ln2; + float32x4_t shift; int32x4_t exponent_bias; #if WANT_SIMD_EXCEPT uint32x4_t thresh; @@ -24,9 +25,8 @@ static const struct data /* Generated using fpminimax with degree=5 in [-log(2)/2, log(2)/2]. */ .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5), V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) }, - .invln2 = V4 (0x1.715476p+0f), - .ln2_hi = V4 (0x1.62e4p-1f), - .ln2_lo = V4 (0x1.7f7d1cp-20f), + /* Stores constants: invln2, ln2_hi, ln2_lo, 0. */ + .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 }, .shift = V4 (0x1.8p23f), .exponent_bias = V4 (0x3f800000), #if !WANT_SIMD_EXCEPT @@ -78,10 +78,11 @@ float32x4_t VPCS_ATTR V_NAME_F1 (expm1) (float32x4_t x) and f = x - i * ln2, then f is in [-ln2/2, ln2/2]. exp(x) - 1 = 2^i * (expm1(f) + 1) - 1 where 2^i is exact because i is an integer. */ - float32x4_t j = vsubq_f32 (vfmaq_f32 (d->shift, d->invln2, x), d->shift); + float32x4_t j = vsubq_f32 ( + vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0), d->shift); int32x4_t i = vcvtq_s32_f32 (j); - float32x4_t f = vfmsq_f32 (x, j, d->ln2_hi); - f = vfmsq_f32 (f, j, d->ln2_lo); + float32x4_t f = vfmsq_laneq_f32 (x, j, d->invln2_and_ln2, 1); + f = vfmsq_laneq_f32 (f, j, d->invln2_and_ln2, 2); /* Approximate expm1(f) using polynomial. Taylor expansion for expm1(x) has the form: diff --git a/pl/math/v_expm1f_inline.h b/pl/math/v_expm1f_inline.h index 5a706f5..6ae94c4 100644 --- a/pl/math/v_expm1f_inline.h +++ b/pl/math/v_expm1f_inline.h @@ -16,20 +16,19 @@ struct v_expm1f_data { float32x4_t poly[5]; - float32x4_t invln2, ln2_lo, ln2_hi, shift; + float32x4_t invln2_and_ln2, shift; int32x4_t exponent_bias; }; /* Coefficients generated using fpminimax with degree=5 in [-log(2)/2, - * log(2)/2]. Exponent bias is asuint(1.0f). */ + log(2)/2]. Exponent bias is asuint(1.0f). + invln2_and_ln2 Stores constants: invln2, ln2_lo, ln2_hi, 0. */ #define V_EXPM1F_DATA \ { \ .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5), \ V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) }, \ - \ - .invln2 = V4 (0x1.715476p+0f), .ln2_hi = V4 (0x1.62e4p-1f), \ - .ln2_lo = V4 (0x1.7f7d1cp-20f), .shift = V4 (0x1.8p23f), \ - .exponent_bias = V4 (0x3f800000), \ + .shift = V4 (0x1.8p23f), .exponent_bias = V4 (0x3f800000), \ + .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 }, \ } static inline float32x4_t @@ -40,10 +39,11 @@ expm1f_inline (float32x4_t x, const struct v_expm1f_data *d) calling routine should handle special values if required. */ /* Reduce argument: f in [-ln2/2, ln2/2], i is exact. */ - float32x4_t j = vsubq_f32 (vfmaq_f32 (d->shift, d->invln2, x), d->shift); + float32x4_t j = vsubq_f32 ( + vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0), d->shift); int32x4_t i = vcvtq_s32_f32 (j); - float32x4_t f = vfmsq_f32 (x, j, d->ln2_hi); - f = vfmsq_f32 (f, j, d->ln2_lo); + float32x4_t f = vfmsq_laneq_f32 (x, j, d->invln2_and_ln2, 1); + f = vfmsq_laneq_f32 (f, j, d->invln2_and_ln2, 2); /* Approximate expm1(f) with polynomial P, expm1(f) ~= f + f^2 * P(f). Uses Estrin scheme, where the main _ZGVnN4v_expm1f routine uses diff --git a/pl/math/v_sinh_3u.c b/pl/math/v_sinh_3u.c index b99d6b9..393efa7 100644 --- a/pl/math/v_sinh_3u.c +++ b/pl/math/v_sinh_3u.c @@ -13,7 +13,7 @@ static const struct data { float64x2_t poly[11]; - float64x2_t inv_ln2, m_ln2_hi, m_ln2_lo, shift; + float64x2_t inv_ln2, m_ln2, shift; uint64x2_t halff; int64x2_t onef; #if WANT_SIMD_EXCEPT @@ -30,8 +30,7 @@ static const struct data V2 (0x1.af5eedae67435p-26), V2 (0x1.1f143d060a28ap-29), }, .inv_ln2 = V2 (0x1.71547652b82fep0), - .m_ln2_hi = V2 (-0x1.62e42fefa39efp-1), - .m_ln2_lo = V2 (-0x1.abc9e3b39803fp-56), + .m_ln2 = (float64x2_t) {-0x1.62e42fefa39efp-1, -0x1.abc9e3b39803fp-56}, .shift = V2 (0x1.8p52), .halff = V2 (0x3fe0000000000000), @@ -58,8 +57,8 @@ expm1_inline (float64x2_t x) and f = x - i * ln2 (f in [-ln2/2, ln2/2]). */ float64x2_t j = vsubq_f64 (vfmaq_f64 (d->shift, d->inv_ln2, x), d->shift); int64x2_t i = vcvtq_s64_f64 (j); - float64x2_t f = vfmaq_f64 (x, j, d->m_ln2_hi); - f = vfmaq_f64 (f, j, d->m_ln2_lo); + float64x2_t f = vfmaq_laneq_f64 (x, j, d->m_ln2, 0); + f = vfmaq_laneq_f64 (f, j, d->m_ln2, 1); /* Approximate expm1(f) using polynomial. */ float64x2_t f2 = vmulq_f64 (f, f); float64x2_t f4 = vmulq_f64 (f2, f2); diff --git a/pl/math/v_tan_3u5.c b/pl/math/v_tan_3u5.c index 90612f1..b390c92 100644 --- a/pl/math/v_tan_3u5.c +++ b/pl/math/v_tan_3u5.c @@ -13,7 +13,7 @@ static const struct data { float64x2_t poly[9]; - float64x2_t half_pi_hi, half_pi_lo, two_over_pi, shift; + float64x2_t half_pi, two_over_pi, shift; #if !WANT_SIMD_EXCEPT float64x2_t range_val; #endif @@ -24,8 +24,7 @@ static const struct data V2 (0x1.226e5e5ecdfa3p-7), V2 (0x1.d6c7ddbf87047p-9), V2 (0x1.7ea75d05b583ep-10), V2 (0x1.289f22964a03cp-11), V2 (0x1.4e4fd14147622p-12) }, - .half_pi_hi = V2 (0x1.921fb54442d18p0), - .half_pi_lo = V2 (0x1.1a62633145c07p-54), + .half_pi = { 0x1.921fb54442d18p0, 0x1.1a62633145c07p-54 }, .two_over_pi = V2 (0x1.45f306dc9c883p-1), .shift = V2 (0x1.8p52), #if !WANT_SIMD_EXCEPT @@ -51,10 +50,10 @@ special_case (float64x2_t x) float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x) { const struct data *dat = ptr_barrier (&data); - /* Our argument reduction cannot calculate q with sufficient accuracy for very - large inputs. Fall back to scalar routine for all lanes if any are too - large, or Inf/NaN. If fenv exceptions are expected, also fall back for tiny - input to avoid underflow. */ + /* Our argument reduction cannot calculate q with sufficient accuracy for + very large inputs. Fall back to scalar routine for all lanes if any are + too large, or Inf/NaN. If fenv exceptions are expected, also fall back for + tiny input to avoid underflow. */ #if WANT_SIMD_EXCEPT uint64x2_t iax = vreinterpretq_u64_f64 (vabsq_f64 (x)); /* iax - tiny_bound > range_val - tiny_bound. */ @@ -72,8 +71,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x) /* Use q to reduce x to r in [-pi/4, pi/4], by: r = x - q * pi/2, in extended precision. */ float64x2_t r = x; - r = vfmsq_f64 (r, q, dat->half_pi_hi); - r = vfmsq_f64 (r, q, dat->half_pi_lo); + r = vfmsq_laneq_f64 (r, q, dat->half_pi, 0); + r = vfmsq_laneq_f64 (r, q, dat->half_pi, 1); /* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle formula. */ r = vmulq_n_f64 (r, 0.5); @@ -96,7 +95,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x) and reciprocity around pi/2: tan(x) = 1 / (tan(pi/2 - x)) to assemble result using change-of-sign and conditional selection of - numerator/denominator, dependent on odd/even-ness of q (hence quadrant). */ + numerator/denominator, dependent on odd/even-ness of q (hence quadrant). + */ float64x2_t n = vfmaq_f64 (v_f64 (-1), p, p); float64x2_t d = vaddq_f64 (p, p); diff --git a/pl/math/v_tanf_3u5.c b/pl/math/v_tanf_3u5.c index e1f0a30..c3dd616 100644 --- a/pl/math/v_tanf_3u5.c +++ b/pl/math/v_tanf_3u5.c @@ -13,7 +13,8 @@ static const struct data { float32x4_t poly[6]; - float32x4_t neg_half_pi_1, neg_half_pi_2, neg_half_pi_3, two_over_pi, shift; + float32x4_t pi_consts; + float32x4_t shift; #if !WANT_SIMD_EXCEPT float32x4_t range_val; #endif @@ -21,10 +22,9 @@ static const struct data /* Coefficients generated using FPMinimax. */ .poly = { V4 (0x1.55555p-2f), V4 (0x1.11166p-3f), V4 (0x1.b88a78p-5f), V4 (0x1.7b5756p-6f), V4 (0x1.4ef4cep-8f), V4 (0x1.0e1e74p-7f) }, - .neg_half_pi_1 = V4 (-0x1.921fb6p+0f), - .neg_half_pi_2 = V4 (0x1.777a5cp-25f), - .neg_half_pi_3 = V4 (0x1.ee59dap-50f), - .two_over_pi = V4 (0x1.45f306p-1f), + /* Stores constants: (-pi/2)_high, (-pi/2)_mid, (-pi/2)_low, and 2/pi. */ + .pi_consts + = { -0x1.921fb6p+0f, 0x1.777a5cp-25f, 0x1.ee59dap-50f, 0x1.45f306p-1f }, .shift = V4 (0x1.8p+23f), #if !WANT_SIMD_EXCEPT .range_val = V4 (0x1p15f), @@ -48,10 +48,11 @@ eval_poly (float32x4_t z, const struct data *d) { float32x4_t z2 = vmulq_f32 (z, z); #if WANT_SIMD_EXCEPT - /* Tiny z (<= 0x1p-31) will underflow when calculating z^4. If fp exceptions - are to be triggered correctly, sidestep this by fixing such lanes to 0. */ + /* Tiny z (<= 0x1p-31) will underflow when calculating z^4. + If fp exceptions are to be triggered correctly, + sidestep this by fixing such lanes to 0. */ uint32x4_t will_uflow - = vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z)), TinyBound); + = vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z)), TinyBound); if (unlikely (v_any_u32 (will_uflow))) z2 = vbslq_f32 (will_uflow, v_f32 (0), z2); #endif @@ -84,16 +85,16 @@ float32x4_t VPCS_ATTR V_NAME_F1 (tan) (float32x4_t x) #endif /* n = rint(x/(pi/2)). */ - float32x4_t q = vfmaq_f32 (d->shift, d->two_over_pi, x); + float32x4_t q = vfmaq_laneq_f32 (d->shift, x, d->pi_consts, 3); float32x4_t n = vsubq_f32 (q, d->shift); /* Determine if x lives in an interval, where |tan(x)| grows to infinity. */ uint32x4_t pred_alt = vtstq_u32 (vreinterpretq_u32_f32 (q), v_u32 (1)); /* r = x - n * (pi/2) (range reduction into -pi./4 .. pi/4). */ float32x4_t r; - r = vfmaq_f32 (x, d->neg_half_pi_1, n); - r = vfmaq_f32 (r, d->neg_half_pi_2, n); - r = vfmaq_f32 (r, d->neg_half_pi_3, n); + r = vfmaq_laneq_f32 (x, n, d->pi_consts, 0); + r = vfmaq_laneq_f32 (r, n, d->pi_consts, 1); + r = vfmaq_laneq_f32 (r, n, d->pi_consts, 2); /* If x lives in an interval, where |tan(x)| - is finite, then use a polynomial approximation of the form -- cgit v1.2.3