1 /* 2 * Single-precision SVE cbrt(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 "pl_sig.h" 10 #include "pl_test.h" 11 #include "poly_sve_f32.h" 12 13 const static struct data 14 { 15 float32_t poly[4]; 16 float32_t table[5]; 17 float32_t one_third, two_thirds; 18 } data = { 19 /* Very rough approximation of cbrt(x) in [0.5, 1], generated with FPMinimax. 20 */ 21 .poly = { 0x1.c14e96p-2, 0x1.dd2d3p-1, -0x1.08e81ap-1, 22 0x1.2c74c2p-3, }, 23 /* table[i] = 2^((i - 2) / 3). */ 24 .table = { 0x1.428a3p-1, 0x1.965feap-1, 0x1p0, 0x1.428a3p0, 0x1.965feap0 }, 25 .one_third = 0x1.555556p-2f, 26 .two_thirds = 0x1.555556p-1f, 27 }; 28 29 #define SmallestNormal 0x00800000 30 #define Thresh 0x7f000000 /* asuint(INFINITY) - SmallestNormal. */ 31 #define MantissaMask 0x007fffff 32 #define HalfExp 0x3f000000 33 34 static svfloat32_t NOINLINE 35 special_case (svfloat32_t x, svfloat32_t y, svbool_t special) 36 { 37 return sv_call_f32 (cbrtf, x, y, special); 38 } 39 40 static inline svfloat32_t 41 shifted_lookup (const svbool_t pg, const float32_t *table, svint32_t i) 42 { 43 return svld1_gather_index (pg, table, svadd_x (pg, i, 2)); 44 } 45 46 /* Approximation for vector single-precision cbrt(x) using Newton iteration 47 with initial guess obtained by a low-order polynomial. Greatest error 48 is 1.64 ULP. This is observed for every value where the mantissa is 49 0x1.85a2aa and the exponent is a multiple of 3, for example: 50 _ZGVsMxv_cbrtf (0x1.85a2aap+3) got 0x1.267936p+1 51 want 0x1.267932p+1. */ 52 svfloat32_t SV_NAME_F1 (cbrt) (svfloat32_t x, const svbool_t pg) 53 { 54 const struct data *d = ptr_barrier (&data); 55 56 svfloat32_t ax = svabs_x (pg, x); 57 svuint32_t iax = svreinterpret_u32 (ax); 58 svuint32_t sign = sveor_x (pg, svreinterpret_u32 (x), iax); 59 60 /* Subnormal, +/-0 and special values. */ 61 svbool_t special = svcmpge (pg, svsub_x (pg, iax, SmallestNormal), Thresh); 62 63 /* Decompose |x| into m * 2^e, where m is in [0.5, 1.0]. This is a vector 64 version of frexpf, which gets subnormal values wrong - these have to be 65 special-cased as a result. */ 66 svfloat32_t m = svreinterpret_f32 (svorr_x ( 67 pg, svand_x (pg, svreinterpret_u32 (x), MantissaMask), HalfExp)); 68 svint32_t e = svsub_x (pg, svreinterpret_s32 (svlsr_x (pg, iax, 23)), 126); 69 70 /* p is a rough approximation for cbrt(m) in [0.5, 1.0]. The better this is, 71 the less accurate the next stage of the algorithm needs to be. An order-4 72 polynomial is enough for one Newton iteration. */ 73 svfloat32_t p 74 = sv_pairwise_poly_3_f32_x (pg, m, svmul_x (pg, m, m), d->poly); 75 76 /* One iteration of Newton's method for iteratively approximating cbrt. */ 77 svfloat32_t m_by_3 = svmul_x (pg, m, d->one_third); 78 svfloat32_t a = svmla_x (pg, svdiv_x (pg, m_by_3, svmul_x (pg, p, p)), p, 79 d->two_thirds); 80 81 /* Assemble the result by the following: 82 83 cbrt(x) = cbrt(m) * 2 ^ (e / 3). 84 85 We can get 2 ^ round(e / 3) using ldexp and integer divide, but since e is 86 not necessarily a multiple of 3 we lose some information. 87 88 Let q = 2 ^ round(e / 3), then t = 2 ^ (e / 3) / q. 89 90 Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3, which 91 is an integer in [-2, 2], and can be looked up in the table T. Hence the 92 result is assembled as: 93 94 cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign. */ 95 svfloat32_t ef = svmul_x (pg, svcvt_f32_x (pg, e), d->one_third); 96 svint32_t ey = svcvt_s32_x (pg, ef); 97 svint32_t em3 = svmls_x (pg, e, ey, 3); 98 99 svfloat32_t my = shifted_lookup (pg, d->table, em3); 100 my = svmul_x (pg, my, a); 101 102 /* Vector version of ldexpf. */ 103 svfloat32_t y = svscale_x (pg, my, ey); 104 105 if (unlikely (svptest_any (pg, special))) 106 return special_case ( 107 x, svreinterpret_f32 (svorr_x (pg, svreinterpret_u32 (y), sign)), 108 special); 109 110 /* Copy sign. */ 111 return svreinterpret_f32 (svorr_x (pg, svreinterpret_u32 (y), sign)); 112 } 113 114 PL_SIG (SV, F, 1, cbrt, -10.0, 10.0) 115 PL_TEST_ULP (SV_NAME_F1 (cbrt), 1.15) 116 PL_TEST_SYM_INTERVAL (SV_NAME_F1 (cbrt), 0, inf, 1000000) 117