diff options
Diffstat (limited to 'sysdeps/aarch64/fpu/sinh_sve.c')
| -rw-r--r-- | sysdeps/aarch64/fpu/sinh_sve.c | 165 |
1 files changed, 114 insertions, 51 deletions
diff --git a/sysdeps/aarch64/fpu/sinh_sve.c b/sysdeps/aarch64/fpu/sinh_sve.c index 963453f..072ba8f 100644 --- a/sysdeps/aarch64/fpu/sinh_sve.c +++ b/sysdeps/aarch64/fpu/sinh_sve.c @@ -18,90 +18,153 @@ <https://www.gnu.org/licenses/>. */ #include "sv_math.h" -#include "poly_sve_f64.h" static const struct data { - float64_t poly[11]; - float64_t inv_ln2, m_ln2_hi, m_ln2_lo, shift; uint64_t halff; - int64_t onef; - uint64_t large_bound; + double c2, c4; + double inv_ln2; + double ln2_hi, ln2_lo; + double c0, c1, c3; + double shift, special_bound, bound; + uint64_t expm1_data[20]; } data = { - /* Generated using Remez, deg=12 in [-log(2)/2, log(2)/2]. */ - .poly = { 0x1p-1, 0x1.5555555555559p-3, 0x1.555555555554bp-5, - 0x1.111111110f663p-7, 0x1.6c16c16c1b5f3p-10, - 0x1.a01a01affa35dp-13, 0x1.a01a018b4ecbbp-16, - 0x1.71ddf82db5bb4p-19, 0x1.27e517fc0d54bp-22, - 0x1.af5eedae67435p-26, 0x1.1f143d060a28ap-29, }, - - .inv_ln2 = 0x1.71547652b82fep0, - .m_ln2_hi = -0x1.62e42fefa39efp-1, - .m_ln2_lo = -0x1.abc9e3b39803fp-56, - .shift = 0x1.8p52, - + /* Table lookup of 2^(i/64) - 1, for values of i from 0..19. */ + .expm1_data = { + 0x0000000000000000, 0x3f864d1f3bc03077, 0x3f966c34c5615d0f, 0x3fa0e8a30eb37901, + 0x3fa6ab0d9f3121ec, 0x3fac7d865a7a3440, 0x3fb1301d0125b50a, 0x3fb429aaea92ddfb, + 0x3fb72b83c7d517ae, 0x3fba35beb6fcb754, 0x3fbd4873168b9aa8, 0x3fc031dc431466b2, + 0x3fc1c3d373ab11c3, 0x3fc35a2b2f13e6e9, 0x3fc4f4efa8fef709, 0x3fc6942d3720185a, + 0x3fc837f0518db8a9, 0x3fc9e0459320b7fa, 0x3fcb8d39b9d54e55, 0x3fcd3ed9a72cffb7, + }, + + /* Generated using Remez, in [-log(2)/128, log(2)/128]. */ + .c0 = 0x1p-1, + .c1 = 0x1.55555555548f9p-3, + .c2 = 0x1.5555555554c22p-5, + .c3 = 0x1.111123aaa2fb2p-7, + .c4 = 0x1.6c16d77d98e5bp-10, + .ln2_hi = 0x1.62e42fefa3800p-1, + .ln2_lo = 0x1.ef35793c76730p-45, + .inv_ln2 = 0x1.71547652b82fep+0, + .shift = 0x1.800000000ffc0p+46, /* 1.5*2^46+1023. */ .halff = 0x3fe0000000000000, - .onef = 0x3ff0000000000000, - /* 2^9. expm1 helper overflows for large input. */ - .large_bound = 0x4080000000000000, + .special_bound = 0x1.62e37e7d8ba72p+9, /* ln(2^(1024 - 1/128)). */ + .bound = 0x1.a56ef8ec924ccp-3 /* 19*ln2/64. */ }; +/* A specialised FEXPA expm1 that is only valid for positive inputs and + has no special cases. Based off the full FEXPA expm1 implementated for + _ZGVsMxv_expm1, with a slightly modified file to keep sinh under 3.5ULP. */ static inline svfloat64_t -expm1_inline (svfloat64_t x, svbool_t pg) +expm1_inline (svbool_t pg, svfloat64_t x) { const struct data *d = ptr_barrier (&data); - /* Reduce argument: - exp(x) - 1 = 2^i * (expm1(f) + 1) - 1 - where i = round(x / ln2) - and f = x - i * ln2 (f in [-ln2/2, ln2/2]). */ - svfloat64_t j - = svsub_x (pg, svmla_x (pg, sv_f64 (d->shift), x, d->inv_ln2), d->shift); - svint64_t i = svcvt_s64_x (pg, j); - svfloat64_t f = svmla_x (pg, x, j, d->m_ln2_hi); - f = svmla_x (pg, f, j, d->m_ln2_lo); - /* Approximate expm1(f) using polynomial. */ - svfloat64_t f2 = svmul_x (pg, f, f); - svfloat64_t f4 = svmul_x (pg, f2, f2); - svfloat64_t f8 = svmul_x (pg, f4, f4); - svfloat64_t p - = svmla_x (pg, f, f2, sv_estrin_10_f64_x (pg, f, f2, f4, f8, d->poly)); - /* t = 2^i. */ - svfloat64_t t = svscale_x (pg, sv_f64 (1), i); - /* expm1(x) ~= p * t + (t - 1). */ - return svmla_x (pg, svsub_x (pg, t, 1.0), p, t); + 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 ln2 = svld1rq (svptrue_b64 (), &d->ln2_hi); + svfloat64_t c24 = svld1rq (svptrue_b64 (), &d->c2); + + svfloat64_t r = x; + r = svmls_lane (r, n, ln2, 0); + r = svmls_lane (r, n, ln2, 1); + + svfloat64_t r2 = svmul_x (svptrue_b64 (), r, r); + + svfloat64_t p; + svfloat64_t c12 = svmla_lane (sv_f64 (d->c1), r, c24, 0); + svfloat64_t c34 = svmla_lane (sv_f64 (d->c3), r, c24, 1); + p = svmad_x (pg, c34, r2, c12); + p = svmad_x (pg, p, r, sv_f64 (d->c0)); + p = svmad_x (pg, p, r2, r); + + svfloat64_t scale = svexpa (u); + + /* We want to construct expm1(x) = (scale - 1) + scale * poly. + However, for values of scale close to 1, scale-1 causes large ULP errors + due to cancellation. + + This can be circumvented by using a small lookup for scale-1 + when our input is below a certain bound, otherwise we can use FEXPA. */ + svbool_t is_small = svaclt (pg, x, d->bound); + + /* Index via the input of FEXPA, but we only care about the lower 5 bits. */ + svuint64_t base_idx = svand_x (pg, u, 0x1f); + + /* Compute scale - 1 from FEXPA, and lookup values where this fails. */ + svfloat64_t scalem1_estimate = svsub_x (pg, scale, sv_f64 (1.0)); + svuint64_t scalem1_lookup + = svld1_gather_index (is_small, d->expm1_data, base_idx); + + /* Select the appropriate scale - 1 value based on x. */ + svfloat64_t scalem1 + = svsel (is_small, svreinterpret_f64 (scalem1_lookup), scalem1_estimate); + + /* return expm1 = scale - 1 + (scale * poly). */ + return svmla_x (pg, scalem1, scale, p); } +/* Vectorised special case to handle values past where exp_inline overflows. + Halves the input value and uses the identity exp(x) = exp(x/2)^2 to double + the valid range of inputs, and returns inf for anything past that. */ static svfloat64_t NOINLINE -special_case (svfloat64_t x, svbool_t pg) +special_case (svbool_t pg, svbool_t special, svfloat64_t ax, + svfloat64_t halfsign, const struct data *d) { - return sv_call_f64 (sinh, x, x, pg); + /* Halves input value, and then check if any cases + are still going to overflow. */ + ax = svmul_x (special, ax, 0.5); + svbool_t is_safe = svaclt (special, ax, d->special_bound); + + svfloat64_t t = expm1_inline (pg, ax); + + /* Finish fastpass to compute values for non-special cases. */ + svfloat64_t y = svadd_x (pg, t, svdiv_x (pg, t, svadd_x (pg, t, 1.0))); + y = svmul_x (pg, y, halfsign); + + /* Computes special lane, and set remaining overflow lanes to inf. */ + svfloat64_t half_special_y = svmul_x (svptrue_b64 (), t, halfsign); + svfloat64_t special_y = svmul_x (svptrue_b64 (), half_special_y, t); + + svuint64_t signed_inf + = svorr_x (svptrue_b64 (), svreinterpret_u64 (halfsign), + sv_u64 (0x7ff0000000000000)); + special_y = svsel (is_safe, special_y, svreinterpret_f64 (signed_inf)); + + /* Join resulting vectors together and return. */ + return svsel (special, special_y, y); } -/* Approximation for SVE double-precision sinh(x) using expm1. - sinh(x) = (exp(x) - exp(-x)) / 2. - The greatest observed error is 2.57 ULP: - _ZGVsMxv_sinh (0x1.a008538399931p-2) got 0x1.ab929fc64bd66p-2 - want 0x1.ab929fc64bd63p-2. */ +/* Approximation for SVE double-precision sinh(x) using FEXPA expm1. + Uses sinh(x) = e^2x - 1 / 2e^x, rewritten for accuracy. + The greatest observed error in the non-special region is 2.63 + 0.5 ULP: + _ZGVsMxv_sinh (0x1.b5e0e13ba88aep-2) got 0x1.c3587faf97b0cp-2 + want 0x1.c3587faf97b09p-2 + + The greatest observed error in the special region is 2.65 + 0.5 ULP: + _ZGVsMxv_sinh (0x1.633ce847dab1ap+9) got 0x1.fffd30eea0066p+1023 + want 0x1.fffd30eea0063p+1023. */ svfloat64_t SV_NAME_D1 (sinh) (svfloat64_t x, svbool_t pg) { const struct data *d = ptr_barrier (&data); + svbool_t special = svacge (pg, x, d->special_bound); svfloat64_t ax = svabs_x (pg, x); svuint64_t sign = sveor_x (pg, svreinterpret_u64 (x), svreinterpret_u64 (ax)); svfloat64_t halfsign = svreinterpret_f64 (svorr_x (pg, sign, d->halff)); - svbool_t special = svcmpge (pg, svreinterpret_u64 (ax), d->large_bound); - /* Fall back to scalar variant for all lanes if any are special. */ if (__glibc_unlikely (svptest_any (pg, special))) - return special_case (x, pg); + return special_case (pg, special, ax, halfsign, d); /* Up to the point that expm1 overflows, we can use it to calculate sinh using a slight rearrangement of the definition of sinh. This allows us to retain acceptable accuracy for very small inputs. */ - svfloat64_t t = expm1_inline (ax, pg); + svfloat64_t t = expm1_inline (pg, ax); t = svadd_x (pg, t, svdiv_x (pg, t, svadd_x (pg, t, 1.0))); return svmul_x (pg, t, halfsign); } |