1 | /* Double-precision vector (SVE) exp function. |
2 | |
3 | Copyright (C) 2023-2024 Free Software Foundation, Inc. |
4 | This file is part of the GNU C Library. |
5 | |
6 | The GNU C Library is free software; you can redistribute it and/or |
7 | modify it under the terms of the GNU Lesser General Public |
8 | License as published by the Free Software Foundation; either |
9 | version 2.1 of the License, or (at your option) any later version. |
10 | |
11 | The GNU C Library is distributed in the hope that it will be useful, |
12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
14 | Lesser General Public License for more details. |
15 | |
16 | You should have received a copy of the GNU Lesser General Public |
17 | License along with the GNU C Library; if not, see |
18 | <https://www.gnu.org/licenses/>. */ |
19 | |
20 | #include "sv_math.h" |
21 | |
22 | static const struct data |
23 | { |
24 | double poly[4]; |
25 | double ln2_hi, ln2_lo, inv_ln2, shift, thres; |
26 | } data = { |
27 | .poly = { /* ulp error: 0.53. */ |
28 | 0x1.fffffffffdbcdp-2, 0x1.555555555444cp-3, 0x1.555573c6a9f7dp-5, |
29 | 0x1.1111266d28935p-7 }, |
30 | .ln2_hi = 0x1.62e42fefa3800p-1, |
31 | .ln2_lo = 0x1.ef35793c76730p-45, |
32 | /* 1/ln2. */ |
33 | .inv_ln2 = 0x1.71547652b82fep+0, |
34 | /* 1.5*2^46+1023. This value is further explained below. */ |
35 | .shift = 0x1.800000000ffc0p+46, |
36 | .thres = 704.0, |
37 | }; |
38 | |
39 | #define C(i) sv_f64 (d->poly[i]) |
40 | #define SpecialOffset 0x6000000000000000 /* 0x1p513. */ |
41 | /* SpecialBias1 + SpecialBias1 = asuint(1.0). */ |
42 | #define SpecialBias1 0x7000000000000000 /* 0x1p769. */ |
43 | #define SpecialBias2 0x3010000000000000 /* 0x1p-254. */ |
44 | |
45 | /* Update of both special and non-special cases, if any special case is |
46 | detected. */ |
47 | static inline svfloat64_t |
48 | special_case (svbool_t pg, svfloat64_t s, svfloat64_t y, svfloat64_t n) |
49 | { |
50 | /* s=2^n may overflow, break it up into s=s1*s2, |
51 | such that exp = s + s*y can be computed as s1*(s2+s2*y) |
52 | and s1*s1 overflows only if n>0. */ |
53 | |
54 | /* If n<=0 then set b to 0x6, 0 otherwise. */ |
55 | svbool_t p_sign = svcmple (pg, n, 0.0); /* n <= 0. */ |
56 | svuint64_t b |
57 | = svdup_u64_z (p_sign, SpecialOffset); /* Inactive lanes set to 0. */ |
58 | |
59 | /* Set s1 to generate overflow depending on sign of exponent n. */ |
60 | svfloat64_t s1 = svreinterpret_f64 ( |
61 | svsubr_x (pg, b, SpecialBias1)); /* 0x70...0 - b. */ |
62 | /* Offset s to avoid overflow in final result if n is below threshold. */ |
63 | svfloat64_t s2 = svreinterpret_f64 ( |
64 | svadd_x (pg, svsub_x (pg, svreinterpret_u64 (s), SpecialBias2), |
65 | b)); /* as_u64 (s) - 0x3010...0 + b. */ |
66 | |
67 | /* |n| > 1280 => 2^(n) overflows. */ |
68 | svbool_t p_cmp = svacgt (pg, n, 1280.0); |
69 | |
70 | svfloat64_t r1 = svmul_x (pg, s1, s1); |
71 | svfloat64_t r2 = svmla_x (pg, s2, s2, y); |
72 | svfloat64_t r0 = svmul_x (pg, r2, s1); |
73 | |
74 | return svsel (p_cmp, r1, r0); |
75 | } |
76 | |
77 | /* SVE exp algorithm. Maximum measured error is 1.01ulps: |
78 | SV_NAME_D1 (exp)(0x1.4619d7b04da41p+6) got 0x1.885d9acc41da7p+117 |
79 | want 0x1.885d9acc41da6p+117. */ |
80 | svfloat64_t SV_NAME_D1 (exp) (svfloat64_t x, const svbool_t pg) |
81 | { |
82 | const struct data *d = ptr_barrier (&data); |
83 | |
84 | svbool_t special = svacgt (pg, x, d->thres); |
85 | |
86 | /* Use a modifed version of the shift used for flooring, such that x/ln2 is |
87 | rounded to a multiple of 2^-6=1/64, shift = 1.5 * 2^52 * 2^-6 = 1.5 * |
88 | 2^46. |
89 | |
90 | n is not an integer but can be written as n = m + i/64, with i and m |
91 | integer, 0 <= i < 64 and m <= n. |
92 | |
93 | Bits 5:0 of z will be null every time x/ln2 reaches a new integer value |
94 | (n=m, i=0), and is incremented every time z (or n) is incremented by 1/64. |
95 | FEXPA expects i in bits 5:0 of the input so it can be used as index into |
96 | FEXPA hardwired table T[i] = 2^(i/64) for i = 0:63, that will in turn |
97 | populate the mantissa of the output. Therefore, we use u=asuint(z) as |
98 | input to FEXPA. |
99 | |
100 | We add 1023 to the modified shift value in order to set bits 16:6 of u to |
101 | 1, such that once these bits are moved to the exponent of the output of |
102 | FEXPA, we get the exponent of 2^n right, i.e. we get 2^m. */ |
103 | svfloat64_t z = svmla_x (pg, sv_f64 (x: d->shift), x, d->inv_ln2); |
104 | svuint64_t u = svreinterpret_u64 (z); |
105 | svfloat64_t n = svsub_x (pg, z, d->shift); |
106 | |
107 | /* r = x - n * ln2, r is in [-ln2/(2N), ln2/(2N)]. */ |
108 | svfloat64_t ln2 = svld1rq (svptrue_b64 (), &d->ln2_hi); |
109 | svfloat64_t r = svmls_lane (x, n, ln2, 0); |
110 | r = svmls_lane (r, n, ln2, 1); |
111 | |
112 | /* y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5. */ |
113 | svfloat64_t r2 = svmul_x (pg, r, r); |
114 | svfloat64_t p01 = svmla_x (pg, C (0), C (1), r); |
115 | svfloat64_t p23 = svmla_x (pg, C (2), C (3), r); |
116 | svfloat64_t p04 = svmla_x (pg, p01, p23, r2); |
117 | svfloat64_t y = svmla_x (pg, r, p04, r2); |
118 | |
119 | /* s = 2^n, computed using FEXPA. FEXPA does not propagate NaNs, so for |
120 | consistent NaN handling we have to manually propagate them. This comes at |
121 | significant performance cost. */ |
122 | svfloat64_t s = svexpa (u); |
123 | |
124 | /* Assemble result as exp(x) = 2^n * exp(r). If |x| > Thresh the |
125 | multiplication may overflow, so use special case routine. */ |
126 | |
127 | if (__glibc_unlikely (svptest_any (pg, special))) |
128 | { |
129 | /* FEXPA zeroes the sign bit, however the sign is meaningful to the |
130 | special case function so needs to be copied. |
131 | e = sign bit of u << 46. */ |
132 | svuint64_t e = svand_x (pg, svlsl_x (pg, u, 46), 0x8000000000000000); |
133 | /* Copy sign to s. */ |
134 | s = svreinterpret_f64 (svadd_x (pg, e, svreinterpret_u64 (s))); |
135 | return special_case (pg, s, y, n); |
136 | } |
137 | |
138 | /* No special case. */ |
139 | return svmla_x (pg, s, s, y); |
140 | } |
141 | |