1 /* 2 * Single-precision vector e^x function. 3 * 4 * Copyright (c) 2019, Arm Limited. 5 * SPDX-License-Identifier: MIT 6 */ 7 8 #include "mathlib.h" 9 #include "v_math.h" 10 #if V_SUPPORTED 11 12 static const float Poly[] = { 13 /* maxerr: 0.36565 +0.5 ulp. */ 14 0x1.6a6000p-10f, 15 0x1.12718ep-7f, 16 0x1.555af0p-5f, 17 0x1.555430p-3f, 18 0x1.fffff4p-2f, 19 }; 20 #define C0 v_f32 (Poly[0]) 21 #define C1 v_f32 (Poly[1]) 22 #define C2 v_f32 (Poly[2]) 23 #define C3 v_f32 (Poly[3]) 24 #define C4 v_f32 (Poly[4]) 25 26 #define Shift v_f32 (0x1.8p23f) 27 #define InvLn2 v_f32 (0x1.715476p+0f) 28 #define Ln2hi v_f32 (0x1.62e4p-1f) 29 #define Ln2lo v_f32 (0x1.7f7d1cp-20f) 30 31 VPCS_ATTR 32 static v_f32_t 33 specialcase (v_f32_t poly, v_f32_t n, v_u32_t e, v_f32_t absn) 34 { 35 /* 2^n may overflow, break it up into s1*s2. */ 36 v_u32_t b = v_cond_u32 (n <= v_f32 (0.0f)) & v_u32 (0x83000000); 37 v_f32_t s1 = v_as_f32_u32 (v_u32 (0x7f000000) + b); 38 v_f32_t s2 = v_as_f32_u32 (e - b); 39 v_u32_t cmp = v_cond_u32 (absn > v_f32 (192.0f)); 40 v_f32_t r1 = s1 * s1; 41 v_f32_t r0 = poly * s1 * s2; 42 return v_as_f32_u32 ((cmp & v_as_u32_f32 (r1)) | (~cmp & v_as_u32_f32 (r0))); 43 } 44 45 VPCS_ATTR 46 v_f32_t 47 V_NAME(expf_1u) (v_f32_t x) 48 { 49 v_f32_t n, r, scale, poly, absn, z; 50 v_u32_t cmp, e; 51 52 /* exp(x) = 2^n * poly(r), with poly(r) in [1/sqrt(2),sqrt(2)] 53 x = ln2*n + r, with r in [-ln2/2, ln2/2]. */ 54 #if 1 55 z = v_fma_f32 (x, InvLn2, Shift); 56 n = z - Shift; 57 r = v_fma_f32 (n, -Ln2hi, x); 58 r = v_fma_f32 (n, -Ln2lo, r); 59 e = v_as_u32_f32 (z) << 23; 60 #else 61 z = x * InvLn2; 62 n = v_round_f32 (z); 63 r = v_fma_f32 (n, -Ln2hi, x); 64 r = v_fma_f32 (n, -Ln2lo, r); 65 e = v_as_u32_s32 (v_round_s32 (z)) << 23; 66 #endif 67 scale = v_as_f32_u32 (e + v_u32 (0x3f800000)); 68 absn = v_abs_f32 (n); 69 cmp = v_cond_u32 (absn > v_f32 (126.0f)); 70 poly = v_fma_f32 (C0, r, C1); 71 poly = v_fma_f32 (poly, r, C2); 72 poly = v_fma_f32 (poly, r, C3); 73 poly = v_fma_f32 (poly, r, C4); 74 poly = v_fma_f32 (poly, r, v_f32 (1.0f)); 75 poly = v_fma_f32 (poly, r, v_f32 (1.0f)); 76 if (unlikely (v_any_u32 (cmp))) 77 return specialcase (poly, n, e, absn); 78 return scale * poly; 79 } 80 #endif 81