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