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