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