1 /*
2 * Double-precision vector e^(x+tail) function.
3 *
4 * Copyright (c) 2019-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7 #ifndef PL_MATH_V_EXP_TAIL_INLINE_H
8 #define PL_MATH_V_EXP_TAIL_INLINE_H
9
10 #include "v_math.h"
11 #include "poly_advsimd_f64.h"
12
13 #ifndef WANT_V_EXP_TAIL_SPECIALCASE
14 #error \
15 "Cannot use v_exp_tail_inline.h without specifying whether you need the special case computation."
16 #endif
17
18 #define N (1 << V_EXP_TAIL_TABLE_BITS)
19
20 static const struct data
21 {
22 float64x2_t poly[4];
23 #if WANT_V_EXP_TAIL_SPECIALCASE
24 float64x2_t big_bound, huge_bound;
25 #endif
26 float64x2_t shift, invln2, ln2_hi, ln2_lo;
27 } data = {
28 #if WANT_V_EXP_TAIL_SPECIALCASE
29 .big_bound = V2 (704.0),
30 .huge_bound = V2 (1280.0 * N),
31 #endif
32 .shift = V2 (0x1.8p52),
33 .invln2 = V2 (0x1.71547652b82fep8), /* N/ln2. */
34 .ln2_hi = V2 (0x1.62e42fefa39efp-9), /* ln2/N. */
35 .ln2_lo = V2 (0x1.abc9e3b39803f3p-64),
36 .poly = { V2 (1.0), V2 (0x1.fffffffffffd4p-2), V2 (0x1.5555571d6b68cp-3),
37 V2 (0x1.5555576a59599p-5) },
38 };
39
40 static inline uint64x2_t
lookup_sbits(uint64x2_t i)41 lookup_sbits (uint64x2_t i)
42 {
43 return (uint64x2_t){__v_exp_tail_data[i[0]], __v_exp_tail_data[i[1]]};
44 }
45
46 #if WANT_V_EXP_TAIL_SPECIALCASE
47 #define SpecialOffset v_u64 (0x6000000000000000) /* 0x1p513. */
48 /* The following 2 bias when combined form the exponent bias:
49 SpecialBias1 - SpecialBias2 = asuint64(1.0). */
50 #define SpecialBias1 v_u64 (0x7000000000000000) /* 0x1p769. */
51 #define SpecialBias2 v_u64 (0x3010000000000000) /* 0x1p-254. */
52 static float64x2_t VPCS_ATTR
v_exp_tail_special_case(float64x2_t s,float64x2_t y,float64x2_t n,const struct data * d)53 v_exp_tail_special_case (float64x2_t s, float64x2_t y, float64x2_t n,
54 const struct data *d)
55 {
56 /* 2^(n/N) may overflow, break it up into s1*s2. */
57 uint64x2_t b = vandq_u64 (vclezq_f64 (n), SpecialOffset);
58 float64x2_t s1 = vreinterpretq_f64_u64 (vsubq_u64 (SpecialBias1, b));
59 float64x2_t s2 = vreinterpretq_f64_u64 (
60 vaddq_u64 (vsubq_u64 (vreinterpretq_u64_f64 (s), SpecialBias2), b));
61 uint64x2_t oflow = vcagtq_f64 (n, d->huge_bound);
62 float64x2_t r0 = vmulq_f64 (vfmaq_f64 (s2, y, s2), s1);
63 float64x2_t r1 = vmulq_f64 (s1, s1);
64 return vbslq_f64 (oflow, r1, r0);
65 }
66 #endif
67
68 static inline float64x2_t VPCS_ATTR
v_exp_tail_inline(float64x2_t x,float64x2_t xtail)69 v_exp_tail_inline (float64x2_t x, float64x2_t xtail)
70 {
71 const struct data *d = ptr_barrier (&data);
72 #if WANT_V_EXP_TAIL_SPECIALCASE
73 uint64x2_t special = vcgtq_f64 (vabsq_f64 (x), d->big_bound);
74 #endif
75 /* n = round(x/(ln2/N)). */
76 float64x2_t z = vfmaq_f64 (d->shift, x, d->invln2);
77 uint64x2_t u = vreinterpretq_u64_f64 (z);
78 float64x2_t n = vsubq_f64 (z, d->shift);
79
80 /* r = x - n*ln2/N. */
81 float64x2_t r = x;
82 r = vfmsq_f64 (r, d->ln2_hi, n);
83 r = vfmsq_f64 (r, d->ln2_lo, n);
84
85 uint64x2_t e = vshlq_n_u64 (u, 52 - V_EXP_TAIL_TABLE_BITS);
86 uint64x2_t i = vandq_u64 (u, v_u64 (N - 1));
87
88 /* y = tail + exp(r) - 1 ~= r + C1 r^2 + C2 r^3 + C3 r^4, using Horner. */
89 float64x2_t y = v_horner_3_f64 (r, d->poly);
90 y = vfmaq_f64 (xtail, y, r);
91
92 /* s = 2^(n/N). */
93 u = lookup_sbits (i);
94 float64x2_t s = vreinterpretq_f64_u64 (vaddq_u64 (u, e));
95
96 #if WANT_V_EXP_TAIL_SPECIALCASE
97 if (unlikely (v_any_u64 (special)))
98 return v_exp_tail_special_case (s, y, n, d);
99 #endif
100 return vfmaq_f64 (s, y, s);
101 }
102 #endif // PL_MATH_V_EXP_TAIL_INLINE_H
103