1 /* 2 * Helper for single-precision routines which calculate exp(x) - 1 and do not 3 * need special-case handling 4 * 5 * Copyright (c) 2022-2023, Arm Limited. 6 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 7 */ 8 9 #ifndef PL_MATH_V_EXPM1F_INLINE_H 10 #define PL_MATH_V_EXPM1F_INLINE_H 11 12 #include "v_math.h" 13 #include "math_config.h" 14 #include "estrinf.h" 15 16 #define One 0x3f800000 17 #define Shift v_f32 (0x1.8p23f) 18 #define InvLn2 v_f32 (0x1.715476p+0f) 19 #define MLn2hi v_f32 (-0x1.62e4p-1f) 20 #define MLn2lo v_f32 (-0x1.7f7d1cp-20f) 21 22 #define C(i) v_f32 (__expm1f_poly[i]) 23 24 static inline v_f32_t 25 expm1f_inline (v_f32_t x) 26 { 27 /* Helper routine for calculating exp(x) - 1. 28 Copied from v_expm1f_1u6.c, with all special-case handling removed - the 29 calling routine should handle special values if required. */ 30 31 /* Reduce argument: f in [-ln2/2, ln2/2], i is exact. */ 32 v_f32_t j = v_fma_f32 (InvLn2, x, Shift) - Shift; 33 v_s32_t i = v_to_s32_f32 (j); 34 v_f32_t f = v_fma_f32 (j, MLn2hi, x); 35 f = v_fma_f32 (j, MLn2lo, f); 36 37 /* Approximate expm1(f) with polynomial P, expm1(f) ~= f + f^2 * P(f). 38 Uses Estrin scheme, where the main __v_expm1f routine uses Horner. */ 39 v_f32_t f2 = f * f; 40 v_f32_t p = ESTRIN_4 (f, f2, f2 * f2, C); 41 p = v_fma_f32 (f2, p, f); 42 43 /* t = 2^i. */ 44 v_f32_t t = v_as_f32_u32 (v_as_u32_s32 (i << 23) + One); 45 /* expm1(x) ~= p * t + (t - 1). */ 46 return v_fma_f32 (p, t, t - 1); 47 } 48 49 #endif // PL_MATH_V_EXPM1F_INLINE_H 50