1 /* 2 * Single-precision vector e^x function. 3 * 4 * Copyright (c) 2019-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "sv_math.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 #if SV_SUPPORTED 13 14 #define C(i) __sv_expf_poly[i] 15 16 #define InvLn2 (0x1.715476p+0f) 17 #define Ln2hi (0x1.62e4p-1f) 18 #define Ln2lo (0x1.7f7d1cp-20f) 19 20 #if SV_EXPF_USE_FEXPA 21 22 #define Shift (0x1.903f8p17f) /* 1.5*2^17 + 127. */ 23 #define Thres \ 24 (0x1.5d5e2ap+6f) /* Roughly 87.3. For x < -Thres, the result is subnormal \ 25 and not handled correctly by FEXPA. */ 26 27 static NOINLINE sv_f32_t 28 special_case (sv_f32_t x, sv_f32_t y, svbool_t special) 29 { 30 /* The special-case handler from the Neon routine does not handle subnormals 31 in a way that is compatible with FEXPA. For the FEXPA variant we just fall 32 back to scalar expf. */ 33 return sv_call_f32 (expf, x, y, special); 34 } 35 36 #else 37 38 #define Shift (0x1.8p23f) /* 1.5 * 2^23. */ 39 #define Thres (126.0f) 40 41 /* Special-case handler adapted from Neon variant. Uses s, y and n to produce 42 the final result (normal cases included). It performs an update of all lanes! 43 Therefore: 44 - all previous computation need to be done on all lanes indicated by input 45 pg 46 - we cannot simply apply the special case to the special-case-activated 47 lanes. Besides it is likely that this would not increase performance (no 48 scatter/gather). */ 49 static inline sv_f32_t 50 specialcase (svbool_t pg, sv_f32_t poly, sv_f32_t n, sv_u32_t e, 51 svbool_t p_cmp1, sv_f32_t scale) 52 { 53 /* s=2^(n/N) may overflow, break it up into s=s1*s2, 54 such that exp = s + s*y can be computed as s1*(s2+s2*y) 55 and s1*s1 overflows only if n>0. */ 56 57 /* If n<=0 then set b to 0x820...0, 0 otherwise. */ 58 svbool_t p_sign = svcmple_n_f32 (pg, n, 0.0f); /* n <= 0. */ 59 sv_u32_t b 60 = svdup_n_u32_z (p_sign, 0x82000000); /* Inactive lanes set to 0. */ 61 62 /* Set s1 to generate overflow depending on sign of exponent n. */ 63 sv_f32_t s1 64 = sv_as_f32_u32 (svadd_n_u32_x (pg, b, 0x7f000000)); /* b + 0x7f000000. */ 65 /* Offset s to avoid overflow in final result if n is below threshold. */ 66 sv_f32_t s2 = sv_as_f32_u32 ( 67 svsub_u32_x (pg, e, b)); /* as_u32 (s) - 0x3010...0 + b. */ 68 69 /* |n| > 192 => 2^(n/N) overflows. */ 70 svbool_t p_cmp2 = svacgt_n_f32 (pg, n, 192.0f); 71 72 sv_f32_t r2 = svmul_f32_x (pg, s1, s1); 73 sv_f32_t r1 = sv_fma_f32_x (pg, poly, s2, s2); 74 r1 = svmul_f32_x (pg, r1, s1); 75 sv_f32_t r0 = sv_fma_f32_x (pg, poly, scale, scale); 76 77 /* Apply condition 1 then 2. 78 Returns r2 if cond2 is true, otherwise 79 if cond1 is true then return r1, otherwise return r0. */ 80 sv_f32_t r = svsel_f32 (p_cmp1, r1, r0); 81 82 return svsel_f32 (p_cmp2, r2, r); 83 } 84 85 #endif 86 87 /* Optimised single-precision SVE exp function. By default this is an SVE port 88 of the Neon algorithm from math/. Alternatively, enable a modification of 89 that algorithm that looks up scale using SVE FEXPA instruction with 90 SV_EXPF_USE_FEXPA. 91 92 Worst-case error of the default algorithm is 1.95 ulp: 93 __sv_expf(-0x1.4cb74ap+2) got 0x1.6a022cp-8 94 want 0x1.6a023p-8. 95 96 Worst-case error when using FEXPA is 1.04 ulp: 97 __sv_expf(0x1.a8eda4p+1) got 0x1.ba74bcp+4 98 want 0x1.ba74bap+4. */ 99 sv_f32_t 100 __sv_expf_x (sv_f32_t x, const svbool_t pg) 101 { 102 /* exp(x) = 2^n (1 + poly(r)), with 1 + poly(r) in [1/sqrt(2),sqrt(2)] 103 x = ln2*n + r, with r in [-ln2/2, ln2/2]. */ 104 105 /* n = round(x/(ln2/N)). */ 106 sv_f32_t z = sv_fma_n_f32_x (pg, InvLn2, x, sv_f32 (Shift)); 107 sv_f32_t n = svsub_n_f32_x (pg, z, Shift); 108 109 /* r = x - n*ln2/N. */ 110 sv_f32_t r = sv_fma_n_f32_x (pg, -Ln2hi, n, x); 111 r = sv_fma_n_f32_x (pg, -Ln2lo, n, r); 112 113 /* scale = 2^(n/N). */ 114 #if SV_EXPF_USE_FEXPA 115 /* NaNs also need special handling with FEXPA. */ 116 svbool_t is_special_case 117 = svorr_b_z (pg, svacgt_n_f32 (pg, x, Thres), svcmpne_f32 (pg, x, x)); 118 sv_f32_t scale = svexpa_f32 (sv_as_u32_f32 (z)); 119 #else 120 sv_u32_t e = svlsl_n_u32_x (pg, sv_as_u32_f32 (z), 23); 121 svbool_t is_special_case = svacgt_n_f32 (pg, n, Thres); 122 sv_f32_t scale = sv_as_f32_u32 (svadd_n_u32_x (pg, e, 0x3f800000)); 123 #endif 124 125 /* y = exp(r) - 1 ~= r + C1 r^2 + C2 r^3 + C3 r^4. */ 126 sv_f32_t r2 = svmul_f32_x (pg, r, r); 127 sv_f32_t p = sv_fma_n_f32_x (pg, C (0), r, sv_f32 (C (1))); 128 sv_f32_t q = sv_fma_n_f32_x (pg, C (2), r, sv_f32 (C (3))); 129 q = sv_fma_f32_x (pg, p, r2, q); 130 p = svmul_n_f32_x (pg, r, C (4)); 131 sv_f32_t poly = sv_fma_f32_x (pg, q, r2, p); 132 133 if (unlikely (svptest_any (pg, is_special_case))) 134 #if SV_EXPF_USE_FEXPA 135 return special_case (x, sv_fma_f32_x (pg, poly, scale, scale), 136 is_special_case); 137 #else 138 return specialcase (pg, poly, n, e, is_special_case, scale); 139 #endif 140 141 return sv_fma_f32_x (pg, poly, scale, scale); 142 } 143 144 PL_ALIAS (__sv_expf_x, _ZGVsMxv_expf) 145 146 PL_SIG (SV, F, 1, exp, -9.9, 9.9) 147 PL_TEST_ULP (__sv_expf, 1.46) 148 PL_TEST_INTERVAL (__sv_expf, 0, 0x1p-23, 40000) 149 PL_TEST_INTERVAL (__sv_expf, 0x1p-23, 1, 50000) 150 PL_TEST_INTERVAL (__sv_expf, 1, 0x1p23, 50000) 151 PL_TEST_INTERVAL (__sv_expf, 0x1p23, inf, 50000) 152 PL_TEST_INTERVAL (__sv_expf, -0, -0x1p-23, 40000) 153 PL_TEST_INTERVAL (__sv_expf, -0x1p-23, -1, 50000) 154 PL_TEST_INTERVAL (__sv_expf, -1, -0x1p23, 50000) 155 PL_TEST_INTERVAL (__sv_expf, -0x1p23, -inf, 50000) 156 #endif // SV_SUPPORTED 157