diff options
| -rw-r--r-- | pl/math/include/mathlib.h | 2 | ||||
| -rw-r--r-- | pl/math/sv_erff_1u3.c | 91 | ||||
| -rw-r--r-- | pl/math/test/mathbench_funcs.h | 3 | ||||
| -rwxr-xr-x | pl/math/test/runulp.sh | 14 | ||||
| -rw-r--r-- | pl/math/test/ulp_funcs.h | 2 | ||||
| -rw-r--r-- | pl/math/test/ulp_wrappers.h | 1 |
6 files changed, 113 insertions, 0 deletions
diff --git a/pl/math/include/mathlib.h b/pl/math/include/mathlib.h index 4c906b4..6979b96 100644 --- a/pl/math/include/mathlib.h +++ b/pl/math/include/mathlib.h @@ -116,6 +116,7 @@ svfloat64_t __sv_atan_x (svfloat64_t, svbool_t); svfloat64_t __sv_atan2_x (svfloat64_t, svfloat64_t, svbool_t); svfloat32_t __sv_cosf_x (svfloat32_t, svbool_t); svfloat64_t __sv_cos_x (svfloat64_t, svbool_t); +svfloat32_t __sv_erff_x (svfloat32_t, svbool_t); svfloat32_t __sv_expf_x (svfloat32_t, svbool_t); svfloat32_t __sv_logf_x (svfloat32_t, svbool_t); svfloat64_t __sv_log_x (svfloat64_t, svbool_t); @@ -130,6 +131,7 @@ svfloat64_t _ZGVsMxv_atan (svfloat64_t, svbool_t); svfloat64_t _ZGVsMxvv_atan2 (svfloat64_t, svfloat64_t, svbool_t); svfloat32_t _ZGVsMxv_cosf (svfloat32_t, svbool_t); svfloat64_t _ZGVsMxv_cos (svfloat64_t, svbool_t); +svfloat32_t _ZGVsMxv_erff (svfloat32_t, svbool_t); svfloat32_t _ZGVsMxv_expf (svfloat32_t, svbool_t); svfloat32_t _ZGVsMxv_logf (svfloat32_t, svbool_t); svfloat64_t _ZGVsMxv_log (svfloat64_t, svbool_t); diff --git a/pl/math/sv_erff_1u3.c b/pl/math/sv_erff_1u3.c new file mode 100644 index 0000000..f0af98e --- /dev/null +++ b/pl/math/sv_erff_1u3.c @@ -0,0 +1,91 @@ +/* + * Single-precision vector erf(x) function. + * + * Copyright (c) 2020-2022, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "sv_math.h" +#if SV_SUPPORTED + +#define AbsMask (0x7fffffff) + +static NOINLINE sv_f32_t +__sv_erff_specialcase (sv_f32_t x, sv_f32_t y, svbool_t cmp) +{ + return sv_call_f32 (erff, x, y, cmp); +} + +sv_f32_t __sv_expf_x (svbool_t, sv_f32_t); + +/* Optimized single precision vector erf. Worst-case error is 1.25 ULP: + __sv_erff(0x1.dc59fap-1) got 0x1.9f9c88p-1 + want 0x1.9f9c8ap-1. */ +sv_f32_t +__sv_erff_x (sv_f32_t x, const svbool_t pg) +{ + sv_u32_t ix = sv_as_u32_f32 (x); + sv_u32_t atop = svand_n_u32_x (pg, svlsr_n_u32_x (pg, ix, 16), 0x7fff); + /* Handle both inf/nan as well as small values (|x|<2^-28). */ + svbool_t cmp + = svcmpge_n_u32 (pg, svsub_n_u32_x (pg, atop, 0x3180), 0x7ff0 - 0x3180); + + sv_u32_t sign = svand_n_u32_x (pg, ix, ~AbsMask); + /* |x| < 0.921875. */ + svbool_t red = svaclt_n_f32 (pg, x, 0.921875f); + /* |x| > 4.0. */ + svbool_t bor = svacgt_n_f32 (pg, x, 4.0f); + + /* Load polynomial coefficients. */ + sv_u32_t idx_lo = svsel (red, sv_u32 (0), sv_u32 (1)); + sv_u32_t idx_hi = svadd_n_u32_x (pg, idx_lo, 2); + + const float *base = (float *) __v_erff_data.coeffs; + sv_f32_t c_2_5 = svld1rq (svptrue_b32 (), base + 2); + sv_f32_t c_6_9 = svld1rq (svptrue_b32 (), base + 6); + sv_f32_t c_10_13 = svld1rq (svptrue_b32 (), base + 10); + + /* Do not need to store elem 0 of __v_erff_data as it is not used. */ + sv_f32_t p1 = svtbl (c_2_5, idx_lo); + sv_f32_t p2 = svtbl (c_2_5, idx_hi); + sv_f32_t p3 = svtbl (c_6_9, idx_lo); + sv_f32_t p4 = svtbl (c_6_9, idx_hi); + sv_f32_t p5 = svtbl (c_10_13, idx_lo); + sv_f32_t p6 = svtbl (c_10_13, idx_hi); + + sv_f32_t a = svabs_f32_x (pg, x); + /* Square with merging mul - z is x^2 for reduced, |x| otherwise. */ + sv_f32_t z = svmul_f32_m (red, a, a); + + /* Evaluate polynomial on |x| or x^2. */ + sv_f32_t r = sv_fma_f32_x (pg, z, p6, p5); + r = sv_fma_f32_x (pg, z, r, p4); + r = sv_fma_f32_x (pg, z, r, p3); + r = sv_fma_f32_x (pg, z, r, p2); + r = sv_fma_f32_x (pg, z, r, p1); + /* Use merging svmad for last operation - apply first coefficient if not + reduced, otherwise r is propagated unchanged. This is because the reduced + polynomial has lower order than the non-reduced. */ + r = svmad_n_f32_m (svnot_b_z (pg, red), r, z, base[1]); + r = sv_fma_f32_x (pg, a, r, a); + + /* y = |x| + |x| * P(x^2) if |x| < 0.921875 + y = 1 - exp (-(|x| + |x| * P(|x|))) otherwise. */ + sv_f32_t y = __sv_expf_x (pg, svneg_f32_x (pg, r)); + y = svsel_f32 (red, r, svsubr_n_f32_x (pg, y, 1.0)); + + /* Boring domain (absolute value is required to get the sign of erf(-nan) + right). */ + y = svsel_f32 (bor, sv_f32 (1.0f), svabs_f32_x (pg, y)); + + /* y = erf(x) if x>0, -erf(-x) otherwise. */ + y = sv_as_f32_u32 (sveor_u32_x (pg, sv_as_u32_f32 (y), sign)); + + if (unlikely (svptest_any (pg, cmp))) + return __sv_erff_specialcase (x, y, cmp); + return y; +} + +strong_alias (__sv_erff_x, _ZGVsMxv_erff) + +#endif diff --git a/pl/math/test/mathbench_funcs.h b/pl/math/test/mathbench_funcs.h index 134ae01..dad6ad6 100644 --- a/pl/math/test/mathbench_funcs.h +++ b/pl/math/test/mathbench_funcs.h @@ -119,6 +119,9 @@ SVD (_ZGVsMxv_atan, -3.1, 3.1) {"__sv_atan2", 'd', 'n', -10.0, 10.0, {.svd = __sv_atan2_wrap}}, {"_ZGVsM2vv_atan2", 'd', 'n', -10.0, 10.0, {.svd = _Z_sv_atan2_wrap}}, +SVF (__sv_erff_x, -4.0, 4.0) +SVF (_ZGVsMxv_erff, -4.0, 4.0) + SVF (__sv_expf_x, -9.9, 9.9) SVF (_ZGVsMxv_expf, -9.9, 9.9) diff --git a/pl/math/test/runulp.sh b/pl/math/test/runulp.sh index 32fa32f..caa49d7 100755 --- a/pl/math/test/runulp.sh +++ b/pl/math/test/runulp.sh @@ -395,6 +395,17 @@ range_sve_expf=' -0x1p23 -inf 50000 ' +range_sve_erff=' + 0 0x1p-28 20000 + 0x1p-28 1 60000 + 1 0x1p28 60000 + 0x1p28 inf 20000 + -0 -0x1p-28 20000 + -0x1p-28 -1 60000 + -1 -0x1p28 60000 + -0x1p28 -inf 20000 +' + # error limits L_erfc=3.14 L_erfcf=0.26 @@ -427,6 +438,7 @@ L_sve_log10f=2.82 L_sve_logf=2.85 L_sve_log=1.68 L_sve_expf=1.46 +L_sve_erff=0.76 while read G F R do @@ -522,6 +534,8 @@ sve_logf __sv_logf $runsv sve_logf _ZGVsMxv_logf $runsv sve_expf __sv_expf $runsv sve_expf _ZGVsMxv_expf $runsv +sve_erff __sv_erff $runsv +sve_erff _ZGVsMxv_erff $runsv sve_cos __sv_cos $runsv sve_cos _ZGVsMxv_cos $runsv diff --git a/pl/math/test/ulp_funcs.h b/pl/math/test/ulp_funcs.h index 3aa9d3a..452950f 100644 --- a/pl/math/test/ulp_funcs.h +++ b/pl/math/test/ulp_funcs.h @@ -80,6 +80,8 @@ SVF1 (cos) ZSVF1 (cos) SVD1 (cos) ZSVD1 (cos) +SVF1 (erf) +ZSVF1 (erf) SVF1 (exp) ZSVF1 (exp) SVF1 (log) diff --git a/pl/math/test/ulp_wrappers.h b/pl/math/test/ulp_wrappers.h index 623f945..16647e4 100644 --- a/pl/math/test/ulp_wrappers.h +++ b/pl/math/test/ulp_wrappers.h @@ -90,6 +90,7 @@ ZVND1_WRAP(log2) ZSVNF2_WRAP(atan2) ZSVNF1_WRAP(atan) ZSVNF1_WRAP(cos) +ZSVNF1_WRAP(erf) ZSVNF1_WRAP(exp) ZSVNF1_WRAP(log) ZSVNF1_WRAP(log10) |