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 "estrinf.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 C(i) __cbrtf_data.poly[i] 18 #define T(i) __cbrtf_data.table[i] 19 20 /* Approximation for single-precision cbrt(x), using low-order polynomial and 21 one Newton iteration on a reduced interval. Greatest error is 1.5 ULP. This 22 is observed for every value where the mantissa is 0x1.81410e and the exponent 23 is a multiple of 3, for example: 24 cbrtf(0x1.81410ep+30) got 0x1.255d96p+10 25 want 0x1.255d92p+10. */ 26 float 27 cbrtf (float x) 28 { 29 uint32_t ix = asuint (x); 30 uint32_t iax = ix & AbsMask; 31 uint32_t sign = ix & SignMask; 32 33 if (unlikely (iax == 0 || iax == 0x7f800000)) 34 return x; 35 36 /* |x| = m * 2^e, where m is in [0.5, 1.0]. 37 We can easily decompose x into m and e using frexpf. */ 38 int e; 39 float m = frexpf (asfloat (iax), &e); 40 41 /* p is a rough approximation for cbrt(m) in [0.5, 1.0]. The better this is, 42 the less accurate the next stage of the algorithm needs to be. An order-4 43 polynomial is enough for one Newton iteration. */ 44 float p = ESTRIN_3 (m, m * m, C); 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_INTERVAL (cbrtf, 0, inf, 1000000) 67 PL_TEST_INTERVAL (cbrtf, -0, -inf, 1000000) 68