aarch64/fpu: Add vector variants of cosh

Reviewed-by: Szabolcs Nagy <szabolcs.nagy@arm.com>

18 files changed, 648 insertions, 1 deletions

diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile
index 320b6ed43a..019c3a5188 100644
--- a/sysdeps/aarch64/fpu/Makefile
+++ b/sysdeps/aarch64/fpu/Makefile
@@ -3,6 +3,7 @@ libmvec-supported-funcs = acos \
                           atan \
                           atan2 \
                           cos \
+                          cosh \
                           erf \
                           exp \
                           exp10 \
@@ -32,7 +33,8 @@ libmvec-support = $(addsuffix f_advsimd,$(float-advsimd-funcs)) \
                   erf_data \
                   erff_data \
                   sv_erf_data \
-                  sv_erff_data
+                  sv_erff_data \
+                  v_exp_tail_data
 endif
 
 sve-cflags = -march=armv8-a+sve
diff --git a/sysdeps/aarch64/fpu/Versions b/sysdeps/aarch64/fpu/Versions
index d7b1e87191..884b4b57f0 100644
--- a/sysdeps/aarch64/fpu/Versions
+++ b/sysdeps/aarch64/fpu/Versions
@@ -79,6 +79,11 @@ libmvec {
     _ZGVsMxv_tan;
   }
   GLIBC_2.40 {
+    _ZGVnN2v_cosh;
+    _ZGVnN2v_coshf;
+    _ZGVnN4v_coshf;
+    _ZGVsMxv_cosh;
+    _ZGVsMxv_coshf;
     _ZGVnN2v_erf;
     _ZGVnN2v_erff;
     _ZGVnN4v_erff;
diff --git a/sysdeps/aarch64/fpu/advsimd_f32_protos.h b/sysdeps/aarch64/fpu/advsimd_f32_protos.h
index d8d88de218..c63b2948d4 100644
--- a/sysdeps/aarch64/fpu/advsimd_f32_protos.h
+++ b/sysdeps/aarch64/fpu/advsimd_f32_protos.h
@@ -21,6 +21,7 @@ libmvec_hidden_proto (V_NAME_F1(acos));
 libmvec_hidden_proto (V_NAME_F1(asin));
 libmvec_hidden_proto (V_NAME_F1(atan));
 libmvec_hidden_proto (V_NAME_F1(cos));
+libmvec_hidden_proto (V_NAME_F1(cosh));
 libmvec_hidden_proto (V_NAME_F1(erf));
 libmvec_hidden_proto (V_NAME_F1(exp10));
 libmvec_hidden_proto (V_NAME_F1(exp2));
diff --git a/sysdeps/aarch64/fpu/bits/math-vector.h b/sysdeps/aarch64/fpu/bits/math-vector.h
index 71f53363a0..8ca5509870 100644
--- a/sysdeps/aarch64/fpu/bits/math-vector.h
+++ b/sysdeps/aarch64/fpu/bits/math-vector.h
@@ -49,6 +49,10 @@
 # define __DECL_SIMD_cos __DECL_SIMD_aarch64
 # undef __DECL_SIMD_cosf
 # define __DECL_SIMD_cosf __DECL_SIMD_aarch64
+# undef __DECL_SIMD_cosh
+# define __DECL_SIMD_cosh __DECL_SIMD_aarch64
+# undef __DECL_SIMD_coshf
+# define __DECL_SIMD_coshf __DECL_SIMD_aarch64
 # undef __DECL_SIMD_erf
 # define __DECL_SIMD_erf __DECL_SIMD_aarch64
 # undef __DECL_SIMD_erff
@@ -124,6 +128,7 @@ __vpcs __f32x4_t _ZGVnN4v_acosf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_asinf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_atanf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_cosf (__f32x4_t);
+__vpcs __f32x4_t _ZGVnN4v_coshf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_erff (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_expf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_exp10f (__f32x4_t);
@@ -141,6 +146,7 @@ __vpcs __f64x2_t _ZGVnN2v_acos (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_asin (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_atan (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_cos (__f64x2_t);
+__vpcs __f64x2_t _ZGVnN2v_cosh (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_erf (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_exp (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_exp10 (__f64x2_t);
@@ -163,6 +169,7 @@ __sv_f32_t _ZGVsMxv_acosf (__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);
+__sv_f32_t _ZGVsMxv_coshf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_erff (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_expf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_exp10f (__sv_f32_t, __sv_bool_t);
@@ -180,6 +187,7 @@ __sv_f64_t _ZGVsMxv_acos (__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);
+__sv_f64_t _ZGVsMxv_cosh (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_erf (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_exp (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_exp10 (__sv_f64_t, __sv_bool_t);
diff --git a/sysdeps/aarch64/fpu/cosh_advsimd.c b/sysdeps/aarch64/fpu/cosh_advsimd.c
new file mode 100644
index 0000000000..ec7b59637e
--- /dev/null
+++ b/sysdeps/aarch64/fpu/cosh_advsimd.c
@@ -0,0 +1,108 @@
+/* Double-precision vector (AdvSIMD) cosh 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_math.h"
+
+static const struct data
+{
+  float64x2_t poly[3];
+  float64x2_t inv_ln2, ln2, shift, thres;
+  uint64x2_t index_mask, special_bound;
+} data = {
+  .poly = { V2 (0x1.fffffffffffd4p-2), V2 (0x1.5555571d6b68cp-3),
+	    V2 (0x1.5555576a59599p-5), },
+
+  .inv_ln2 = V2 (0x1.71547652b82fep8), /* N/ln2.  */
+  /* -ln2/N.  */
+  .ln2 = {-0x1.62e42fefa39efp-9, -0x1.abc9e3b39803f3p-64},
+  .shift = V2 (0x1.8p+52),
+  .thres = V2 (704.0),
+
+  .index_mask = V2 (0xff),
+  /* 0x1.6p9, above which exp overflows.  */
+  .special_bound = V2 (0x4086000000000000),
+};
+
+static float64x2_t NOINLINE VPCS_ATTR
+special_case (float64x2_t x, float64x2_t y, uint64x2_t special)
+{
+  return v_call_f64 (cosh, x, y, special);
+}
+
+/* Helper for approximating exp(x). Copied from v_exp_tail, with no
+   special-case handling or tail.  */
+static inline float64x2_t
+exp_inline (float64x2_t x)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  /* n = round(x/(ln2/N)).  */
+  float64x2_t z = vfmaq_f64 (d->shift, x, d->inv_ln2);
+  uint64x2_t u = vreinterpretq_u64_f64 (z);
+  float64x2_t n = vsubq_f64 (z, d->shift);
+
+  /* r = x - n*ln2/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);
+
+  /* y = tail + exp(r) - 1 ~= r + C1 r^2 + C2 r^3 + C3 r^4.  */
+  float64x2_t y = vfmaq_f64 (d->poly[1], d->poly[2], r);
+  y = vfmaq_f64 (d->poly[0], y, r);
+  y = vmulq_f64 (vfmaq_f64 (v_f64 (1), y, r), r);
+
+  /* s = 2^(n/N).  */
+  u = v_lookup_u64 (__v_exp_tail_data, i);
+  float64x2_t s = vreinterpretq_f64_u64 (vaddq_u64 (u, e));
+
+  return vfmaq_f64 (s, y, s);
+}
+
+/* Approximation for vector double-precision cosh(x) using exp_inline.
+   cosh(x) = (exp(x) + exp(-x)) / 2.
+   The greatest observed error is in the scalar fall-back region, so is the
+   same as the scalar routine, 1.93 ULP:
+   _ZGVnN2v_cosh (0x1.628af341989dap+9) got 0x1.fdf28623ef921p+1021
+				       want 0x1.fdf28623ef923p+1021.
+
+   The greatest observed error in the non-special region is 1.54 ULP:
+   _ZGVnN2v_cosh (0x1.8e205b6ecacf7p+2) got 0x1.f711dcb0c77afp+7
+				       want 0x1.f711dcb0c77b1p+7.  */
+float64x2_t VPCS_ATTR V_NAME_D1 (cosh) (float64x2_t x)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  float64x2_t ax = vabsq_f64 (x);
+  uint64x2_t special
+      = vcgtq_u64 (vreinterpretq_u64_f64 (ax), d->special_bound);
+
+  /* Up to the point that exp overflows, we can use it to calculate cosh by
+     exp(|x|) / 2 + 1 / (2 * exp(|x|)).  */
+  float64x2_t t = exp_inline (ax);
+  float64x2_t half_t = vmulq_n_f64 (t, 0.5);
+  float64x2_t half_over_t = vdivq_f64 (v_f64 (0.5), t);
+
+  /* Fall back to scalar for any special cases.  */
+  if (__glibc_unlikely (v_any_u64 (special)))
+    return special_case (x, vaddq_f64 (half_t, half_over_t), special);
+
+  return vaddq_f64 (half_t, half_over_t);
+}
diff --git a/sysdeps/aarch64/fpu/cosh_sve.c b/sysdeps/aarch64/fpu/cosh_sve.c
new file mode 100644
index 0000000000..919f34604a
--- /dev/null
+++ b/sysdeps/aarch64/fpu/cosh_sve.c
@@ -0,0 +1,105 @@
+/* Double-precision vector (SVE) cosh 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 "sv_math.h"
+
+static const struct data
+{
+  float64_t poly[3];
+  float64_t inv_ln2, ln2_hi, ln2_lo, shift, thres;
+  uint64_t index_mask, special_bound;
+} data = {
+  .poly = { 0x1.fffffffffffd4p-2, 0x1.5555571d6b68cp-3,
+	    0x1.5555576a59599p-5, },
+
+  .inv_ln2 = 0x1.71547652b82fep8, /* N/ln2.  */
+  /* -ln2/N.  */
+  .ln2_hi = -0x1.62e42fefa39efp-9,
+  .ln2_lo = -0x1.abc9e3b39803f3p-64,
+  .shift = 0x1.8p+52,
+  .thres = 704.0,
+
+  .index_mask = 0xff,
+  /* 0x1.6p9, above which exp overflows.  */
+  .special_bound = 0x4086000000000000,
+};
+
+static svfloat64_t NOINLINE
+special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
+{
+  return sv_call_f64 (cosh, x, y, special);
+}
+
+/* Helper for approximating exp(x). Copied from sv_exp_tail, with no
+   special-case handling or tail.  */
+static inline svfloat64_t
+exp_inline (svfloat64_t x, const svbool_t pg, const struct data *d)
+{
+  /* Calculate exp(x).  */
+  svfloat64_t z = svmla_x (pg, sv_f64 (d->shift), x, d->inv_ln2);
+  svfloat64_t n = svsub_x (pg, z, d->shift);
+
+  svfloat64_t r = svmla_x (pg, x, n, d->ln2_hi);
+  r = svmla_x (pg, r, n, d->ln2_lo);
+
+  svuint64_t u = svreinterpret_u64 (z);
+  svuint64_t e = svlsl_x (pg, u, 52 - V_EXP_TAIL_TABLE_BITS);
+  svuint64_t i = svand_x (pg, u, d->index_mask);
+
+  svfloat64_t y = svmla_x (pg, sv_f64 (d->poly[1]), r, d->poly[2]);
+  y = svmla_x (pg, sv_f64 (d->poly[0]), r, y);
+  y = svmla_x (pg, sv_f64 (1.0), r, y);
+  y = svmul_x (pg, r, y);
+
+  /* s = 2^(n/N).  */
+  u = svld1_gather_index (pg, __v_exp_tail_data, i);
+  svfloat64_t s = svreinterpret_f64 (svadd_x (pg, u, e));
+
+  return svmla_x (pg, s, s, y);
+}
+
+/* Approximation for SVE double-precision cosh(x) using exp_inline.
+   cosh(x) = (exp(x) + exp(-x)) / 2.
+   The greatest observed error is in the scalar fall-back region, so is the
+   same as the scalar routine, 1.93 ULP:
+   _ZGVsMxv_cosh (0x1.628ad45039d2fp+9) got 0x1.fd774e958236dp+1021
+				       want 0x1.fd774e958236fp+1021.
+
+   The greatest observed error in the non-special region is 1.54 ULP:
+   _ZGVsMxv_cosh (0x1.ba5651dd4486bp+2) got 0x1.f5e2bb8d5c98fp+8
+				       want 0x1.f5e2bb8d5c991p+8.  */
+svfloat64_t SV_NAME_D1 (cosh) (svfloat64_t x, const svbool_t pg)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  svfloat64_t ax = svabs_x (pg, x);
+  svbool_t special = svcmpgt (pg, svreinterpret_u64 (ax), d->special_bound);
+
+  /* Up to the point that exp overflows, we can use it to calculate cosh by
+     exp(|x|) / 2 + 1 / (2 * exp(|x|)).  */
+  svfloat64_t t = exp_inline (ax, pg, d);
+  svfloat64_t half_t = svmul_x (pg, t, 0.5);
+  svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
+
+  /* Fall back to scalar for any special cases.  */
+  if (__glibc_unlikely (svptest_any (pg, special)))
+    return special_case (x, svadd_x (pg, half_t, half_over_t), special);
+
+  return svadd_x (pg, half_t, half_over_t);
+}
diff --git a/sysdeps/aarch64/fpu/coshf_advsimd.c b/sysdeps/aarch64/fpu/coshf_advsimd.c
new file mode 100644
index 0000000000..c1ab4923b8
--- /dev/null
+++ b/sysdeps/aarch64/fpu/coshf_advsimd.c
@@ -0,0 +1,84 @@
+/* Single-precision vector (AdvSIMD) cosh 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_expf_inline.h"
+#include "v_math.h"
+
+static const struct data
+{
+  struct v_expf_data expf_consts;
+  uint32x4_t tiny_bound, special_bound;
+} data = {
+  .expf_consts = V_EXPF_DATA,
+  .tiny_bound = V4 (0x20000000), /* 0x1p-63: Round to 1 below this.  */
+  /* 0x1.5a92d8p+6: expf overflows above this, so have to use special case.  */
+  .special_bound = V4 (0x42ad496c),
+};
+
+#if !WANT_SIMD_EXCEPT
+static float32x4_t NOINLINE VPCS_ATTR
+special_case (float32x4_t x, float32x4_t y, uint32x4_t special)
+{
+  return v_call_f32 (coshf, x, y, special);
+}
+#endif
+
+/* Single-precision vector cosh, using vector expf.
+   Maximum error is 2.38 ULP:
+   _ZGVnN4v_coshf (0x1.e8001ep+1) got 0x1.6a491ep+4
+				 want 0x1.6a4922p+4.  */
+float32x4_t VPCS_ATTR V_NAME_F1 (cosh) (float32x4_t x)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  float32x4_t ax = vabsq_f32 (x);
+  uint32x4_t iax = vreinterpretq_u32_f32 (ax);
+  uint32x4_t special = vcgeq_u32 (iax, d->special_bound);
+
+#if WANT_SIMD_EXCEPT
+  /* If fp exceptions are to be triggered correctly, fall back to the scalar
+     variant for all inputs if any input is a special value or above the bound
+     at which expf overflows.  */
+  if (__glibc_unlikely (v_any_u32 (special)))
+    return v_call_f32 (coshf, x, x, v_u32 (-1));
+
+  uint32x4_t tiny = vcleq_u32 (iax, d->tiny_bound);
+  /* If any input is tiny, avoid underflow exception by fixing tiny lanes of
+     input to 0, which will generate no exceptions.  */
+  if (__glibc_unlikely (v_any_u32 (tiny)))
+    ax = v_zerofy_f32 (ax, tiny);
+#endif
+
+  /* Calculate cosh by exp(x) / 2 + exp(-x) / 2.  */
+  float32x4_t t = v_expf_inline (ax, &d->expf_consts);
+  float32x4_t half_t = vmulq_n_f32 (t, 0.5);
+  float32x4_t half_over_t = vdivq_f32 (v_f32 (0.5), t);
+
+#if WANT_SIMD_EXCEPT
+  if (__glibc_unlikely (v_any_u32 (tiny)))
+    return vbslq_f32 (tiny, v_f32 (1), vaddq_f32 (half_t, half_over_t));
+#else
+  if (__glibc_unlikely (v_any_u32 (special)))
+    return special_case (x, vaddq_f32 (half_t, half_over_t), special);
+#endif
+
+  return vaddq_f32 (half_t, half_over_t);
+}
+libmvec_hidden_def (V_NAME_F1 (cosh))
+HALF_WIDTH_ALIAS_F1 (cosh)
diff --git a/sysdeps/aarch64/fpu/coshf_sve.c b/sysdeps/aarch64/fpu/coshf_sve.c
new file mode 100644
index 0000000000..e5d8a299c6
--- /dev/null
+++ b/sysdeps/aarch64/fpu/coshf_sve.c
@@ -0,0 +1,59 @@
+/* Single-precision vector (SVE) cosh 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 "sv_math.h"
+#include "sv_expf_inline.h"
+
+static const struct data
+{
+  struct sv_expf_data expf_consts;
+  uint32_t special_bound;
+} data = {
+  .expf_consts = SV_EXPF_DATA,
+  /* 0x1.5a92d8p+6: expf overflows above this, so have to use special case.  */
+  .special_bound = 0x42ad496c,
+};
+
+static svfloat32_t NOINLINE
+special_case (svfloat32_t x, svfloat32_t y, svbool_t pg)
+{
+  return sv_call_f32 (coshf, x, y, pg);
+}
+
+/* Single-precision vector cosh, using vector expf.
+   Maximum error is 1.89 ULP:
+   _ZGVsMxv_coshf (-0x1.65898cp+6) got 0x1.f00aep+127
+				  want 0x1.f00adcp+127.  */
+svfloat32_t SV_NAME_F1 (cosh) (svfloat32_t x, svbool_t pg)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  svfloat32_t ax = svabs_x (pg, x);
+  svbool_t special = svcmpge (pg, svreinterpret_u32 (ax), d->special_bound);
+
+  /* Calculate cosh by exp(x) / 2 + exp(-x) / 2.  */
+  svfloat32_t t = expf_inline (ax, pg, &d->expf_consts);
+  svfloat32_t half_t = svmul_x (pg, t, 0.5);
+  svfloat32_t half_over_t = svdivr_x (pg, t, 0.5);
+
+  if (__glibc_unlikely (svptest_any (pg, special)))
+    return special_case (x, svadd_x (pg, half_t, half_over_t), special);
+
+  return svadd_x (pg, half_t, half_over_t);
+}
diff --git a/sysdeps/aarch64/fpu/sv_expf_inline.h b/sysdeps/aarch64/fpu/sv_expf_inline.h
new file mode 100644
index 0000000000..23963b5f8e
--- /dev/null
+++ b/sysdeps/aarch64/fpu/sv_expf_inline.h
@@ -0,0 +1,75 @@
+/* SVE helper for single-precision routines which depend on exp
+
+   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_EXPF_INLINE_H
+#define AARCH64_FPU_SV_EXPF_INLINE_H
+
+#include "sv_math.h"
+
+struct sv_expf_data
+{
+  float poly[5];
+  float inv_ln2, ln2_hi, ln2_lo, shift;
+};
+
+/* Coefficients copied from the polynomial in AdvSIMD variant, reversed for
+   compatibility with polynomial helpers. Shift is 1.5*2^17 + 127.  */
+#define SV_EXPF_DATA                                                          \
+  {                                                                           \
+    .poly = { 0x1.ffffecp-1f, 0x1.fffdb6p-2f, 0x1.555e66p-3f, 0x1.573e2ep-5f, \
+	      0x1.0e4020p-7f },                                               \
+                                                                              \
+    .inv_ln2 = 0x1.715476p+0f, .ln2_hi = 0x1.62e4p-1f,                        \
+    .ln2_lo = 0x1.7f7d1cp-20f, .shift = 0x1.803f8p17f,                        \
+