1072a4ba8SAndrew Turner /*
2072a4ba8SAndrew Turner * Single-precision cbrt(x) function.
3072a4ba8SAndrew Turner *
4072a4ba8SAndrew Turner * Copyright (c) 2022-2023, Arm Limited.
5072a4ba8SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6072a4ba8SAndrew Turner */
7072a4ba8SAndrew Turner
8*5a02ffc3SAndrew Turner #include "poly_scalar_f32.h"
9072a4ba8SAndrew Turner #include "math_config.h"
10072a4ba8SAndrew Turner #include "pl_sig.h"
11072a4ba8SAndrew Turner #include "pl_test.h"
12072a4ba8SAndrew Turner
13072a4ba8SAndrew Turner #define AbsMask 0x7fffffff
14072a4ba8SAndrew Turner #define SignMask 0x80000000
15072a4ba8SAndrew Turner #define TwoThirds 0x1.555556p-1f
16072a4ba8SAndrew Turner
17072a4ba8SAndrew Turner #define T(i) __cbrtf_data.table[i]
18072a4ba8SAndrew Turner
19072a4ba8SAndrew Turner /* Approximation for single-precision cbrt(x), using low-order polynomial and
20072a4ba8SAndrew Turner one Newton iteration on a reduced interval. Greatest error is 1.5 ULP. This
21072a4ba8SAndrew Turner is observed for every value where the mantissa is 0x1.81410e and the exponent
22072a4ba8SAndrew Turner is a multiple of 3, for example:
23072a4ba8SAndrew Turner cbrtf(0x1.81410ep+30) got 0x1.255d96p+10
24072a4ba8SAndrew Turner want 0x1.255d92p+10. */
25072a4ba8SAndrew Turner float
cbrtf(float x)26072a4ba8SAndrew Turner cbrtf (float x)
27072a4ba8SAndrew Turner {
28072a4ba8SAndrew Turner uint32_t ix = asuint (x);
29072a4ba8SAndrew Turner uint32_t iax = ix & AbsMask;
30072a4ba8SAndrew Turner uint32_t sign = ix & SignMask;
31072a4ba8SAndrew Turner
32072a4ba8SAndrew Turner if (unlikely (iax == 0 || iax == 0x7f800000))
33072a4ba8SAndrew Turner return x;
34072a4ba8SAndrew Turner
35072a4ba8SAndrew Turner /* |x| = m * 2^e, where m is in [0.5, 1.0].
36072a4ba8SAndrew Turner We can easily decompose x into m and e using frexpf. */
37072a4ba8SAndrew Turner int e;
38072a4ba8SAndrew Turner float m = frexpf (asfloat (iax), &e);
39072a4ba8SAndrew Turner
40072a4ba8SAndrew Turner /* p is a rough approximation for cbrt(m) in [0.5, 1.0]. The better this is,
41072a4ba8SAndrew Turner the less accurate the next stage of the algorithm needs to be. An order-4
42072a4ba8SAndrew Turner polynomial is enough for one Newton iteration. */
43*5a02ffc3SAndrew Turner float p = pairwise_poly_3_f32 (m, m * m, __cbrtf_data.poly);
44*5a02ffc3SAndrew Turner
45072a4ba8SAndrew Turner /* One iteration of Newton's method for iteratively approximating cbrt. */
46072a4ba8SAndrew Turner float m_by_3 = m / 3;
47072a4ba8SAndrew Turner float a = fmaf (TwoThirds, p, m_by_3 / (p * p));
48072a4ba8SAndrew Turner
49072a4ba8SAndrew Turner /* Assemble the result by the following:
50072a4ba8SAndrew Turner
51072a4ba8SAndrew Turner cbrt(x) = cbrt(m) * 2 ^ (e / 3).
52072a4ba8SAndrew Turner
53072a4ba8SAndrew Turner Let t = (2 ^ (e / 3)) / (2 ^ round(e / 3)).
54072a4ba8SAndrew Turner
55072a4ba8SAndrew Turner Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3.
56072a4ba8SAndrew Turner i is an integer in [-2, 2], so t can be looked up in the table T.
57072a4ba8SAndrew Turner Hence the result is assembled as:
58072a4ba8SAndrew Turner
59072a4ba8SAndrew Turner cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign.
60072a4ba8SAndrew Turner Which can be done easily using ldexpf. */
61072a4ba8SAndrew Turner return asfloat (asuint (ldexpf (a * T (2 + e % 3), e / 3)) | sign);
62072a4ba8SAndrew Turner }
63072a4ba8SAndrew Turner
64072a4ba8SAndrew Turner PL_SIG (S, F, 1, cbrt, -10.0, 10.0)
65072a4ba8SAndrew Turner PL_TEST_ULP (cbrtf, 1.03)
66*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (cbrtf, 0, inf, 1000000)
67