1 /* 2 * Single-precision SVE 2^x function. 3 * 4 * Copyright (c) 2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "sv_math.h" 9 #include "include/mathlib.h" 10 #include "pl_sig.h" 11 #include "pl_test.h" 12 #include "poly_sve_f32.h" 13 14 /* For x < -SpecialBound, the result is subnormal and not handled correctly by 15 FEXPA. */ 16 #define SpecialBound 37.9 17 18 static const struct data 19 { 20 float poly[5]; 21 float shift, log10_2, log2_10_hi, log2_10_lo, special_bound; 22 } data = { 23 /* Coefficients generated using Remez algorithm with minimisation of relative 24 error. 25 rel error: 0x1.89dafa3p-24 26 abs error: 0x1.167d55p-23 in [-log10(2)/2, log10(2)/2] 27 maxerr: 0.52 +0.5 ulp. */ 28 .poly = { 0x1.26bb16p+1f, 0x1.5350d2p+1f, 0x1.04744ap+1f, 0x1.2d8176p+0f, 29 0x1.12b41ap-1f }, 30 /* 1.5*2^17 + 127, a shift value suitable for FEXPA. */ 31 .shift = 0x1.903f8p17f, 32 .log10_2 = 0x1.a934fp+1, 33 .log2_10_hi = 0x1.344136p-2, 34 .log2_10_lo = -0x1.ec10cp-27, 35 .special_bound = SpecialBound, 36 }; 37 38 static svfloat32_t NOINLINE 39 special_case (svfloat32_t x, svfloat32_t y, svbool_t special) 40 { 41 return sv_call_f32 (exp10f, x, y, special); 42 } 43 44 /* Single-precision SVE exp10f routine. Implements the same algorithm 45 as AdvSIMD exp10f. 46 Worst case error is 1.02 ULPs. 47 _ZGVsMxv_exp10f(-0x1.040488p-4) got 0x1.ba5f9ep-1 48 want 0x1.ba5f9cp-1. */ 49 svfloat32_t SV_NAME_F1 (exp10) (svfloat32_t x, const svbool_t pg) 50 { 51 const struct data *d = ptr_barrier (&data); 52 /* exp10(x) = 2^(n/N) * 10^r = 2^n * (1 + poly (r)), 53 with poly(r) in [1/sqrt(2), sqrt(2)] and 54 x = r + n * log10(2) / N, with r in [-log10(2)/2N, log10(2)/2N]. */ 55 56 /* Load some constants in quad-word chunks to minimise memory access (last 57 lane is wasted). */ 58 svfloat32_t log10_2_and_inv = svld1rq (svptrue_b32 (), &d->log10_2); 59 60 /* n = round(x/(log10(2)/N)). */ 61 svfloat32_t shift = sv_f32 (d->shift); 62 svfloat32_t z = svmla_lane (shift, x, log10_2_and_inv, 0); 63 svfloat32_t n = svsub_x (pg, z, shift); 64 65 /* r = x - n*log10(2)/N. */ 66 svfloat32_t r = svmls_lane (x, n, log10_2_and_inv, 1); 67 r = svmls_lane (r, n, log10_2_and_inv, 2); 68 69 svbool_t special = svacgt (pg, x, d->special_bound); 70 svfloat32_t scale = svexpa (svreinterpret_u32 (z)); 71 72 /* Polynomial evaluation: poly(r) ~ exp10(r)-1. */ 73 svfloat32_t r2 = svmul_x (pg, r, r); 74 svfloat32_t poly 75 = svmla_x (pg, svmul_x (pg, r, d->poly[0]), 76 sv_pairwise_poly_3_f32_x (pg, r, r2, d->poly + 1), r2); 77 78 if (unlikely (svptest_any (pg, special))) 79 return special_case (x, svmla_x (pg, scale, scale, poly), special); 80 81 return svmla_x (pg, scale, scale, poly); 82 } 83 84 PL_SIG (SV, F, 1, exp10, -9.9, 9.9) 85 PL_TEST_ULP (SV_NAME_F1 (exp10), 0.52) 86 PL_TEST_SYM_INTERVAL (SV_NAME_F1 (exp10), 0, SpecialBound, 50000) 87 PL_TEST_SYM_INTERVAL (SV_NAME_F1 (exp10), SpecialBound, inf, 50000) 88