1 /*
2 * Single-precision vector exp(x) - 1 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 "v_math.h"
9 #include "poly_advsimd_f32.h"
10 #include "pl_sig.h"
11 #include "pl_test.h"
12
13 static const struct data
14 {
15 float32x4_t poly[5];
16 float invln2_and_ln2[4];
17 float32x4_t shift;
18 int32x4_t exponent_bias;
19 #if WANT_SIMD_EXCEPT
20 uint32x4_t thresh;
21 #else
22 float32x4_t oflow_bound;
23 #endif
24 } data = {
25 /* Generated using fpminimax with degree=5 in [-log(2)/2, log(2)/2]. */
26 .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5),
27 V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) },
28 /* Stores constants: invln2, ln2_hi, ln2_lo, 0. */
29 .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 },
30 .shift = V4 (0x1.8p23f),
31 .exponent_bias = V4 (0x3f800000),
32 #if !WANT_SIMD_EXCEPT
33 /* Value above which expm1f(x) should overflow. Absolute value of the
34 underflow bound is greater than this, so it catches both cases - there is
35 a small window where fallbacks are triggered unnecessarily. */
36 .oflow_bound = V4 (0x1.5ebc4p+6),
37 #else
38 /* asuint(oflow_bound) - asuint(0x1p-23), shifted left by 1 for absolute
39 compare. */
40 .thresh = V4 (0x1d5ebc40),
41 #endif
42 };
43
44 /* asuint(0x1p-23), shifted by 1 for abs compare. */
45 #define TinyBound v_u32 (0x34000000 << 1)
46
47 static float32x4_t VPCS_ATTR NOINLINE
special_case(float32x4_t x,float32x4_t y,uint32x4_t special)48 special_case (float32x4_t x, float32x4_t y, uint32x4_t special)
49 {
50 return v_call_f32 (expm1f, x, y, special);
51 }
52
53 /* Single-precision vector exp(x) - 1 function.
54 The maximum error is 1.51 ULP:
55 _ZGVnN4v_expm1f (0x1.8baa96p-2) got 0x1.e2fb9p-2
56 want 0x1.e2fb94p-2. */
V_NAME_F1(expm1)57 float32x4_t VPCS_ATTR V_NAME_F1 (expm1) (float32x4_t x)
58 {
59 const struct data *d = ptr_barrier (&data);
60 uint32x4_t ix = vreinterpretq_u32_f32 (x);
61
62 #if WANT_SIMD_EXCEPT
63 /* If fp exceptions are to be triggered correctly, fall back to scalar for
64 |x| < 2^-23, |x| > oflow_bound, Inf & NaN. Add ix to itself for
65 shift-left by 1, and compare with thresh which was left-shifted offline -
66 this is effectively an absolute compare. */
67 uint32x4_t special
68 = vcgeq_u32 (vsubq_u32 (vaddq_u32 (ix, ix), TinyBound), d->thresh);
69 if (unlikely (v_any_u32 (special)))
70 x = v_zerofy_f32 (x, special);
71 #else
72 /* Handles very large values (+ve and -ve), +/-NaN, +/-Inf. */
73 uint32x4_t special = vcagtq_f32 (x, d->oflow_bound);
74 #endif
75
76 /* Reduce argument to smaller range:
77 Let i = round(x / ln2)
78 and f = x - i * ln2, then f is in [-ln2/2, ln2/2].
79 exp(x) - 1 = 2^i * (expm1(f) + 1) - 1
80 where 2^i is exact because i is an integer. */
81 float32x4_t invln2_and_ln2 = vld1q_f32 (d->invln2_and_ln2);
82 float32x4_t j
83 = vsubq_f32 (vfmaq_laneq_f32 (d->shift, x, invln2_and_ln2, 0), d->shift);
84 int32x4_t i = vcvtq_s32_f32 (j);
85 float32x4_t f = vfmsq_laneq_f32 (x, j, invln2_and_ln2, 1);
86 f = vfmsq_laneq_f32 (f, j, invln2_and_ln2, 2);
87
88 /* Approximate expm1(f) using polynomial.
89 Taylor expansion for expm1(x) has the form:
90 x + ax^2 + bx^3 + cx^4 ....
91 So we calculate the polynomial P(f) = a + bf + cf^2 + ...
92 and assemble the approximation expm1(f) ~= f + f^2 * P(f). */
93 float32x4_t p = v_horner_4_f32 (f, d->poly);
94 p = vfmaq_f32 (f, vmulq_f32 (f, f), p);
95
96 /* Assemble the result.
97 expm1(x) ~= 2^i * (p + 1) - 1
98 Let t = 2^i. */
99 int32x4_t u = vaddq_s32 (vshlq_n_s32 (i, 23), d->exponent_bias);
100 float32x4_t t = vreinterpretq_f32_s32 (u);
101
102 if (unlikely (v_any_u32 (special)))
103 return special_case (vreinterpretq_f32_u32 (ix),
104 vfmaq_f32 (vsubq_f32 (t, v_f32 (1.0f)), p, t),
105 special);
106
107 /* expm1(x) ~= p * t + (t - 1). */
108 return vfmaq_f32 (vsubq_f32 (t, v_f32 (1.0f)), p, t);
109 }
110
111 PL_SIG (V, F, 1, expm1, -9.9, 9.9)
112 PL_TEST_ULP (V_NAME_F1 (expm1), 1.02)
113 PL_TEST_EXPECT_FENV (V_NAME_F1 (expm1), WANT_SIMD_EXCEPT)
114 PL_TEST_SYM_INTERVAL (V_NAME_F1 (expm1), 0, 0x1p-23, 1000)
115 PL_TEST_INTERVAL (V_NAME_F1 (expm1), -0x1p-23, 0x1.5ebc4p+6, 1000000)
116 PL_TEST_INTERVAL (V_NAME_F1 (expm1), -0x1p-23, -0x1.9bbabcp+6, 1000000)
117 PL_TEST_INTERVAL (V_NAME_F1 (expm1), 0x1.5ebc4p+6, inf, 1000)
118 PL_TEST_INTERVAL (V_NAME_F1 (expm1), -0x1.9bbabcp+6, -inf, 1000)
119