1 /*
2 * Double-precision SVE asinh(x) function.
3 *
4 * Copyright (c) 2022-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "sv_math.h"
9 #include "poly_sve_f64.h"
10 #include "pl_sig.h"
11 #include "pl_test.h"
12
13 #define SignMask (0x8000000000000000)
14 #define One (0x3ff0000000000000)
15 #define Thres (0x5fe0000000000000) /* asuint64 (0x1p511). */
16
17 static const struct data
18 {
19 double poly[18];
20 double ln2, p3, p1, p4, p0, p2;
21 uint64_t n;
22 uint64_t off;
23
24 } data = {
25 /* Polynomial generated using Remez on [2^-26, 1]. */
26 .poly
27 = { -0x1.55555555554a7p-3, 0x1.3333333326c7p-4, -0x1.6db6db68332e6p-5,
28 0x1.f1c71b26fb40dp-6, -0x1.6e8b8b654a621p-6, 0x1.1c4daa9e67871p-6,
29 -0x1.c9871d10885afp-7, 0x1.7a16e8d9d2ecfp-7, -0x1.3ddca533e9f54p-7,
30 0x1.0becef748dafcp-7, -0x1.b90c7099dd397p-8, 0x1.541f2bb1ffe51p-8,
31 -0x1.d217026a669ecp-9, 0x1.0b5c7977aaf7p-9, -0x1.e0f37daef9127p-11,
32 0x1.388b5fe542a6p-12, -0x1.021a48685e287p-14, 0x1.93d4ba83d34dap-18 },
33 .ln2 = 0x1.62e42fefa39efp-1,
34 .p0 = -0x1.ffffffffffff7p-2,
35 .p1 = 0x1.55555555170d4p-2,
36 .p2 = -0x1.0000000399c27p-2,
37 .p3 = 0x1.999b2e90e94cap-3,
38 .p4 = -0x1.554e550bd501ep-3,
39 .n = 1 << V_LOG_TABLE_BITS,
40 .off = 0x3fe6900900000000
41 };
42
43 static svfloat64_t NOINLINE
special_case(svfloat64_t x,svfloat64_t y,svbool_t special)44 special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
45 {
46 return sv_call_f64 (asinh, x, y, special);
47 }
48
49 static inline svfloat64_t
__sv_log_inline(svfloat64_t x,const struct data * d,const svbool_t pg)50 __sv_log_inline (svfloat64_t x, const struct data *d, const svbool_t pg)
51 {
52 /* Double-precision SVE log, copied from SVE log implementation with some
53 cosmetic modification and special-cases removed. See that file for details
54 of the algorithm used. */
55
56 svuint64_t ix = svreinterpret_u64 (x);
57 svuint64_t tmp = svsub_x (pg, ix, d->off);
58 svuint64_t i = svand_x (pg, svlsr_x (pg, tmp, (51 - V_LOG_TABLE_BITS)),
59 (d->n - 1) << 1);
60 svint64_t k = svasr_x (pg, svreinterpret_s64 (tmp), 52);
61 svuint64_t iz = svsub_x (pg, ix, svand_x (pg, tmp, 0xfffULL << 52));
62 svfloat64_t z = svreinterpret_f64 (iz);
63
64 svfloat64_t invc = svld1_gather_index (pg, &__v_log_data.table[0].invc, i);
65 svfloat64_t logc = svld1_gather_index (pg, &__v_log_data.table[0].logc, i);
66
67 svfloat64_t ln2_p3 = svld1rq (svptrue_b64 (), &d->ln2);
68 svfloat64_t p1_p4 = svld1rq (svptrue_b64 (), &d->p1);
69
70 svfloat64_t r = svmla_x (pg, sv_f64 (-1.0), invc, z);
71 svfloat64_t kd = svcvt_f64_x (pg, k);
72
73 svfloat64_t hi = svmla_lane (svadd_x (pg, logc, r), kd, ln2_p3, 0);
74 svfloat64_t r2 = svmul_x (pg, r, r);
75
76 svfloat64_t y = svmla_lane (sv_f64 (d->p2), r, ln2_p3, 1);
77
78 svfloat64_t p = svmla_lane (sv_f64 (d->p0), r, p1_p4, 0);
79 y = svmla_lane (y, r2, p1_p4, 1);
80 y = svmla_x (pg, p, r2, y);
81 y = svmla_x (pg, hi, r2, y);
82 return y;
83 }
84
85 /* Double-precision implementation of SVE asinh(x).
86 asinh is very sensitive around 1, so it is impractical to devise a single
87 low-cost algorithm which is sufficiently accurate on a wide range of input.
88 Instead we use two different algorithms:
89 asinh(x) = sign(x) * log(|x| + sqrt(x^2 + 1) if |x| >= 1
90 = sign(x) * (|x| + |x|^3 * P(x^2)) otherwise
91 where log(x) is an optimized log approximation, and P(x) is a polynomial
92 shared with the scalar routine. The greatest observed error 2.51 ULP, in
93 |x| >= 1:
94 _ZGVsMxv_asinh(0x1.170469d024505p+0) got 0x1.e3181c43b0f36p-1
95 want 0x1.e3181c43b0f39p-1. */
SV_NAME_D1(asinh)96 svfloat64_t SV_NAME_D1 (asinh) (svfloat64_t x, const svbool_t pg)
97 {
98 const struct data *d = ptr_barrier (&data);
99
100 svuint64_t ix = svreinterpret_u64 (x);
101 svuint64_t iax = svbic_x (pg, ix, SignMask);
102 svuint64_t sign = svand_x (pg, ix, SignMask);
103 svfloat64_t ax = svreinterpret_f64 (iax);
104
105 svbool_t ge1 = svcmpge (pg, iax, One);
106 svbool_t special = svcmpge (pg, iax, Thres);
107
108 /* Option 1: |x| >= 1.
109 Compute asinh(x) according by asinh(x) = log(x + sqrt(x^2 + 1)). */
110 svfloat64_t option_1 = sv_f64 (0);
111 if (likely (svptest_any (pg, ge1)))
112 {
113 svfloat64_t x2 = svmul_x (pg, ax, ax);
114 option_1 = __sv_log_inline (
115 svadd_x (pg, ax, svsqrt_x (pg, svadd_x (pg, x2, 1))), d, pg);
116 }
117
118 /* Option 2: |x| < 1.
119 Compute asinh(x) using a polynomial.
120 The largest observed error in this region is 1.51 ULPs:
121 _ZGVsMxv_asinh(0x1.fe12bf8c616a2p-1) got 0x1.c1e649ee2681bp-1
122 want 0x1.c1e649ee2681dp-1. */
123 svfloat64_t option_2 = sv_f64 (0);
124 if (likely (svptest_any (pg, svnot_z (pg, ge1))))
125 {
126 svfloat64_t x2 = svmul_x (pg, ax, ax);
127 svfloat64_t x4 = svmul_x (pg, x2, x2);
128 svfloat64_t p = sv_pw_horner_17_f64_x (pg, x2, x4, d->poly);
129 option_2 = svmla_x (pg, ax, p, svmul_x (pg, x2, ax));
130 }
131
132 /* Choose the right option for each lane. */
133 svfloat64_t y = svsel (ge1, option_1, option_2);
134
135 if (unlikely (svptest_any (pg, special)))
136 return special_case (
137 x, svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (y), sign)),
138 special);
139 return svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (y), sign));
140 }
141
142 PL_SIG (SV, D, 1, asinh, -10.0, 10.0)
143 PL_TEST_ULP (SV_NAME_D1 (asinh), 2.52)
144 /* Test vector asinh 3 times, with control lane < 1, > 1 and special.
145 Ensures the svsel is choosing the right option in all cases. */
146 #define SV_ASINH_INTERVAL(lo, hi, n) \
147 PL_TEST_SYM_INTERVAL_C (SV_NAME_D1 (asinh), lo, hi, n, 0.5) \
148 PL_TEST_SYM_INTERVAL_C (SV_NAME_D1 (asinh), lo, hi, n, 2) \
149 PL_TEST_SYM_INTERVAL_C (SV_NAME_D1 (asinh), lo, hi, n, 0x1p600)
150 SV_ASINH_INTERVAL (0, 0x1p-26, 50000)
151 SV_ASINH_INTERVAL (0x1p-26, 1, 50000)
152 SV_ASINH_INTERVAL (1, 0x1p511, 50000)
153 SV_ASINH_INTERVAL (0x1p511, inf, 40000)
154