| /* |
| * Single-precision vector e^x function. |
| * |
| * Copyright (c) 2019-2023, Arm Limited. |
| * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception |
| */ |
| |
| #include "sv_math.h" |
| #include "pl_sig.h" |
| #include "pl_test.h" |
| |
| #if SV_SUPPORTED |
| |
| #define C(i) __sv_expf_poly[i] |
| |
| #define InvLn2 (0x1.715476p+0f) |
| #define Ln2hi (0x1.62e4p-1f) |
| #define Ln2lo (0x1.7f7d1cp-20f) |
| |
| #if SV_EXPF_USE_FEXPA |
| |
| #define Shift (0x1.903f8p17f) /* 1.5*2^17 + 127. */ |
| #define Thres \ |
| (0x1.5d5e2ap+6f) /* Roughly 87.3. For x < -Thres, the result is subnormal \ |
| and not handled correctly by FEXPA. */ |
| |
| static NOINLINE sv_f32_t |
| special_case (sv_f32_t x, sv_f32_t y, svbool_t special) |
| { |
| /* The special-case handler from the Neon routine does not handle subnormals |
| in a way that is compatible with FEXPA. For the FEXPA variant we just fall |
| back to scalar expf. */ |
| return sv_call_f32 (expf, x, y, special); |
| } |
| |
| #else |
| |
| #define Shift (0x1.8p23f) /* 1.5 * 2^23. */ |
| #define Thres (126.0f) |
| |
| /* Special-case handler adapted from Neon variant. Uses s, y and n to produce |
| the final result (normal cases included). It performs an update of all lanes! |
| Therefore: |
| - all previous computation need to be done on all lanes indicated by input |
| pg |
| - we cannot simply apply the special case to the special-case-activated |
| lanes. Besides it is likely that this would not increase performance (no |
| scatter/gather). */ |
| static inline sv_f32_t |
| specialcase (svbool_t pg, sv_f32_t poly, sv_f32_t n, sv_u32_t e, |
| svbool_t p_cmp1, sv_f32_t scale) |
| { |
| /* s=2^(n/N) may overflow, break it up into s=s1*s2, |
| such that exp = s + s*y can be computed as s1*(s2+s2*y) |
| and s1*s1 overflows only if n>0. */ |
| |
| /* If n<=0 then set b to 0x820...0, 0 otherwise. */ |
| svbool_t p_sign = svcmple_n_f32 (pg, n, 0.0f); /* n <= 0. */ |
| sv_u32_t b |
| = svdup_n_u32_z (p_sign, 0x82000000); /* Inactive lanes set to 0. */ |
| |
| /* Set s1 to generate overflow depending on sign of exponent n. */ |
| sv_f32_t s1 |
| = sv_as_f32_u32 (svadd_n_u32_x (pg, b, 0x7f000000)); /* b + 0x7f000000. */ |
| /* Offset s to avoid overflow in final result if n is below threshold. */ |
| sv_f32_t s2 = sv_as_f32_u32 ( |
| svsub_u32_x (pg, e, b)); /* as_u32 (s) - 0x3010...0 + b. */ |
| |
| /* |n| > 192 => 2^(n/N) overflows. */ |
| svbool_t p_cmp2 = svacgt_n_f32 (pg, n, 192.0f); |
| |
| sv_f32_t r2 = svmul_f32_x (pg, s1, s1); |
| sv_f32_t r1 = sv_fma_f32_x (pg, poly, s2, s2); |
| r1 = svmul_f32_x (pg, r1, s1); |
| sv_f32_t r0 = sv_fma_f32_x (pg, poly, scale, scale); |
| |
| /* Apply condition 1 then 2. |
| Returns r2 if cond2 is true, otherwise |
| if cond1 is true then return r1, otherwise return r0. */ |
| sv_f32_t r = svsel_f32 (p_cmp1, r1, r0); |
| |
| return svsel_f32 (p_cmp2, r2, r); |
| } |
| |
| #endif |
| |
| /* Optimised single-precision SVE exp function. By default this is an SVE port |
| of the Neon algorithm from math/. Alternatively, enable a modification of |
| that algorithm that looks up scale using SVE FEXPA instruction with |
| SV_EXPF_USE_FEXPA. |
| |
| Worst-case error of the default algorithm is 1.95 ulp: |
| __sv_expf(-0x1.4cb74ap+2) got 0x1.6a022cp-8 |
| want 0x1.6a023p-8. |
| |
| Worst-case error when using FEXPA is 1.04 ulp: |
| __sv_expf(0x1.a8eda4p+1) got 0x1.ba74bcp+4 |
| want 0x1.ba74bap+4. */ |
| sv_f32_t |
| __sv_expf_x (sv_f32_t x, const svbool_t pg) |
| { |
| /* exp(x) = 2^n (1 + poly(r)), with 1 + poly(r) in [1/sqrt(2),sqrt(2)] |
| x = ln2*n + r, with r in [-ln2/2, ln2/2]. */ |
| |
| /* n = round(x/(ln2/N)). */ |
| sv_f32_t z = sv_fma_n_f32_x (pg, InvLn2, x, sv_f32 (Shift)); |
| sv_f32_t n = svsub_n_f32_x (pg, z, Shift); |
| |
| /* r = x - n*ln2/N. */ |
| sv_f32_t r = sv_fma_n_f32_x (pg, -Ln2hi, n, x); |
| r = sv_fma_n_f32_x (pg, -Ln2lo, n, r); |
| |
| /* scale = 2^(n/N). */ |
| #if SV_EXPF_USE_FEXPA |
| /* NaNs also need special handling with FEXPA. */ |
| svbool_t is_special_case |
| = svorr_b_z (pg, svacgt_n_f32 (pg, x, Thres), svcmpne_f32 (pg, x, x)); |
| sv_f32_t scale = svexpa_f32 (sv_as_u32_f32 (z)); |
| #else |
| sv_u32_t e = svlsl_n_u32_x (pg, sv_as_u32_f32 (z), 23); |
| svbool_t is_special_case = svacgt_n_f32 (pg, n, Thres); |
| sv_f32_t scale = sv_as_f32_u32 (svadd_n_u32_x (pg, e, 0x3f800000)); |
| #endif |
| |
| /* y = exp(r) - 1 ~= r + C1 r^2 + C2 r^3 + C3 r^4. */ |
| sv_f32_t r2 = svmul_f32_x (pg, r, r); |
| sv_f32_t p = sv_fma_n_f32_x (pg, C (0), r, sv_f32 (C (1))); |
| sv_f32_t q = sv_fma_n_f32_x (pg, C (2), r, sv_f32 (C (3))); |
| q = sv_fma_f32_x (pg, p, r2, q); |
| p = svmul_n_f32_x (pg, r, C (4)); |
| sv_f32_t poly = sv_fma_f32_x (pg, q, r2, p); |
| |
| if (unlikely (svptest_any (pg, is_special_case))) |
| #if SV_EXPF_USE_FEXPA |
| return special_case (x, sv_fma_f32_x (pg, poly, scale, scale), |
| is_special_case); |
| #else |
| return specialcase (pg, poly, n, e, is_special_case, scale); |
| #endif |
| |
| return sv_fma_f32_x (pg, poly, scale, scale); |
| } |
| |
| PL_ALIAS (__sv_expf_x, _ZGVsMxv_expf) |
| |
| PL_SIG (SV, F, 1, exp, -9.9, 9.9) |
| PL_TEST_ULP (__sv_expf, 1.46) |
| PL_TEST_INTERVAL (__sv_expf, 0, 0x1p-23, 40000) |
| PL_TEST_INTERVAL (__sv_expf, 0x1p-23, 1, 50000) |
| PL_TEST_INTERVAL (__sv_expf, 1, 0x1p23, 50000) |
| PL_TEST_INTERVAL (__sv_expf, 0x1p23, inf, 50000) |
| PL_TEST_INTERVAL (__sv_expf, -0, -0x1p-23, 40000) |
| PL_TEST_INTERVAL (__sv_expf, -0x1p-23, -1, 50000) |
| PL_TEST_INTERVAL (__sv_expf, -1, -0x1p23, 50000) |
| PL_TEST_INTERVAL (__sv_expf, -0x1p23, -inf, 50000) |
| #endif // SV_SUPPORTED |