xref: /freebsd/contrib/arm-optimized-routines/math/aarch64/experimental/cbrt_2u.c (revision f3087bef11543b42e0d69b708f367097a4118d24)
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
cbrt(double x)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