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