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 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 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 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