1 /* 2 * Single-precision 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 "poly_scalar_f32.h" 9 #include "math_config.h" 10 #include "pl_sig.h" 11 #include "pl_test.h" 12 13 #define AbsMask 0x7fffffff 14 #define SignMask 0x80000000 15 #define TwoThirds 0x1.555556p-1f 16 17 #define T(i) __cbrtf_data.table[i] 18 19 /* Approximation for single-precision cbrt(x), using low-order polynomial and 20 one Newton iteration on a reduced interval. Greatest error is 1.5 ULP. This 21 is observed for every value where the mantissa is 0x1.81410e and the exponent 22 is a multiple of 3, for example: 23 cbrtf(0x1.81410ep+30) got 0x1.255d96p+10 24 want 0x1.255d92p+10. */ 25 float 26 cbrtf (float x) 27 { 28 uint32_t ix = asuint (x); 29 uint32_t iax = ix & AbsMask; 30 uint32_t sign = ix & SignMask; 31 32 if (unlikely (iax == 0 || iax == 0x7f800000)) 33 return x; 34 35 /* |x| = m * 2^e, where m is in [0.5, 1.0]. 36 We can easily decompose x into m and e using frexpf. */ 37 int e; 38 float m = frexpf (asfloat (iax), &e); 39 40 /* p is a rough approximation for cbrt(m) in [0.5, 1.0]. The better this is, 41 the less accurate the next stage of the algorithm needs to be. An order-4 42 polynomial is enough for one Newton iteration. */ 43 float p = pairwise_poly_3_f32 (m, m * m, __cbrtf_data.poly); 44 45 /* One iteration of Newton's method for iteratively approximating cbrt. */ 46 float m_by_3 = m / 3; 47 float a = fmaf (TwoThirds, p, m_by_3 / (p * p)); 48 49 /* Assemble the result by the following: 50 51 cbrt(x) = cbrt(m) * 2 ^ (e / 3). 52 53 Let t = (2 ^ (e / 3)) / (2 ^ round(e / 3)). 54 55 Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3. 56 i is an integer in [-2, 2], so t can be looked up in the table T. 57 Hence the result is assembled as: 58 59 cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign. 60 Which can be done easily using ldexpf. */ 61 return asfloat (asuint (ldexpf (a * T (2 + e % 3), e / 3)) | sign); 62 } 63 64 PL_SIG (S, F, 1, cbrt, -10.0, 10.0) 65 PL_TEST_ULP (cbrtf, 1.03) 66 PL_TEST_SYM_INTERVAL (cbrtf, 0, inf, 1000000) 67