1 /* 2 * Single-precision vector cbrt(x) function. 3 * 4 * Copyright (c) 2022-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "v_math.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 #include "poly_advsimd_f32.h" 12 13 const static struct data 14 { 15 float32x4_t poly[4], one_third; 16 float table[5]; 17 } data = { 18 .poly = { /* Very rough approximation of cbrt(x) in [0.5, 1], generated with 19 FPMinimax. */ 20 V4 (0x1.c14e96p-2), V4 (0x1.dd2d3p-1), V4 (-0x1.08e81ap-1), 21 V4 (0x1.2c74c2p-3) }, 22 .table = { /* table[i] = 2^((i - 2) / 3). */ 23 0x1.428a3p-1, 0x1.965feap-1, 0x1p0, 0x1.428a3p0, 0x1.965feap0 }, 24 .one_third = V4 (0x1.555556p-2f), 25 }; 26 27 #define SignMask v_u32 (0x80000000) 28 #define SmallestNormal v_u32 (0x00800000) 29 #define Thresh vdup_n_u16 (0x7f00) /* asuint(INFINITY) - SmallestNormal. */ 30 #define MantissaMask v_u32 (0x007fffff) 31 #define HalfExp v_u32 (0x3f000000) 32 33 static float32x4_t VPCS_ATTR NOINLINE 34 special_case (float32x4_t x, float32x4_t y, uint16x4_t special) 35 { 36 return v_call_f32 (cbrtf, x, y, vmovl_u16 (special)); 37 } 38 39 static inline float32x4_t 40 shifted_lookup (const float *table, int32x4_t i) 41 { 42 return (float32x4_t){ table[i[0] + 2], table[i[1] + 2], table[i[2] + 2], 43 table[i[3] + 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 _ZGVnN4v_cbrtf(0x1.85a2aap+3) got 0x1.267936p+1 51 want 0x1.267932p+1. */ 52 VPCS_ATTR float32x4_t V_NAME_F1 (cbrt) (float32x4_t x) 53 { 54 const struct data *d = ptr_barrier (&data); 55 uint32x4_t iax = vreinterpretq_u32_f32 (vabsq_f32 (x)); 56 57 /* Subnormal, +/-0 and special values. */ 58 uint16x4_t special = vcge_u16 (vsubhn_u32 (iax, SmallestNormal), Thresh); 59 60 /* Decompose |x| into m * 2^e, where m is in [0.5, 1.0]. This is a vector 61 version of frexpf, which gets subnormal values wrong - these have to be 62 special-cased as a result. */ 63 float32x4_t m = vbslq_f32 (MantissaMask, x, v_f32 (0.5)); 64 int32x4_t e 65 = vsubq_s32 (vreinterpretq_s32_u32 (vshrq_n_u32 (iax, 23)), v_s32 (126)); 66 67 /* p is a rough approximation for cbrt(m) in [0.5, 1.0]. The better this is, 68 the less accurate the next stage of the algorithm needs to be. An order-4 69 polynomial is enough for one Newton iteration. */ 70 float32x4_t p = v_pairwise_poly_3_f32 (m, vmulq_f32 (m, m), d->poly); 71 72 float32x4_t one_third = d->one_third; 73 float32x4_t two_thirds = vaddq_f32 (one_third, one_third); 74 75 /* One iteration of Newton's method for iteratively approximating cbrt. */ 76 float32x4_t m_by_3 = vmulq_f32 (m, one_third); 77 float32x4_t a 78 = vfmaq_f32 (vdivq_f32 (m_by_3, vmulq_f32 (p, p)), two_thirds, p); 79 80 /* Assemble the result by the following: 81 82 cbrt(x) = cbrt(m) * 2 ^ (e / 3). 83 84 We can get 2 ^ round(e / 3) using ldexp and integer divide, but since e is 85 not necessarily a multiple of 3 we lose some information. 86 87 Let q = 2 ^ round(e / 3), then t = 2 ^ (e / 3) / q. 88 89 Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3, which 90 is an integer in [-2, 2], and can be looked up in the table T. Hence the 91 result is assembled as: 92 93 cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign. */ 94 float32x4_t ef = vmulq_f32 (vcvtq_f32_s32 (e), one_third); 95 int32x4_t ey = vcvtq_s32_f32 (ef); 96 int32x4_t em3 = vsubq_s32 (e, vmulq_s32 (ey, v_s32 (3))); 97 98 float32x4_t my = shifted_lookup (d->table, em3); 99 my = vmulq_f32 (my, a); 100 101 /* Vector version of ldexpf. */ 102 float32x4_t y 103 = vreinterpretq_f32_s32 (vshlq_n_s32 (vaddq_s32 (ey, v_s32 (127)), 23)); 104 y = vmulq_f32 (y, my); 105 106 if (unlikely (v_any_u16h (special))) 107 return special_case (x, vbslq_f32 (SignMask, x, y), special); 108 109 /* Copy sign. */ 110 return vbslq_f32 (SignMask, x, y); 111 } 112 113 PL_SIG (V, F, 1, cbrt, -10.0, 10.0) 114 PL_TEST_ULP (V_NAME_F1 (cbrt), 1.15) 115 PL_TEST_EXPECT_FENV_ALWAYS (V_NAME_F1 (cbrt)) 116 PL_TEST_SYM_INTERVAL (V_NAME_F1 (cbrt), 0, inf, 1000000) 117