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