[PATCH] Inverse Hyperbolic
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Thu Feb 26 12:02:27 GMT 2026
What does this patch aim to accomplish? Is it a performance improvement, correctness
fix or some other improvement?
On 25/02/26 12:51, Cosmina.Dunca@arm.com wrote:
> ---
> sysdeps/aarch64/fpu/acosh_advsimd.c | 13 ++++-
> sysdeps/aarch64/fpu/acosh_sve.c | 32 ++++++++++--
> sysdeps/aarch64/fpu/acoshf_advsimd.c | 73 +++++++++++++++-------------
> sysdeps/aarch64/fpu/asinh_advsimd.c | 68 +++++++++++++++++---------
> sysdeps/aarch64/fpu/asinh_sve.c | 48 ++++++++++++------
> sysdeps/aarch64/fpu/asinhf_advsimd.c | 63 ++++++++++++++----------
> sysdeps/aarch64/fpu/atanh_advsimd.c | 27 +++++-----
> sysdeps/aarch64/fpu/atanh_sve.c | 24 ++++++---
> 8 files changed, 224 insertions(+), 124 deletions(-)
>
> diff --git a/sysdeps/aarch64/fpu/acosh_advsimd.c b/sysdeps/aarch64/fpu/acosh_advsimd.c
> index ac0db3ed15..b5efe9222f 100644
> --- a/sysdeps/aarch64/fpu/acosh_advsimd.c
> +++ b/sysdeps/aarch64/fpu/acosh_advsimd.c
> @@ -34,7 +34,18 @@ 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);
> + /* Acosh is define on [1, inf). Its formula can be re-written knowing that 1
> + becomes negligible when x is a very large number. So for special numbers,
> + where x >= 2^511, acosh ~= ln(2x). But, ln(2x) = ln(2) + ln(x) and below we
> + calculate ln(x) and then add ln(2) to the result.
> +
> + Right before returning we check if x is infinity or if x is lower than 1,
> + in which case we return infinity or NaN. */
> + float64x2_t one = v_f64 (1.0);
> + float64x2_t t = log1p_inline (vbslq_f64 (special, x, y), d);
> + t = vbslq_f64 (special, vaddq_f64 (t, v_f64 (d->ln2[0])), t);
> + t = vbslq_f64 (vceqq_f64 (x, v_f64 (INFINITY)), v_f64 (INFINITY), t);
> + return vbslq_f64 (vcltq_f64 (x, one), v_f64 (NAN), t);
> }
>
> /* Vector approximation for double-precision acosh, based on log1p.
> diff --git a/sysdeps/aarch64/fpu/acosh_sve.c b/sysdeps/aarch64/fpu/acosh_sve.c
> index 1f6b8fc19a..6d996e3d36 100644
> --- a/sysdeps/aarch64/fpu/acosh_sve.c
> +++ b/sysdeps/aarch64/fpu/acosh_sve.c
> @@ -20,13 +20,33 @@
> #define WANT_SV_LOG1P_K0_SHORTCUT 1
> #include "sv_log1p_inline.h"
>
> -#define One (0x3ff0000000000000)
> -#define Thres (0x1ff0000000000000) /* asuint64 (0x1p511) - One. */
> +#define One (0x3ff0000000000000ULL)
> +#define Thres (0x1ff0000000000000ULL)
>
> +const static struct data
> +{
> + double ln2;
> +} data = { .ln2 = 0x1.62e42fefa39efp-1 };
> +
> +/* Acosh is define on [1, inf). Its formula can be re-written knowing that 1
> + becomes negligible when x is a very large number. So for special numbers,
> + where x >= 2^511, acosh ~= ln(2x). But, ln(2x) = ln(2) + ln(x) and below we
> + calculate ln(x) and then add ln(2) to the result.
> +
> + Right before returning we check if x is infinity or if x is lower than 1,
> + in which case we return infinity or NaN. */
> static svfloat64_t NOINLINE
> -special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
> +special_case (svfloat64_t x, svfloat64_t xm1, svfloat64_t y, svbool_t special,
> + svbool_t pg, const struct data *d)
> {
> - return sv_call_f64 (acosh, x, y, special);
> + svfloat64_t inf = sv_f64 (INFINITY);
> + svfloat64_t nan = sv_f64 (NAN);
> + svfloat64_t log = sv_log1p_inline (svsel (special, xm1, y), pg);
> + svfloat64_t result = svadd_m (special, log, sv_f64 (d->ln2));
> + svbool_t is_inf_nan = svorr_b_z (special, svcmpge (special, x, inf),
> + svcmplt (special, x, sv_f64 (1.0)));
> + svfloat64_t inf_or_nan = svsel (svcmpge (special, x, inf), inf, nan);
> + return svsel (is_inf_nan, inf_or_nan, result);
> }
>
> /* SVE approximation for double-precision acosh, based on log1p.
> @@ -36,6 +56,8 @@ special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
> want 0x1.ef0cee7c33ce4p-2. */
> svfloat64_t SV_NAME_D1 (acosh) (svfloat64_t x, const svbool_t pg)
> {
> + const struct data *d = ptr_barrier (&data);
> +
> /* (ix - One) >= (BigBound - One). */
> svuint64_t ix = svreinterpret_u64 (x);
> svbool_t special = svcmpge (pg, svsub_x (pg, ix, One), Thres);
> @@ -46,6 +68,6 @@ svfloat64_t SV_NAME_D1 (acosh) (svfloat64_t x, const svbool_t pg)
>
> /* 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 special_case (x, xm1, y, special, pg, d);
> return sv_log1p_inline (y, pg);
> }
> diff --git a/sysdeps/aarch64/fpu/acoshf_advsimd.c b/sysdeps/aarch64/fpu/acoshf_advsimd.c
> index 6071190938..275de77e1b 100644
> --- a/sysdeps/aarch64/fpu/acoshf_advsimd.c
> +++ b/sysdeps/aarch64/fpu/acoshf_advsimd.c
> @@ -24,15 +24,13 @@ const static struct data
> struct v_log1pf_data log1pf_consts;
> uint32x4_t one;
> uint16x4_t special_bound_u16;
> - uint32x4_t special_bound_u32;
> - float32x4_t pinf, nan;
> + float32x4_t inf, nan;
> } data = {
> .log1pf_consts = V_LOG1PF_CONSTANTS_TABLE,
> .one = V4 (0x3f800000),
> - .special_bound_u16 = V4 (0x2000),
> /* asuint(sqrt(FLT_MAX)) - asuint(1). */
> - .special_bound_u32 = V4 (0x20000000),
> - .pinf = V4 (INFINITY),
> + .special_bound_u16 = V4 (0x2000),
> + .inf = V4 (INFINITY),
> .nan = V4 (NAN),
> };
>
> @@ -51,34 +49,42 @@ inline_acoshf (float32x4_t x, const struct data *d)
> }
>
> static float32x4_t VPCS_ATTR NOINLINE
> -special_case (float32x4_t x, const struct data *d)
> +special_case (float32x4_t x, uint16x4_t special_u16, const struct data *d)
> {
> - uint32x4_t special = vcgeq_u32 (
> - vsubq_u32 (vreinterpretq_u32_f32 (x), d->one), d->special_bound_u32);
> -
> - /* To avoid the overflow in x^2 (so the x < sqrt(FLT_MAX) constraint), we
> - reduce the input of acosh to a narrower interval by relying on the identity
> - acosh(t) = 1/2acosh(2t^2 - 1) for t>=1.
> - If we set t=sqrt((x+1)/2), since x>=1 then t>=sqrt(2/2)=1, and therefore
> - acosh(x) = 2acosh(sqrt((x+1)/2)). */
> - float32x4_t r = vaddq_f32 (x, vreinterpretq_f32_u32 (d->one));
> - r = vmulq_f32 (r, v_f32 (0.5f));
> - r = vbslq_f32 (special, vsqrtq_f32 (r), x);
> -
> - float32x4_t y = inline_acoshf (r, d);
> -
> - y = vbslq_f32 (special, vmulq_f32 (y, v_f32 (2.0f)), y);
> -
> - /* Check whether x is less than 1, or x is inf or nan. */
> - uint32x4_t inf_minus_one
> - = vsubq_u32 (vreinterpretq_u32_f32 (d->pinf), d->one);
> - uint32x4_t is_infnan = vcgeq_u32 (
> - vsubq_u32 (vreinterpretq_u32_f32 (x), d->one), inf_minus_one);
> -
> - y = vbslq_f32 (is_infnan, d->nan, y);
> - uint32x4_t ret_pinf = vceqq_f32 (x, d->pinf);
> - y = vbslq_f32 (ret_pinf, d->pinf, y);
> - return y;
> + float32x4_t one = v_f32 (1.0);
> + float32x4_t ln2 = d->log1pf_consts.ln2;
> +
> + /* acosh(x) = ln(x + sqrt[x^2 -1]).
> + So acosh(x) = log1p (x + sqrt[x^2 - 1] - 1). */
> + float32x4_t xm1 = vsubq_f32 (x, one);
> + float32x4_t u = vmulq_f32 (xm1, vaddq_f32 (x, one));
> +
> + float32x4_t y = vaddq_f32 (xm1, vsqrtq_f32 (u));
> +
> + /* special_u16 will be 0x0000ffff for true lanes, which doesn't work with bit
> + selects therefore to make it a 32 bit predicate, we have to add special
> + << 16. */
> + uint32x4_t special = vmovl_u16 (special_u16);
> + special = vaddq_u32 (special, vshlq_n_u32 (special, 16));
> + float32x4_t xy = vbslq_f32 (special, x, y);
> +
> + /* For large inputs, acosh(x) ≈ log(x) + ln(2).
> + We use log1pf-inline log implementation and add ln(2). */
> + float32x4_t log_xy = log1pf_inline (xy, &d->log1pf_consts);
> +
> + /* For acoshf there are three special cases that need considering. Infinity
> + and NaNs, which are also returned unchanged and for cases of x < 1 we'll
> + return all NaNs since acoshf is defined on (1, inf). */
> + uint32x4_t is_lower_one = vcltq_f32 (x, one);
> + uint32x4_t is_finite = vcltq_f32 (x, d->inf);
> +
> + /* This dependency chain of selects can run in parallel with y and log_x
> + being calculated before the last addition. */
> + float32x4_t ln2_inf_nan = vbslq_f32 (is_finite, ln2, x);
> + ln2_inf_nan = vbslq_f32 (is_lower_one, d->nan, ln2_inf_nan);
> + float32x4_t ln2_inf_nan_zero = vbslq_f32 (special, ln2_inf_nan, v_f32 (0));
> +
> + return vaddq_f32 (log_xy, ln2_inf_nan_zero);
> }
>
> /* Vector approximation for single-precision acosh, based on log1p.
> @@ -97,8 +103,9 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F1 (acosh) (float32x4_t x)
> = vcge_u16 (vsubhn_u32 (ix, d->one), d->special_bound_u16);
>
> if (__glibc_unlikely (v_any_u16h (special)))
> - return special_case (x, d);
> + return special_case (x, special, d);
> return inline_acoshf (x, d);
> }
> +
> libmvec_hidden_def (V_NAME_F1 (acosh))
> HALF_WIDTH_ALIAS_F1 (acosh)
> diff --git a/sysdeps/aarch64/fpu/asinh_advsimd.c b/sysdeps/aarch64/fpu/asinh_advsimd.c
> index 51689179ff..650ed46eec 100644
> --- a/sysdeps/aarch64/fpu/asinh_advsimd.c
> +++ b/sysdeps/aarch64/fpu/asinh_advsimd.c
> @@ -27,36 +27,43 @@ const static struct data
> double lc1, lc3, ln2, lc4;
> float64x2_t c0, c2, c4, c6, c8, c10, c12, c14, c16, c17;
> double c1, c3, c5, c7, c9, c11, c13, c15;
> + float64x2_t inf;
> } data = {
>
> /* Even terms of polynomial s.t. asinh(x) is approximated by
> asinh(x) ~= x + x^3 * (C0 + C1 * x + C2 * x^2 + C3 * x^3 + ...).
> Generated using Remez, f = (asinh(sqrt(x)) - sqrt(x))/x^(3/2). */
> - .c0 = V2 (-0x1.55555555554a7p-3), .c1 = 0x1.3333333326c7p-4,
> - .c2 = V2 (-0x1.6db6db68332e6p-5), .c3 = 0x1.f1c71b26fb40dp-6,
> - .c4 = V2 (-0x1.6e8b8b654a621p-6), .c5 = 0x1.1c4daa9e67871p-6,
> - .c6 = V2 (-0x1.c9871d10885afp-7), .c7 = 0x1.7a16e8d9d2ecfp-7,
> - .c8 = V2 (-0x1.3ddca533e9f54p-7), .c9 = 0x1.0becef748dafcp-7,
> - .c10 = V2 (-0x1.b90c7099dd397p-8), .c11 = 0x1.541f2bb1ffe51p-8,
> - .c12 = V2 (-0x1.d217026a669ecp-9), .c13 = 0x1.0b5c7977aaf7p-9,
> - .c14 = V2 (-0x1.e0f37daef9127p-11), .c15 = 0x1.388b5fe542a6p-12,
> - .c16 = V2 (-0x1.021a48685e287p-14), .c17 = V2 (0x1.93d4ba83d34dap-18),
> - .lc0 = V2 (-0x1.ffffffffffff7p-2), .lc1 = 0x1.55555555170d4p-2,
> - .lc2 = V2 (-0x1.0000000399c27p-2), .lc3 = 0x1.999b2e90e94cap-3,
> - .lc4 = -0x1.554e550bd501ep-3, .ln2 = 0x1.62e42fefa39efp-1,
> - .off = V2 (0x3fe6900900000000), .huge_bound = V2 (0x5fe0000000000000),
> - .abs_mask = V2 (0x7fffffffffffffff), .mask = V2 (0xfffULL << 52),
> + .c0 = V2 (-0x1.55555555554a7p-3),
> + .c1 = 0x1.3333333326c7p-4,
> + .c2 = V2 (-0x1.6db6db68332e6p-5),
> + .c3 = 0x1.f1c71b26fb40dp-6,
> + .c4 = V2 (-0x1.6e8b8b654a621p-6),
> + .c5 = 0x1.1c4daa9e67871p-6,
> + .c6 = V2 (-0x1.c9871d10885afp-7),
> + .c7 = 0x1.7a16e8d9d2ecfp-7,
> + .c8 = V2 (-0x1.3ddca533e9f54p-7),
> + .c9 = 0x1.0becef748dafcp-7,
> + .c10 = V2 (-0x1.b90c7099dd397p-8),
> + .c11 = 0x1.541f2bb1ffe51p-8,
> + .c12 = V2 (-0x1.d217026a669ecp-9),
> + .c13 = 0x1.0b5c7977aaf7p-9,
> + .c14 = V2 (-0x1.e0f37daef9127p-11),
> + .c15 = 0x1.388b5fe542a6p-12,
> + .c16 = V2 (-0x1.021a48685e287p-14),
> + .c17 = V2 (0x1.93d4ba83d34dap-18),
> + .lc0 = V2 (-0x1.ffffffffffff7p-2),
> + .lc1 = 0x1.55555555170d4p-2,
> + .lc2 = V2 (-0x1.0000000399c27p-2),
> + .lc3 = 0x1.999b2e90e94cap-3,
> + .lc4 = -0x1.554e550bd501ep-3,
> + .ln2 = 0x1.62e42fefa39efp-1,
> + .off = V2 (0x3fe6900900000000),
> + .huge_bound = V2 (0x5fe0000000000000),
> + .abs_mask = V2 (0x7fffffffffffffff),
> + .mask = V2 (0xfffULL << 52),
> + .inf = V2 (INFINITY)
> };
>
> -static float64x2_t NOINLINE VPCS_ATTR
> -special_case (float64x2_t x, float64x2_t y, uint64x2_t abs_mask,
> - uint64x2_t special)
> -{
> - /* Copy sign. */
> - y = vbslq_f64 (abs_mask, y, x);
> - return v_call_f64 (asinh, x, y, special);
> -}
> -
> #define N (1 << V_LOG_TABLE_BITS)
> #define IndexMask (N - 1)
>
> @@ -110,6 +117,19 @@ log_inline (float64x2_t ax, const struct data *d)
> return vfmaq_f64 (hi, y, r2);
> }
>
> +static float64x2_t NOINLINE VPCS_ATTR
> +special_case (float64x2_t x, uint64x2_t special, float64x2_t y, float64x2_t ax,
> + const struct data *d)
> +{
> + /* For very large inputs (x > 2^511), asinh(x) ≈ ln(2x).
> + In this range the +sqrt(x^2+1) term is negligible, so we compute
> + asinh(x) as ln(x) + ln(2). */
> + float64x2_t res_asinh = vaddq_f64 (log_inline (ax, d), v_f64 (d->ln2));
> + res_asinh = vbslq_f64 (vceqq_f64 (ax, d->inf), d->inf, res_asinh);
> + res_asinh = vbslq_f64 (special, res_asinh, y);
> + return vbslq_f64 (d->abs_mask, res_asinh, x);
> +}
> +
> /* Double-precision implementation of vector asinh(x).
> asinh is very sensitive around 1, so it is impractical to devise a single
> low-cost algorithm which is sufficiently accurate on a wide range of input.
> @@ -179,7 +199,7 @@ VPCS_ATTR float64x2_t V_NAME_D1 (asinh) (float64x2_t x)
> /* Choose the right option for each lane. */
> float64x2_t y = vbslq_f64 (gt1, option_1, option_2);
> if (__glibc_unlikely (v_any_u64 (special)))
> - return special_case (x, y, d->abs_mask, special);
> + return special_case (x, special, y, ax, d);
>
> /* Copy sign. */
> return vbslq_f64 (d->abs_mask, y, x);
> diff --git a/sysdeps/aarch64/fpu/asinh_sve.c b/sysdeps/aarch64/fpu/asinh_sve.c
> index bf9dad9e4b..cb32b60c3d 100644
> --- a/sysdeps/aarch64/fpu/asinh_sve.c
> +++ b/sysdeps/aarch64/fpu/asinh_sve.c
> @@ -29,7 +29,7 @@ static const struct data
> double even_coeffs[9];
> double ln2, p3, p1, p4, p0, p2, c1, c3, c5, c7, c9, c11, c13, c15, c17;
> uint64_t off, mask;
> -
> + double inf;
> } data = {
> /* Polynomial generated using Remez on [2^-26, 1]. */
> .even_coeffs ={
> @@ -61,14 +61,9 @@ static const struct data
> .p4 = -0x1.554e550bd501ep-3,
> .off = 0x3fe6900900000000,
> .mask = 0xfffULL << 52,
> + .inf = INFINITY
> };
>
> -static svfloat64_t NOINLINE
> -special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
> -{
> - return sv_call_f64 (asinh, x, y, special);
> -}
> -
> static inline svfloat64_t
> __sv_log_inline (svfloat64_t x, const struct data *d, const svbool_t pg)
> {
> @@ -104,12 +99,38 @@ __sv_log_inline (svfloat64_t x, const struct data *d, const svbool_t pg)
> return y;
> }
>
> +static svfloat64_t NOINLINE
> +special_case (svfloat64_t ax, svfloat64_t y, svuint64_t sign, svbool_t special,
> + svbool_t pg, const struct data *d)
> +{
> + /* For very large inputs (x > 2^511), asinh(x) ≈ ln(2x).
> + In this range the +sqrt(x^2+1) term is negligible, so we compute
> + asinh(x) as ln(x) + ln(2) later in this function. */
> + svfloat64_t log_ax = __sv_log_inline (ax, d, special);
> +
> + /* The only special cases that need considering are infinity and NaNs since
> + 0 will be handled by other calculations. */
> + svfloat64_t inf = sv_f64 (d->inf);
> + svbool_t is_inf = svcmpeq (special, ax, inf);
> + svbool_t is_nan = svcmpne (special, ax, ax);
> + svfloat64_t inf_ln2 = svsel (is_inf, inf, sv_f64 (d->ln2));
> + svfloat64_t inf_nan_ln2 = svsel (is_nan, sv_f64 (NAN), inf_ln2);
> + svfloat64_t asinh_x_res = svadd_x (special, log_ax, inf_nan_ln2);
> +
> + /* Now select (based on special) between x and y to change the type and,
> + return either the positive or negative value, considering the input and
> + its sign. */
> + svfloat64_t result = svsel (special, asinh_x_res, y);
> + svuint64_t result_uint = svreinterpret_u64 (result);
> + return svreinterpret_f64 (sveor_m (pg, result_uint, sign));
> +}
> +
> /* Double-precision implementation of SVE asinh(x).
> asinh is very sensitive around 1, so it is impractical to devise a single
> low-cost algorithm which is sufficiently accurate on a wide range of input.
> Instead we use two different algorithms:
> asinh(x) = sign(x) * log(|x| + sqrt(x^2 + 1) if |x| >= 1
> - = sign(x) * (|x| + |x|^3 * P(x^2)) otherwise
> + = sign(x) * (|x| + |x|^3 * P(x^2)) otherwise
> where log(x) is an optimized log approximation, and P(x) is a polynomial
> shared with the scalar routine. The greatest observed error 2.51 ULP, in
> |x| >= 1:
> @@ -180,14 +201,11 @@ svfloat64_t SV_NAME_D1 (asinh) (svfloat64_t x, const svbool_t pg)
> option_2 = svmla_x (pg, ax, p, svmul_x (svptrue_b64 (), x2, ax));
> }
>
> - if (__glibc_unlikely (svptest_any (pg, special)))
> - return special_case (
> - x,
> - svreinterpret_f64 (sveor_x (
> - pg, svreinterpret_u64 (svsel (ge1, option_1, option_2)), sign)),
> - special);
> -
> /* Choose the right option for each lane. */
> svfloat64_t y = svsel (ge1, option_1, option_2);
> +
> + if (__glibc_unlikely (svptest_any (pg, special)))
> + return special_case (ax, y, sign, special, pg, d);
> +
> return svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (y), sign));
> }
> diff --git a/sysdeps/aarch64/fpu/asinhf_advsimd.c b/sysdeps/aarch64/fpu/asinhf_advsimd.c
> index f1b1dd7a46..3454e465cb 100644
> --- a/sysdeps/aarch64/fpu/asinhf_advsimd.c
> +++ b/sysdeps/aarch64/fpu/asinhf_advsimd.c
> @@ -25,12 +25,12 @@ const static struct data
> struct v_log1pf_data log1pf_consts;
> float32x4_t one;
> uint32x4_t square_lim;
> - float32x4_t pinf, nan;
> + float32x4_t inf, nan;
> } data = {
> .one = V4 (1),
> .log1pf_consts = V_LOG1PF_CONSTANTS_TABLE,
> .square_lim = V4 (0x5f800000), /* asuint(sqrt(FLT_MAX)). */
> - .pinf = V4 (INFINITY),
> + .inf = V4 (INFINITY),
> .nan = V4 (NAN),
> };
>
> @@ -51,31 +51,39 @@ static float32x4_t VPCS_ATTR NOINLINE
> special_case (float32x4_t ax, uint32x4_t sign, uint32x4_t special,
> const struct data *d)
> {
> - /* To avoid overflow in x^2 (so the x < sqrt(FLT_MAX) constraint), we
> - reduce the input of asinh to a narrower interval by relying on the
> - identity: 2asinh(t) = +-acosh(2t^2 + 1)
> - If we set t=sqrt((x-1)/2), then
> - 2asinh(sqrt((x-1)/2)) = acosh(x).
> - Found that, for a high input x, asinh(x) very closely approximates
> - acosh(x), so implemented it with this function instead. */
> - float32x4_t r = vsubq_f32 (ax, d->one);
> - r = vmulq_f32 (r, v_f32 (0.5f));
> - r = vbslq_f32 (special, vsqrtq_f32 (r), ax);
> -
> - float32x4_t y = inline_asinhf (r, sign, d);
> -
> - y = vbslq_f32 (special, vmulq_f32 (y, v_f32 (2.0f)), y);
> -
> - /* Check whether x is inf or nan. */
> - uint32x4_t ret_inf = vceqq_f32 (ax, d->pinf);
> - uint32x4_t ret_nan = vmvnq_u32 (vcleq_f32 (ax, d->pinf));
> - y = vbslq_f32 (ret_inf, d->pinf, y);
> - y = vbslq_f32 (ret_nan, d->nan, y);
> - /* Put sign back in for minf, as it doesn't happen in log1pf_inline call. */
> - y = vbslq_f32 (
> - ret_inf,
> - vreinterpretq_f32_u32 (veorq_u32 (vreinterpretq_u32_f32 (y), sign)), y);
> - return y;
> + float32x4_t t
> + = vaddq_f32 (v_f32 (1.0f), vsqrtq_f32 (vfmaq_f32 (d->one, ax, ax)));
> + float32x4_t y = vaddq_f32 (ax, vdivq_f32 (vmulq_f32 (ax, ax), t));
> +
> + /* For large inputs (x > 2^64), asinh(x) ≈ ln(2x).
> + 1 becomes negligible in sqrt(x^2+1), so we compute
> + asinh(x) as ln(x) + ln(2). */
> + float32x4_t xy = vbslq_f32 (special, ax, y);
> + float32x4_t log_xy = log1pf_inline (xy, &d->log1pf_consts);
> +
> + /* Infinity and NaNs are the only other special cases that need checking
> + before we return the values. 0 is handled by inline_asinhf as it returns
> + 0. Below we implememt the logic that returns infinity when infinity is
> + passed and NaNs when NaNs are passed. Sometimes due to intrinsics like
> + vdivq_f32, infinity can change into NaNs, so we want to make sure the
> + right result is returned.
> +
> + Since these steps run in parallel with the log, we select between the
> + results we'll want to add to log_x in time, since addition with infinity
> + and NaNs doesn't make a difference. When we get to the adition step,
> + everything is already in the right place. */
> + float32x4_t ln2 = d->log1pf_consts.ln2;
> + uint32x4_t is_finite = vcltq_f32 (ax, d->inf);
> + float32x4_t ln2_inf_nan = vbslq_f32 (is_finite, ln2, ax);
> +
> + /* Before returning the result, the right sign will be assinged to the
> + absolute result. This is because we pass an absoulte x to the function.
> + */
> +
> + float32x4_t asinhf
> + = vbslq_f32 (special, vaddq_f32 (log_xy, ln2_inf_nan), log_xy);
> + return vreinterpretq_f32_u32 (
> + veorq_u32 (sign, vreinterpretq_u32_f32 (asinhf)));
> }
>
> /* Single-precision implementation of vector asinh(x), using vector log1p.
> @@ -99,5 +107,6 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F1 (asinh) (float32x4_t x)
> return special_case (ax, sign, special, d);
> return inline_asinhf (ax, sign, d);
> }
> +
> libmvec_hidden_def (V_NAME_F1 (asinh))
> HALF_WIDTH_ALIAS_F1 (asinh)
> diff --git a/sysdeps/aarch64/fpu/atanh_advsimd.c b/sysdeps/aarch64/fpu/atanh_advsimd.c
> index 711a234b41..a3798efc34 100644
> --- a/sysdeps/aarch64/fpu/atanh_advsimd.c
> +++ b/sysdeps/aarch64/fpu/atanh_advsimd.c
> @@ -23,38 +23,43 @@
> const static struct data
> {
> struct v_log1p_data log1p_consts;
> - uint64x2_t one;
> + float64x2_t one;
> uint64x2_t sign_mask;
> } data = { .log1p_consts = V_LOG1P_CONSTANTS_TABLE,
> - .one = V2 (0x3ff0000000000000),
> + .one = V2 (1.0),
> .sign_mask = V2 (0x8000000000000000) };
>
> static float64x2_t VPCS_ATTR NOINLINE
> -special_case (float64x2_t x, float64x2_t halfsign, float64x2_t y,
> +special_case (float64x2_t ax, float64x2_t halfsign, float64x2_t y,
> uint64x2_t special, const struct data *d)
> {
> y = log1p_inline (y, &d->log1p_consts);
> - return v_call_f64 (atanh, vbslq_f64 (d->sign_mask, halfsign, x),
> - vmulq_f64 (halfsign, y), special);
> + /* Copy the sign bit from the input to inf. */
> + uint64x2_t is_one = vceqq_f64 (ax, d->one);
> + float64x2_t signed_inf
> + = vbslq_f64 (d->sign_mask, halfsign, v_f64 (INFINITY));
> + /* Here we check for all the rest of the cases where |x| > 1 will return a
> + NaN, including if x = NaN. */
> + float64x2_t res = vbslq_f64 (special, v_f64 (NAN), vmulq_f64 (halfsign, y));
> +
> + return vbslq_f64 (is_one, signed_inf, res);
> }
>
> /* Approximation for vector double-precision atanh(x) using modified log1p.
> - The greatest observed error is 3.31 ULP:
> + The greatest observed error is 2.81 + 0.5 ULP:
> _ZGVnN2v_atanh(0x1.ffae6288b601p-6) got 0x1.ffd8ff31b5019p-6
> want 0x1.ffd8ff31b501cp-6. */
> -VPCS_ATTR
> -float64x2_t V_NAME_D1 (atanh) (float64x2_t x)
> +VPCS_ATTR float64x2_t V_NAME_D1 (atanh) (float64x2_t x)
> {
> const struct data *d = ptr_barrier (&data);
>
> float64x2_t halfsign = vbslq_f64 (d->sign_mask, x, v_f64 (0.5));
> float64x2_t ax = vabsq_f64 (x);
> - uint64x2_t ia = vreinterpretq_u64_f64 (ax);
> - uint64x2_t special = vcgeq_u64 (ia, d->one);
> + uint64x2_t special = vcageq_f64 (x, d->one);
>
> float64x2_t y;
> y = vaddq_f64 (ax, ax);
> - y = vdivq_f64 (y, vsubq_f64 (vreinterpretq_f64_u64 (d->one), ax));
> + y = vdivq_f64 (y, vsubq_f64 (d->one, ax));
>
> if (__glibc_unlikely (v_any_u64 (special)))
> return special_case (ax, halfsign, y, special, d);
> diff --git a/sysdeps/aarch64/fpu/atanh_sve.c b/sysdeps/aarch64/fpu/atanh_sve.c
> index fddf4c8706..558984974d 100644
> --- a/sysdeps/aarch64/fpu/atanh_sve.c
> +++ b/sysdeps/aarch64/fpu/atanh_sve.c
> @@ -20,13 +20,20 @@
> #define WANT_SV_LOG1P_K0_SHORTCUT 0
> #include "sv_log1p_inline.h"
>
> -#define One (0x3ff0000000000000)
> -#define Half (0x3fe0000000000000)
> +static const struct data
> +{
> + uint64_t half;
> + double inf;
> + double nan;
> +} data = { .half = 0x3fe0000000000000, .inf = INFINITY, .nan = NAN };
>
> static svfloat64_t NOINLINE
> -special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
> +special_case (svfloat64_t ax, svfloat64_t y, svbool_t pg, svbool_t special,
> + svfloat64_t halfsign, const struct data *d)
> {
> - return sv_call_f64 (atanh, x, y, special);
> + svfloat64_t res = svsel (special, sv_f64 (d->nan), y);
> + res = svsel (svcmpeq (special, ax, sv_f64 (1.0)), sv_f64 (d->inf), res);
> + return svmul_x (pg, res, halfsign);
> }
>
> /* SVE approximation for double-precision atanh, based on log1p.
> @@ -35,24 +42,25 @@ special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
> want 0x1.ffd8ff31b501cp-6. */
> svfloat64_t SV_NAME_D1 (atanh) (svfloat64_t x, const svbool_t pg)
> {
> + const struct data *d = ptr_barrier (&data);
>
> svfloat64_t ax = svabs_x (pg, x);
> svuint64_t iax = svreinterpret_u64 (ax);
> svuint64_t sign = sveor_x (pg, svreinterpret_u64 (x), iax);
> - svfloat64_t halfsign = svreinterpret_f64 (svorr_x (pg, sign, Half));
> + svfloat64_t halfsign = svreinterpret_f64 (svorr_x (pg, sign, d->half));
>
> /* It is special if iax >= 1. */
> - svbool_t special = svacge (pg, x, 1.0);
> + svbool_t special = svacge (pg, ax, 1.0);
>
> /* Computation is performed based on the following sequence of equality:
> (1+x)/(1-x) = 1 + 2x/(1-x). */
> svfloat64_t y;
> y = svadd_x (pg, ax, ax);
> - y = svdiv_x (pg, y, svsub_x (pg, sv_f64 (1), ax));
> + y = svdiv_x (pg, y, svsub_x (pg, sv_f64 (1.0), ax));
> /* ln((1+x)/(1-x)) = ln(1+2x/(1-x)) = ln(1 + y). */
> y = sv_log1p_inline (y, pg);
>
> if (__glibc_unlikely (svptest_any (pg, special)))
> - return special_case (x, svmul_x (pg, halfsign, y), special);
> + return special_case (ax, y, pg, special, halfsign, d);
> return svmul_x (pg, halfsign, y);
> }
More information about the Libc-alpha
mailing list