diff options
Diffstat (limited to 'contrib/arm-optimized-routines/math/aarch64/experimental/sve')
5 files changed, 461 insertions, 0 deletions
diff --git a/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinv_25u.c b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinv_25u.c new file mode 100644 index 000000000000..4de6d08ab80f --- /dev/null +++ b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinv_25u.c @@ -0,0 +1,156 @@ +/* + * Double-precision inverse error function (SVE variant). + * + * Copyright (c) 2024, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ +#include "sv_math.h" +#include "test_defs.h" +#include "math_config.h" +#include "test_sig.h" +#include "sv_poly_f64.h" +#define SV_LOG_INLINE_POLY_ORDER 4 +#include "sv_log_inline.h" + +const static struct data +{ + /* We use P_N and Q_N to refer to arrays of coefficients, where P_N is the + coeffs of the numerator in table N of Blair et al, and Q_N is the coeffs + of the denominator. P is interleaved P_17 and P_37, similar for Q. */ + double P[7][2], Q[7][2]; + double P_57[9], Q_57[9], tailshift, P37_0; + struct sv_log_inline_data log_tbl; +} data = { + .P37_0 = -0x1.f3596123109edp-7, + .tailshift = -0.87890625, + .P = { { 0x1.007ce8f01b2e8p+4, 0x1.60b8fe375999ep-2 }, + { -0x1.6b23cc5c6c6d7p+6, -0x1.779bb9bef7c0fp+1 }, + { 0x1.74e5f6ceb3548p+7, 0x1.786ea384470a2p+3 }, + { -0x1.5200bb15cc6bbp+7, -0x1.6a7c1453c85d3p+4 }, + { 0x1.05d193233a849p+6, 0x1.31f0fc5613142p+4 }, + { -0x1.148c5474ee5e1p+3, -0x1.5ea6c007d4dbbp+2 }, + { 0x1.689181bbafd0cp-3, 0x1.e66f265ce9e5p-3 } }, + .Q = { { 0x1.d8fb0f913bd7bp+3, -0x1.636b2dcf4edbep-7 }, + { -0x1.6d7f25a3f1c24p+6, 0x1.0b5411e2acf29p-2 }, + { 0x1.a450d8e7f4cbbp+7, -0x1.3413109467a0bp+1 }, + { -0x1.bc3480485857p+7, 0x1.563e8136c554ap+3 }, + { 0x1.ae6b0c504ee02p+6, -0x1.7b77aab1dcafbp+4 }, + { -0x1.499dfec1a7f5fp+4, 0x1.8a3e174e05ddcp+4 }, + { 0x1p+0, -0x1.4075c56404eecp+3 } }, + .P_57 = { 0x1.b874f9516f7f1p-14, 0x1.5921f2916c1c4p-7, 0x1.145ae7d5b8fa4p-2, + 0x1.29d6dcc3b2fb7p+1, 0x1.cabe2209a7985p+2, 0x1.11859f0745c4p+3, + 0x1.b7ec7bc6a2ce5p+2, 0x1.d0419e0bb42aep+1, 0x1.c5aa03eef7258p-1 }, + .Q_57 = { 0x1.b8747e12691f1p-14, 0x1.59240d8ed1e0ap-7, 0x1.14aef2b181e2p-2, + 0x1.2cd181bcea52p+1, 0x1.e6e63e0b7aa4cp+2, 0x1.65cf8da94aa3ap+3, + 0x1.7e5c787b10a36p+3, 0x1.0626d68b6cea3p+3, 0x1.065c5f193abf6p+2 }, + .log_tbl = SV_LOG_CONSTANTS +}; + +static inline svfloat64_t +special (svbool_t pg, svfloat64_t x, const struct data *d) +{ + /* Note erfinv(inf) should return NaN, and erfinv(1) should return Inf. + By using log here, instead of log1p, we return finite values for both + these inputs, and values outside [-1, 1]. This is non-compliant, but is an + acceptable optimisation at Ofast. To get correct behaviour for all finite + values use the log1p_inline helper on -abs(x) - note that erfinv(inf) + will still be finite. */ + svfloat64_t ax = svabs_x (pg, x); + svfloat64_t t + = svneg_x (pg, sv_log_inline (pg, svsubr_x (pg, ax, 1), &d->log_tbl)); + t = svdivr_x (pg, svsqrt_x (pg, t), 1); + svuint64_t sign + = sveor_x (pg, svreinterpret_u64 (ax), svreinterpret_u64 (x)); + svfloat64_t ts + = svreinterpret_f64 (svorr_x (pg, sign, svreinterpret_u64 (t))); + + svfloat64_t q = svadd_x (pg, t, d->Q_57[8]); + for (int i = 7; i >= 0; i--) + q = svmad_x (pg, q, t, d->Q_57[i]); + + return svdiv_x (pg, sv_horner_8_f64_x (pg, t, d->P_57), svmul_x (pg, ts, q)); +} + +static inline svfloat64_t +lookup (const double *c, svuint64_t idx) +{ + svfloat64_t x = svld1rq_f64 (svptrue_b64 (), c); + return svtbl (x, idx); +} + +static inline svfloat64_t +notails (svbool_t pg, svfloat64_t x, const struct data *d) +{ + svfloat64_t t = svmad_x (pg, x, x, -0.5625); + svfloat64_t p = svmla_x (pg, sv_f64 (d->P[5][0]), t, d->P[6][0]); + svfloat64_t q = svadd_x (pg, t, d->Q[5][0]); + for (int i = 4; i >= 0; i--) + { + p = svmad_x (pg, t, p, d->P[i][0]); + q = svmad_x (pg, t, q, d->Q[i][0]); + } + p = svmul_x (pg, p, x); + return svdiv_x (pg, p, q); +} + +/* Vector implementation of Blair et al's rational approximation to inverse + error function in double precision. Largest observed error is 24.75 ULP: + _ZGVsMxv_erfinv(0x1.fc861d81c2ba8p-1) got 0x1.ea05472686625p+0 + want 0x1.ea0547268660cp+0. */ +svfloat64_t SV_NAME_D1 (erfinv) (svfloat64_t x, svbool_t pg) +{ + const struct data *d = ptr_barrier (&data); + /* Calculate inverse error using algorithm described in + J. M. Blair, C. A. Edwards, and J. H. Johnson, + "Rational Chebyshev approximations for the inverse of the error function", + Math. Comp. 30, pp. 827--830 (1976). + https://doi.org/10.1090/S0025-5718-1976-0421040-7. + + Algorithm has 3 intervals: + - 'Normal' region [-0.75, 0.75] + - Tail region [0.75, 0.9375] U [-0.9375, -0.75] + - Extreme tail [-1, -0.9375] U [0.9375, 1] + Normal and tail are both rational approximation of similar order on + shifted input - these are typically performed in parallel using gather + loads to obtain correct coefficients depending on interval. */ + + svbool_t no_tail = svacle (pg, x, 0.75); + if (unlikely (!svptest_any (pg, svnot_z (pg, no_tail)))) + return notails (pg, x, d); + + svbool_t is_tail = svnot_z (pg, no_tail); + svbool_t extreme_tail = svacgt (pg, x, 0.9375); + svuint64_t idx = svdup_n_u64_z (is_tail, 1); + + svfloat64_t t = svsel_f64 (is_tail, sv_f64 (d->tailshift), sv_f64 (-0.5625)); + t = svmla_x (pg, t, x, x); + + svfloat64_t p = lookup (&d->P[6][0], idx); + svfloat64_t q + = svmla_x (pg, lookup (&d->Q[6][0], idx), svdup_n_f64_z (is_tail, 1), t); + for (int i = 5; i >= 0; i--) + { + p = svmla_x (pg, lookup (&d->P[i][0], idx), p, t); + q = svmla_x (pg, lookup (&d->Q[i][0], idx), q, t); + } + p = svmad_m (is_tail, p, t, d->P37_0); + p = svmul_x (pg, p, x); + + if (likely (svptest_any (pg, extreme_tail))) + return svsel (extreme_tail, special (pg, x, d), svdiv_x (pg, p, q)); + return svdiv_x (pg, p, q); +} + +#if USE_MPFR +# warning Not generating tests for _ZGVsMxv_erfinv, as MPFR has no suitable reference +#else +TEST_SIG (SV, D, 1, erfinv, -0.99, 0.99) +TEST_ULP (SV_NAME_D1 (erfinv), 24.5) +TEST_DISABLE_FENV (SV_NAME_D1 (erfinv)) +/* Test with control lane in each interval. */ +TEST_SYM_INTERVAL (SV_NAME_F1 (erfinv), 0, 1, 100000) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.5) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.8) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.95) +#endif +CLOSE_SVE_ATTR diff --git a/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinvf_5u.c b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinvf_5u.c new file mode 100644 index 000000000000..2c81c4e0b9a2 --- /dev/null +++ b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/erfinvf_5u.c @@ -0,0 +1,156 @@ +/* + * Single-precision inverse error function (SVE variant). + * + * Copyright (c) 2024, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ +#include "sv_math.h" +#include "test_sig.h" +#include "test_defs.h" +#include "sv_poly_f32.h" +#include "sv_logf_inline.h" + +const static struct data +{ + /* We use P_N and Q_N to refer to arrays of coefficients, where P_N + is the coeffs of the numerator in table N of Blair et al, and + Q_N is the coeffs of the denominator. Coefficients stored in + interleaved format to support lookup scheme. */ + float P10_2, P29_3, Q10_2, Q29_2; + float P10_0, P29_1, P10_1, P29_2; + float Q10_0, Q29_0, Q10_1, Q29_1; + float P29_0, P_50[6], Q_50[2], tailshift; + struct sv_logf_data logf_tbl; +} data = { .P10_0 = -0x1.a31268p+3, + .P10_1 = 0x1.ac9048p+4, + .P10_2 = -0x1.293ff6p+3, + .P29_0 = -0x1.fc0252p-4, + .P29_1 = 0x1.119d44p+0, + .P29_2 = -0x1.f59ee2p+0, + .P29_3 = 0x1.b13626p-2, + .Q10_0 = -0x1.8265eep+3, + .Q10_1 = 0x1.ef5eaep+4, + .Q10_2 = -0x1.12665p+4, + .Q29_0 = -0x1.69952p-4, + .Q29_1 = 0x1.c7b7d2p-1, + .Q29_2 = -0x1.167d7p+1, + .P_50 = { 0x1.3d8948p-3, 0x1.61f9eap+0, 0x1.61c6bcp-1, + -0x1.20c9f2p+0, 0x1.5c704cp-1, -0x1.50c6bep-3 }, + .Q_50 = { 0x1.3d7dacp-3, 0x1.629e5p+0 }, + .tailshift = -0.87890625, + .logf_tbl = SV_LOGF_CONSTANTS }; + +static inline svfloat32_t +special (svbool_t pg, svfloat32_t x, const struct data *d) +{ + svfloat32_t ax = svabs_x (pg, x); + svfloat32_t t = svdivr_x ( + pg, + svsqrt_x (pg, svneg_x (pg, sv_logf_inline (pg, svsubr_x (pg, ax, 1), + &d->logf_tbl))), + 1); + svuint32_t sign + = sveor_x (pg, svreinterpret_u32 (ax), svreinterpret_u32 (x)); + svfloat32_t ts + = svreinterpret_f32 (svorr_x (pg, sign, svreinterpret_u32 (t))); + svfloat32_t q + = svmla_x (pg, sv_f32 (d->Q_50[0]), svadd_x (pg, t, d->Q_50[1]), t); + return svdiv_x (pg, sv_horner_5_f32_x (pg, t, d->P_50), svmul_x (pg, ts, q)); +} + +static inline svfloat32_t +notails (svbool_t pg, svfloat32_t x, const struct data *d) +{ + /* Shortcut when no input is in a tail region - no need to gather shift or + coefficients. */ + svfloat32_t t = svmad_x (pg, x, x, -0.5625); + svfloat32_t q = svadd_x (pg, t, d->Q10_2); + q = svmad_x (pg, t, q, d->Q10_1); + q = svmad_x (pg, t, q, d->Q10_0); + + svfloat32_t p = svmla_x (pg, sv_f32 (d->P10_1), t, d->P10_2); + p = svmad_x (pg, p, t, d->P10_0); + + return svdiv_x (pg, svmul_x (pg, x, p), q); +} + +/* Vector implementation of Blair et al's rational approximation to inverse + error function in single-precision. Worst-case error is 4.71 ULP, in the + tail region: + _ZGVsMxv_erfinvf(0x1.f84e9ap-1) got 0x1.b8326ap+0 + want 0x1.b83274p+0. */ +svfloat32_t SV_NAME_F1 (erfinv) (svfloat32_t x, svbool_t pg) +{ + const struct data *d = ptr_barrier (&data); + + /* Calculate inverse error using algorithm described in + J. M. Blair, C. A. Edwards, and J. H. Johnson, + "Rational Chebyshev approximations for the inverse of the error function", + Math. Comp. 30, pp. 827--830 (1976). + https://doi.org/10.1090/S0025-5718-1976-0421040-7. */ + + /* Algorithm has 3 intervals: + - 'Normal' region [-0.75, 0.75] + - Tail region [0.75, 0.9375] U [-0.9375, -0.75] + - Extreme tail [-1, -0.9375] U [0.9375, 1] + Normal and tail are both rational approximation of similar order on + shifted input - these are typically performed in parallel using gather + loads to obtain correct coefficients depending on interval. */ + svbool_t is_tail = svacge (pg, x, 0.75); + svbool_t extreme_tail = svacge (pg, x, 0.9375); + + if (likely (!svptest_any (pg, is_tail))) + return notails (pg, x, d); + + /* Select requisite shift depending on interval: polynomial is evaluated on + x * x - shift. + Normal shift = 0.5625 + Tail shift = 0.87890625. */ + svfloat32_t t = svmla_x ( + pg, svsel (is_tail, sv_f32 (d->tailshift), sv_f32 (-0.5625)), x, x); + + svuint32_t idx = svdup_u32_z (is_tail, 1); + svuint32_t idxhi = svadd_x (pg, idx, 2); + + /* Load coeffs in quadwords and select them according to interval. */ + svfloat32_t pqhi = svld1rq (svptrue_b32 (), &d->P10_2); + svfloat32_t plo = svld1rq (svptrue_b32 (), &d->P10_0); + svfloat32_t qlo = svld1rq (svptrue_b32 (), &d->Q10_0); + + svfloat32_t p2 = svtbl (pqhi, idx); + svfloat32_t p1 = svtbl (plo, idxhi); + svfloat32_t p0 = svtbl (plo, idx); + svfloat32_t q0 = svtbl (qlo, idx); + svfloat32_t q1 = svtbl (qlo, idxhi); + svfloat32_t q2 = svtbl (pqhi, idxhi); + + svfloat32_t p = svmla_x (pg, p1, p2, t); + p = svmla_x (pg, p0, p, t); + /* Tail polynomial has higher order - merge with normal lanes. */ + p = svmad_m (is_tail, p, t, d->P29_0); + svfloat32_t y = svmul_x (pg, x, p); + + /* Least significant term of both Q polynomials is 1, so no need to generate + it. */ + svfloat32_t q = svadd_x (pg, t, q2); + q = svmla_x (pg, q1, q, t); + q = svmla_x (pg, q0, q, t); + + if (unlikely (svptest_any (pg, extreme_tail))) + return svsel (extreme_tail, special (extreme_tail, x, d), + svdiv_x (pg, y, q)); + return svdiv_x (pg, y, q); +} + +#if USE_MPFR +# warning Not generating tests for _ZGVsMxv_erfinvf, as MPFR has no suitable reference +#else +TEST_SIG (SV, F, 1, erfinv, -0.99, 0.99) +TEST_ULP (SV_NAME_F1 (erfinv), 4.09) +TEST_DISABLE_FENV (SV_NAME_F1 (erfinv)) +TEST_SYM_INTERVAL (SV_NAME_F1 (erfinv), 0, 1, 40000) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.5) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.8) +TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.95) +#endif +CLOSE_SVE_ATTR diff --git a/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powi.c b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powi.c new file mode 100644 index 000000000000..62dd1b114970 --- /dev/null +++ b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powi.c @@ -0,0 +1,49 @@ +/* + * Double-precision SVE powi(x, n) function. + * + * Copyright (c) 2020-2024, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "sv_math.h" + +/* Optimized double-precision vector powi (double base, long integer power). + powi is developed for environments in which accuracy is of much less + importance than performance, hence we provide no estimate for worst-case + error. */ +svfloat64_t +_ZGVsMxvv_powk (svfloat64_t as, svint64_t ns, svbool_t p) +{ + /* Compute powi by successive squaring, right to left. */ + svfloat64_t acc = sv_f64 (1.0); + svbool_t want_recip = svcmplt (p, ns, 0); + svuint64_t ns_abs = svreinterpret_u64 (svabs_x (p, ns)); + + /* We use a max to avoid needing to check whether any lane != 0 on each + iteration. */ + uint64_t max_n = svmaxv (p, ns_abs); + + svfloat64_t c = as; + /* Successively square c, and use merging predication (_m) to determine + whether or not to perform the multiplication or keep the previous + iteration. */ + while (true) + { + svbool_t px = svcmpeq (p, svand_x (p, ns_abs, 1ull), 1ull); + acc = svmul_m (px, acc, c); + max_n >>= 1; + if (max_n == 0) + break; + + ns_abs = svlsr_x (p, ns_abs, 1); + c = svmul_x (p, c, c); + } + + /* Negative powers are handled by computing the abs(n) version and then + taking the reciprocal. */ + if (svptest_any (want_recip, want_recip)) + acc = svdivr_m (want_recip, acc, 1.0); + + return acc; +} +CLOSE_SVE_ATTR diff --git a/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powif.c b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powif.c new file mode 100644 index 000000000000..fd74acf12df7 --- /dev/null +++ b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/powif.c @@ -0,0 +1,49 @@ +/* + * Single-precision SVE powi(x, n) function. + * + * Copyright (c) 2020-2024, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "sv_math.h" + +/* Optimized single-precision vector powi (float base, integer power). + powi is developed for environments in which accuracy is of much less + importance than performance, hence we provide no estimate for worst-case + error. */ +svfloat32_t +_ZGVsMxvv_powi (svfloat32_t as, svint32_t ns, svbool_t p) +{ + /* Compute powi by successive squaring, right to left. */ + svfloat32_t acc = sv_f32 (1.f); + svbool_t want_recip = svcmplt (p, ns, 0); + svuint32_t ns_abs = svreinterpret_u32 (svabs_x (p, ns)); + + /* We use a max to avoid needing to check whether any lane != 0 on each + iteration. */ + uint32_t max_n = svmaxv (p, ns_abs); + + svfloat32_t c = as; + /* Successively square c, and use merging predication (_m) to determine + whether or not to perform the multiplication or keep the previous + iteration. */ + while (true) + { + svbool_t px = svcmpeq (p, svand_x (p, ns_abs, 1), 1); + acc = svmul_m (px, acc, c); + max_n >>= 1; + if (max_n == 0) + break; + + ns_abs = svlsr_x (p, ns_abs, 1); + c = svmul_x (p, c, c); + } + + /* Negative powers are handled by computing the abs(n) version and then + taking the reciprocal. */ + if (svptest_any (want_recip, want_recip)) + acc = svdivr_m (want_recip, acc, 1.0f); + + return acc; +} +CLOSE_SVE_ATTR diff --git a/contrib/arm-optimized-routines/math/aarch64/experimental/sve/sv_logf_inline.h b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/sv_logf_inline.h new file mode 100644 index 000000000000..c317a23f6fc3 --- /dev/null +++ b/contrib/arm-optimized-routines/math/aarch64/experimental/sve/sv_logf_inline.h @@ -0,0 +1,51 @@ +/* + * Single-precision vector log function - inline version + * + * Copyright (c) 2024, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "sv_math.h" + +struct sv_logf_data +{ + float p1, p3, p5, p6, p0, p2, p4; + float ln2; + uint32_t off, mantissa_mask; +}; + +#define SV_LOGF_CONSTANTS \ + { \ + .p0 = -0x1.ffffc8p-2f, .p1 = 0x1.555d7cp-2f, .p2 = -0x1.00187cp-2f, \ + .p3 = 0x1.961348p-3f, .p4 = -0x1.4f9934p-3f, .p5 = 0x1.5a9aa2p-3f, \ + .p6 = -0x1.3e737cp-3f, .ln2 = 0x1.62e43p-1f, .off = 0x3f2aaaab, \ + .mantissa_mask = 0x007fffff \ + } + +static inline svfloat32_t +sv_logf_inline (svbool_t pg, svfloat32_t x, const struct sv_logf_data *d) +{ + svuint32_t u = svreinterpret_u32 (x); + + /* x = 2^n * (1+r), where 2/3 < 1+r < 4/3. */ + u = svsub_x (pg, u, d->off); + svfloat32_t n = svcvt_f32_s32_x ( + pg, svasr_x (pg, svreinterpret_s32_u32 (u), 23)); /* signextend. */ + u = svand_x (pg, u, d->mantissa_mask); + u = svadd_x (pg, u, d->off); + svfloat32_t r = svsub_x (pg, svreinterpret_f32 (u), 1.0f); + + /* y = log(1+r) + n*ln2. */ + svfloat32_t r2 = svmul_x (pg, r, r); + /* n*ln2 + r + r2*(P1 + r*P2 + r2*(P3 + r*P4 + r2*(P5 + r*P6 + r2*P7))). */ + svfloat32_t p1356 = svld1rq_f32 (svptrue_b32 (), &d->p1); + svfloat32_t p = svmla_lane (sv_f32 (d->p4), r, p1356, 2); + svfloat32_t q = svmla_lane (sv_f32 (d->p2), r, p1356, 1); + svfloat32_t y = svmla_lane (sv_f32 (d->p0), r, p1356, 0); + p = svmla_lane (p, r2, p1356, 3); + q = svmla_x (pg, q, p, r2); + y = svmla_x (pg, y, q, r2); + p = svmla_x (pg, r, n, d->ln2); + + return svmla_x (pg, p, y, r2); +} |