1 /* 2 * Double-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 "math_config.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 PL_SIG (S, D, 1, cbrt, -10.0, 10.0) 13 14 #define AbsMask 0x7fffffffffffffff 15 #define TwoThirds 0x1.5555555555555p-1 16 17 #define C(i) __cbrt_data.poly[i] 18 #define T(i) __cbrt_data.table[i] 19 20 /* Approximation for double-precision cbrt(x), using low-order polynomial and 21 two Newton iterations. Greatest observed error is 1.79 ULP. Errors repeat 22 according to the exponent, for instance an error observed for double value 23 m * 2^e will be observed for any input m * 2^(e + 3*i), where i is an 24 integer. 25 cbrt(0x1.fffff403f0bc6p+1) got 0x1.965fe72821e9bp+0 26 want 0x1.965fe72821e99p+0. */ 27 double 28 cbrt (double x) 29 { 30 uint64_t ix = asuint64 (x); 31 uint64_t iax = ix & AbsMask; 32 uint64_t sign = ix & ~AbsMask; 33 34 if (unlikely (iax == 0 || iax == 0x7ff0000000000000)) 35 return x; 36 37 /* |x| = m * 2^e, where m is in [0.5, 1.0]. 38 We can easily decompose x into m and e using frexp. */ 39 int e; 40 double m = frexp (asdouble (iax), &e); 41 42 /* Calculate rough approximation for cbrt(m) in [0.5, 1.0], starting point for 43 Newton iterations. */ 44 double p_01 = fma (C (1), m, C (0)); 45 double p_23 = fma (C (3), m, C (2)); 46 double p = fma (p_23, m * m, p_01); 47 48 /* Two iterations of Newton's method for iteratively approximating cbrt. */ 49 double m_by_3 = m / 3; 50 double a = fma (TwoThirds, p, m_by_3 / (p * p)); 51 a = fma (TwoThirds, a, m_by_3 / (a * a)); 52 53 /* Assemble the result by the following: 54 55 cbrt(x) = cbrt(m) * 2 ^ (e / 3). 56 57 Let t = (2 ^ (e / 3)) / (2 ^ round(e / 3)). 58 59 Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3. 60 i is an integer in [-2, 2], so t can be looked up in the table T. 61 Hence the result is assembled as: 62 63 cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign. 64 Which can be done easily using ldexp. */ 65 return asdouble (asuint64 (ldexp (a * T (2 + e % 3), e / 3)) | sign); 66 } 67 68 PL_TEST_ULP (cbrt, 1.30) 69 PL_TEST_SYM_INTERVAL (cbrt, 0, inf, 1000000) 70