diff options
| author | Joe Ramsay <Joe.Ramsay@arm.com> | 2024-02-20 16:59:40 +0000 |
|---|---|---|
| committer | Szabolcs Nagy <szabolcs.nagy@arm.com> | 2024-04-04 10:32:58 +0100 |
| commit | b09fee1d21650428a6a3335408a46ebe1165d30d (patch) | |
| tree | a1276ae6426feddebe06b59a2183c478eaf93057 | |
| parent | bdb5705b7bab618ed4445f4b17d4b1e4fbbf94a7 (diff) | |
| download | glibc-b09fee1d21650428a6a3335408a46ebe1165d30d.tar.xz glibc-b09fee1d21650428a6a3335408a46ebe1165d30d.zip | |
aarch64/fpu: Add vector variants of acosh
Reviewed-by: Szabolcs Nagy <szabolcs.nagy@arm.com>
| -rw-r--r-- | sysdeps/aarch64/fpu/Makefile | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/Versions | 5 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/acosh_advsimd.c | 67 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/acosh_sve.c | 51 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/acoshf_advsimd.c | 78 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/acoshf_sve.c | 49 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/advsimd_f32_protos.h | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/bits/math-vector.h | 8 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/sv_log1p_inline.h | 109 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/sv_log1pf_inline.h | 76 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/test-double-sve-wrappers.c | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/test-float-sve-wrappers.c | 1 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/v_log1p_inline.h | 103 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/v_log1pf_inline.h | 78 | ||||
| -rw-r--r-- | sysdeps/aarch64/fpu/v_math.h | 10 | ||||
| -rw-r--r-- | sysdeps/aarch64/libm-test-ulps | 8 | ||||
| -rw-r--r-- | sysdeps/unix/sysv/linux/aarch64/libmvec.abilist | 5 |
19 files changed, 653 insertions, 0 deletions
diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile index 019c3a5188..2e5bbb5a07 100644 --- a/sysdeps/aarch64/fpu/Makefile +++ b/sysdeps/aarch64/fpu/Makefile @@ -1,4 +1,5 @@ libmvec-supported-funcs = acos \ + acosh \ asin \ atan \ atan2 \ diff --git a/sysdeps/aarch64/fpu/Versions b/sysdeps/aarch64/fpu/Versions index 884b4b57f0..60e1cdeace 100644 --- a/sysdeps/aarch64/fpu/Versions +++ b/sysdeps/aarch64/fpu/Versions @@ -79,6 +79,11 @@ libmvec { _ZGVsMxv_tan; } GLIBC_2.40 { + _ZGVnN2v_acosh; + _ZGVnN2v_acoshf; + _ZGVnN4v_acoshf; + _ZGVsMxv_acosh; + _ZGVsMxv_acoshf; _ZGVnN2v_cosh; _ZGVnN2v_coshf; _ZGVnN4v_coshf; diff --git a/sysdeps/aarch64/fpu/acosh_advsimd.c b/sysdeps/aarch64/fpu/acosh_advsimd.c new file mode 100644 index 0000000000..c88283cf11 --- /dev/null +++ b/sysdeps/aarch64/fpu/acosh_advsimd.c @@ -0,0 +1,67 @@ +/* Double-precision vector (Advanced SIMD) acosh function + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#define WANT_V_LOG1P_K0_SHORTCUT 1 +#include "v_log1p_inline.h" + +const static struct data +{ + struct v_log1p_data log1p_consts; + uint64x2_t one, thresh; +} data = { + .log1p_consts = V_LOG1P_CONSTANTS_TABLE, + .one = V2 (0x3ff0000000000000), + .thresh = V2 (0x1ff0000000000000) /* asuint64(0x1p511) - asuint64(1). */ +}; + +static float64x2_t NOINLINE VPCS_ATTR +special_case (float64x2_t x, float64x2_t y, uint64x2_t special, + const struct v_log1p_data *d) +{ + return v_call_f64 (acosh, x, log1p_inline (y, d), special); +} + +/* Vector approximation for double-precision acosh, based on log1p. + The largest observed error is 3.02 ULP in the region where the + argument to log1p falls in the k=0 interval, i.e. x close to 1: + _ZGVnN2v_acosh(0x1.00798aaf80739p+0) got 0x1.f2d6d823bc9dfp-5 + want 0x1.f2d6d823bc9e2p-5. */ +VPCS_ATTR float64x2_t V_NAME_D1 (acosh) (float64x2_t x) +{ + const struct data *d = ptr_barrier (&data); + uint64x2_t special + = vcgeq_u64 (vsubq_u64 (vreinterpretq_u64_f64 (x), d->one), d->thresh); + float64x2_t special_arg = x; + +#if WANT_SIMD_EXCEPT + if (__glibc_unlikely (v_any_u64 (special))) + x = vbslq_f64 (special, vreinterpretq_f64_u64 (d->one), x); +#endif + + float64x2_t xm1 = vsubq_f64 (x, v_f64 (1)); + float64x2_t y; + y = vaddq_f64 (x, v_f64 (1)); + y = vmulq_f64 (y, xm1); + y = vsqrtq_f64 (y); + y = vaddq_f64 (xm1, y); + + if (__glibc_unlikely (v_any_u64 (special))) + return special_case (special_arg, y, special, &d->log1p_consts); + return log1p_inline (y, &d->log1p_consts); +} diff --git a/sysdeps/aarch64/fpu/acosh_sve.c b/sysdeps/aarch64/fpu/acosh_sve.c new file mode 100644 index 0000000000..3e4faaa5ca --- /dev/null +++ b/sysdeps/aarch64/fpu/acosh_sve.c @@ -0,0 +1,51 @@ +/* Double-precision vector (SVE) acosh function + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#define WANT_SV_LOG1P_K0_SHORTCUT 1 +#include "sv_log1p_inline.h" + +#define One (0x3ff0000000000000) +#define Thres (0x1ff0000000000000) /* asuint64 (0x1p511) - One. */ + +static svfloat64_t NOINLINE +special_case (svfloat64_t x, svfloat64_t y, svbool_t special) +{ + return sv_call_f64 (acosh, x, y, special); +} + +/* SVE approximation for double-precision acosh, based on log1p. + The largest observed error is 3.19 ULP in the region where the + argument to log1p falls in the k=0 interval, i.e. x close to 1: + SV_NAME_D1 (acosh)(0x1.1e4388d4ca821p+0) got 0x1.ed23399f5137p-2 + want 0x1.ed23399f51373p-2. */ +svfloat64_t SV_NAME_D1 (acosh) (svfloat64_t x, const svbool_t pg) +{ + /* (ix - One) >= (BigBound - One). */ + svuint64_t ix = svreinterpret_u64 (x); + svbool_t special = svcmpge (pg, svsub_x (pg, ix, One), Thres); + + svfloat64_t xm1 = svsub_x (pg, x, 1.0); + svfloat64_t u = svmul_x (pg, xm1, svadd_x (pg, x, 1.0)); + svfloat64_t y = svadd_x (pg, xm1, svsqrt_x (pg, u)); + + /* Fall back to scalar routine for special lanes. */ + if (__glibc_unlikely (svptest_any (pg, special))) + return special_case (x, sv_log1p_inline (y, pg), special); + return sv_log1p_inline (y, pg); +} diff --git a/sysdeps/aarch64/fpu/acoshf_advsimd.c b/sysdeps/aarch64/fpu/acoshf_advsimd.c new file mode 100644 index 0000000000..8916dcbf40 --- /dev/null +++ b/sysdeps/aarch64/fpu/acoshf_advsimd.c @@ -0,0 +1,78 @@ +/* Single-precision vector (Advanced SIMD) acosh function + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#include "v_log1pf_inline.h" + +#define SquareLim 0x1p64 + +const static struct data +{ + struct v_log1pf_data log1pf_consts; + uint32x4_t one; + uint16x4_t thresh; +} data = { + .log1pf_consts = V_LOG1PF_CONSTANTS_TABLE, + .one = V4 (0x3f800000), + .thresh = V4 (0x2000) /* top(asuint(SquareLim) - asuint(1)). */ +}; + +static float32x4_t NOINLINE VPCS_ATTR +special_case (float32x4_t x, float32x4_t y, uint16x4_t special, + const struct v_log1pf_data d) +{ + return v_call_f32 (acoshf, x, log1pf_inline (y, d), vmovl_u16 (special)); +} + +/* Vector approximation for single-precision acosh, based on log1p. Maximum + error depends on WANT_SIMD_EXCEPT. With SIMD fp exceptions enabled, it + is 2.78 ULP: + __v_acoshf(0x1.07887p+0) got 0x1.ef9e9cp-3 + want 0x1.ef9ea2p-3. + With exceptions disabled, we can compute u with a shorter dependency chain, + which gives maximum error of 3.07 ULP: + __v_acoshf(0x1.01f83ep+0) got 0x1.fbc7fap-4 + want 0x1.fbc7f4p-4. */ + +VPCS_ATTR float32x4_t NOINLINE V_NAME_F1 (acosh) (float32x4_t x) +{ + const struct data *d = ptr_barrier (&data); + uint32x4_t ix = vreinterpretq_u32_f32 (x); + uint16x4_t special = vcge_u16 (vsubhn_u32 (ix, d->one), d->thresh); + +#if WANT_SIMD_EXCEPT + /* Mask special lanes with 1 to side-step spurious invalid or overflow. Use + only xm1 to calculate u, as operating on x will trigger invalid for NaN. + Widening sign-extend special predicate in order to mask with it. */ + uint32x4_t p + = vreinterpretq_u32_s32 (vmovl_s16 (vreinterpret_s16_u16 (special))); + float32x4_t xm1 = v_zerofy_f32 (vsubq_f32 (x, v_f32 (1)), p); + float32x4_t u = vfmaq_f32 (vaddq_f32 (xm1, xm1), xm1, xm1); +#else + float32x4_t xm1 = vsubq_f32 (x, v_f32 (1)); + float32x4_t u = vmulq_f32 (xm1, vaddq_f32 (x, v_f32 (1.0f))); +#endif + + float32x4_t y = vaddq_f32 (xm1, vsqrtq_f32 (u)); + + if (__glibc_unlikely (v_any_u16h (special))) + return special_case (x, y, special, d->log1pf_consts); + return log1pf_inline (y, d->log1pf_consts); +} +libmvec_hidden_def (V_NAME_F1 (acosh)) +HALF_WIDTH_ALIAS_F1 (acosh) diff --git a/sysdeps/aarch64/fpu/acoshf_sve.c b/sysdeps/aarch64/fpu/acoshf_sve.c new file mode 100644 index 0000000000..2110894e62 --- /dev/null +++ b/sysdeps/aarch64/fpu/acoshf_sve.c @@ -0,0 +1,49 @@ +/* Single-precision vector (SVE) acosh function + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#define One 0x3f800000 +#define Thres 0x20000000 /* asuint(0x1p64) - One. */ + +#include "sv_log1pf_inline.h" + +static svfloat32_t NOINLINE +special_case (svfloat32_t x, svfloat32_t y, svbool_t special) +{ + return sv_call_f32 (acoshf, x, y, special); +} + +/* Single-precision SVE acosh(x) routine. Implements the same algorithm as + vector acoshf and log1p. + + Maximum error is 2.78 ULPs: + SV_NAME_F1 (acosh) (0x1.01e996p+0) got 0x1.f45b42p-4 + want 0x1.f45b3cp-4. */ +svfloat32_t SV_NAME_F1 (acosh) (svfloat32_t x, const svbool_t pg) +{ + svuint32_t ix = svreinterpret_u32 (x); + svbool_t special = svcmpge (pg, svsub_x (pg, ix, One), Thres); + + svfloat32_t xm1 = svsub_x (pg, x, 1.0f); + svfloat32_t u = svmul_x (pg, xm1, svadd_x (pg, x, 1.0f)); + svfloat32_t y = sv_log1pf_inline (svadd_x (pg, xm1, svsqrt_x (pg, u)), pg); + + if (__glibc_unlikely (svptest_any (pg, special))) + return special_case (x, y, special); + return y; +} diff --git a/sysdeps/aarch64/fpu/advsimd_f32_protos.h b/sysdeps/aarch64/fpu/advsimd_f32_protos.h index c63b2948d4..22fec4de77 100644 --- a/sysdeps/aarch64/fpu/advsimd_f32_protos.h +++ b/sysdeps/aarch64/fpu/advsimd_f32_protos.h @@ -18,6 +18,7 @@ <https://www.gnu.org/licenses/>. */ libmvec_hidden_proto (V_NAME_F1(acos)); +libmvec_hidden_proto (V_NAME_F1(acosh)); libmvec_hidden_proto (V_NAME_F1(asin)); libmvec_hidden_proto (V_NAME_F1(atan)); libmvec_hidden_proto (V_NAME_F1(cos)); diff --git a/sysdeps/aarch64/fpu/bits/math-vector.h b/sysdeps/aarch64/fpu/bits/math-vector.h index 8ca5509870..841330956c 100644 --- a/sysdeps/aarch64/fpu/bits/math-vector.h +++ b/sysdeps/aarch64/fpu/bits/math-vector.h @@ -33,6 +33,10 @@ # define __DECL_SIMD_acos __DECL_SIMD_aarch64 # undef __DECL_SIMD_acosf # define __DECL_SIMD_acosf __DECL_SIMD_aarch64 +# undef __DECL_SIMD_acosh +# define __DECL_SIMD_acosh __DECL_SIMD_aarch64 +# undef __DECL_SIMD_acoshf +# define __DECL_SIMD_acoshf __DECL_SIMD_aarch64 # undef __DECL_SIMD_asin # define __DECL_SIMD_asin __DECL_SIMD_aarch64 # undef __DECL_SIMD_asinf @@ -125,6 +129,7 @@ typedef __SVBool_t __sv_bool_t; __vpcs __f32x4_t _ZGVnN4vv_atan2f (__f32x4_t, __f32x4_t); __vpcs __f32x4_t _ZGVnN4v_acosf (__f32x4_t); +__vpcs __f32x4_t _ZGVnN4v_acoshf (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_asinf (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_atanf (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_cosf (__f32x4_t); @@ -143,6 +148,7 @@ __vpcs __f32x4_t _ZGVnN4v_tanf (__f32x4_t); __vpcs __f64x2_t _ZGVnN2vv_atan2 (__f64x2_t, __f64x2_t); __vpcs __f64x2_t _ZGVnN2v_acos (__f64x2_t); +__vpcs __f64x2_t _ZGVnN2v_acosh (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_asin (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_atan (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_cos (__f64x2_t); @@ -166,6 +172,7 @@ __vpcs __f64x2_t _ZGVnN2v_tan (__f64x2_t); __sv_f32_t _ZGVsMxvv_atan2f (__sv_f32_t, __sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_acosf (__sv_f32_t, __sv_bool_t); +__sv_f32_t _ZGVsMxv_acoshf (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_asinf (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_atanf (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_cosf (__sv_f32_t, __sv_bool_t); @@ -184,6 +191,7 @@ __sv_f32_t _ZGVsMxv_tanf (__sv_f32_t, __sv_bool_t); __sv_f64_t _ZGVsMxvv_atan2 (__sv_f64_t, __sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_acos (__sv_f64_t, __sv_bool_t); +__sv_f64_t _ZGVsMxv_acosh (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_asin (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_atan (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_cos (__sv_f64_t, __sv_bool_t); diff --git a/sysdeps/aarch64/fpu/sv_log1p_inline.h b/sysdeps/aarch64/fpu/sv_log1p_inline.h new file mode 100644 index 0000000000..da019674f9 --- /dev/null +++ b/sysdeps/aarch64/fpu/sv_log1p_inline.h @@ -0,0 +1,109 @@ +/* Helper for double-precision SVE routines which depend on log1p + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#ifndef AARCH64_FPU_SV_LOG1P_INLINE_H +#define AARCH64_FPU_SV_LOG1P_INLINE_H + +#include "sv_math.h" +#include "poly_sve_f64.h" + +static const struct sv_log1p_data +{ + double poly[19], ln2[2]; + uint64_t hf_rt2_top; + uint64_t one_m_hf_rt2_top; + uint32_t bottom_mask; + int64_t one_top; +} sv_log1p_data = { + /* Coefficients generated using Remez, deg=20, in [sqrt(2)/2-1, sqrt(2)-1]. + */ + .poly = { -0x1.ffffffffffffbp-2, 0x1.55555555551a9p-2, -0x1.00000000008e3p-2, + 0x1.9999999a32797p-3, -0x1.555555552fecfp-3, 0x1.249248e071e5ap-3, + -0x1.ffffff8bf8482p-4, 0x1.c71c8f07da57ap-4, -0x1.9999ca4ccb617p-4, + 0x1.7459ad2e1dfa3p-4, -0x1.554d2680a3ff2p-4, 0x1.3b4c54d487455p-4, + -0x1.2548a9ffe80e6p-4, 0x1.0f389a24b2e07p-4, -0x1.eee4db15db335p-5, + 0x1.e95b494d4a5ddp-5, -0x1.15fdf07cb7c73p-4, 0x1.0310b70800fcfp-4, + -0x1.cfa7385bdb37ep-6 }, + .ln2 = { 0x1.62e42fefa3800p-1, 0x1.ef35793c76730p-45 }, + .hf_rt2_top = 0x3fe6a09e00000000, + .one_m_hf_rt2_top = 0x00095f6200000000, + .bottom_mask = 0xffffffff, + .one_top = 0x3ff +}; + +static inline svfloat64_t +sv_log1p_inline (svfloat64_t x, const svbool_t pg) +{ + /* Helper for calculating log(x + 1). Adapted from v_log1p_inline.h, which + differs from v_log1p_2u5.c by: + - No special-case handling - this should be dealt with by the caller. + - Pairwise Horner polynomial evaluation for improved accuracy. + - Optionally simulate the shortcut for k=0, used in the scalar routine, + using svsel, for improved accuracy when the argument to log1p is close + to 0. This feature is enabled by defining WANT_SV_LOG1P_K0_SHORTCUT as 1 + in the source of the caller before including this file. + See sv_log1p_2u1.c for details of the algorithm. */ + const struct sv_log1p_data *d = ptr_barrier (&sv_log1p_data); + svfloat64_t m = svadd_x (pg, x, 1); + svuint64_t mi = svreinterpret_u64 (m); + svuint64_t u = svadd_x (pg, mi, d->one_m_hf_rt2_top); + + svint64_t ki + = svsub_x (pg, svreinterpret_s64 (svlsr_x (pg, u, 52)), d->one_top); + svfloat64_t k = svcvt_f64_x (pg, ki); + + /* Reduce x to f in [sqrt(2)/2, sqrt(2)]. */ + svuint64_t utop + = svadd_x (pg, svand_x (pg, u, 0x000fffff00000000), d->hf_rt2_top); + svuint64_t u_red = svorr_x (pg, utop, svand_x (pg, mi, d->bottom_mask)); + svfloat64_t f = svsub_x (pg, svreinterpret_f64 (u_red), 1); + + /* Correction term c/m. */ + svfloat64_t c = svsub_x (pg, x, svsub_x (pg, m, 1)); + svfloat64_t cm; + +#ifndef WANT_SV_LOG1P_K0_SHORTCUT +#error \ + "Cannot use sv_log1p_inline.h without specifying whether you need the k0 shortcut for greater accuracy close to 0" +#elif WANT_SV_LOG1P_K0_SHORTCUT + /* Shortcut if k is 0 - set correction term to 0 and f to x. The result is + that the approximation is solely the polynomial. */ + svbool_t knot0 = svcmpne (pg, k, 0); + cm = svdiv_z (knot0, c, m); + if (__glibc_likely (!svptest_any (pg, knot0))) + { + f = svsel (knot0, f, x); + } +#else + /* No shortcut. */ + cm = svdiv_x (pg, c, m); +#endif + + /* Approximate log1p(f) on the reduced input using a polynomial. */ + svfloat64_t f2 = svmul_x (pg, f, f); + svfloat64_t p = sv_pw_horner_18_f64_x (pg, f, f2, d->poly); + + /* Assemble log1p(x) = k * log2 + log1p(f) + c/m. */ + svfloat64_t ylo = svmla_x (pg, cm, k, d->ln2[0]); + svfloat64_t yhi = svmla_x (pg, f, k, d->ln2[1]); + + return svmla_x (pg, svadd_x (pg, ylo, yhi), f2, p); +} + +#endif diff --git a/sysdeps/aarch64/fpu/sv_log1pf_inline.h b/sysdeps/aarch64/fpu/sv_log1pf_inline.h new file mode 100644 index 0000000000..b94b2da055 --- /dev/null +++ b/sysdeps/aarch64/fpu/sv_log1pf_inline.h @@ -0,0 +1,76 @@ +/* Helper for single-precision SVE routines which depend on log1p + + Copyright (C) 2024 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#ifndef AARCH64_FPU_SV_LOG1PF_INLINE_H +#define AARCH64_FPU_SV_LOG1PF_INLINE_H + +#include "sv_math.h" +#include "vecmath_config.h" +#include "poly_sve_f32.h" + +static const struct sv_log1pf_data +{ + float32_t poly[9]; + f |
