1 /*
2 * Double-precision vector cbrt(x) function.
3 *
4 * Copyright (c) 2022-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "v_math.h"
9 #include "test_sig.h"
10 #include "test_defs.h"
11 #include "v_poly_f64.h"
12
13 const static struct data
14 {
15 float64x2_t poly[4], one_third, shift;
16 int64x2_t exp_bias;
17 uint64x2_t abs_mask, tiny_bound;
18 uint32x4_t thresh;
19 double table[5];
20 } data = {
21 .shift = V2 (0x1.8p52),
22 .poly = { /* Generated with fpminimax in [0.5, 1]. */
23 V2 (0x1.c14e8ee44767p-2), V2 (0x1.dd2d3f99e4c0ep-1),
24 V2 (-0x1.08e83026b7e74p-1), V2 (0x1.2c74eaa3ba428p-3) },
25 .exp_bias = V2 (1022),
26 .abs_mask = V2(0x7fffffffffffffff),
27 .tiny_bound = V2(0x0010000000000000), /* Smallest normal. */
28 .thresh = V4(0x7fe00000), /* asuint64 (infinity) - tiny_bound. */
29 .one_third = V2(0x1.5555555555555p-2),
30 .table = { /* table[i] = 2^((i - 2) / 3). */
31 0x1.428a2f98d728bp-1, 0x1.965fea53d6e3dp-1, 0x1p0,
32 0x1.428a2f98d728bp0, 0x1.965fea53d6e3dp0 }
33 };
34
35 #define MantissaMask v_u64 (0x000fffffffffffff)
36
37 static float64x2_t NOINLINE VPCS_ATTR
special_case(float64x2_t x,float64x2_t y,uint32x2_t special)38 special_case (float64x2_t x, float64x2_t y, uint32x2_t special)
39 {
40 return v_call_f64 (cbrt, x, y, vmovl_u32 (special));
41 }
42
43 /* Approximation for double-precision vector cbrt(x), using low-order
44 polynomial and two Newton iterations.
45
46 The vector version of frexp does not handle subnormals
47 correctly. As a result these need to be handled by the scalar
48 fallback, where accuracy may be worse than that of the vector code
49 path.
50
51 Greatest observed error in the normal range is 1.79 ULP. Errors repeat
52 according to the exponent, for instance an error observed for double value
53 m * 2^e will be observed for any input m * 2^(e + 3*i), where i is an
54 integer.
55 _ZGVnN2v_cbrt (0x1.fffff403f0bc6p+1) got 0x1.965fe72821e9bp+0
56 want 0x1.965fe72821e99p+0. */
V_NAME_D1(cbrt)57 VPCS_ATTR float64x2_t V_NAME_D1 (cbrt) (float64x2_t x)
58 {
59 const struct data *d = ptr_barrier (&data);
60 uint64x2_t iax = vreinterpretq_u64_f64 (vabsq_f64 (x));
61
62 /* Subnormal, +/-0 and special values. */
63 uint32x2_t special
64 = vcge_u32 (vsubhn_u64 (iax, d->tiny_bound), vget_low_u32 (d->thresh));
65
66 /* Decompose |x| into m * 2^e, where m is in [0.5, 1.0]. This is a vector
67 version of frexp, which gets subnormal values wrong - these have to be
68 special-cased as a result. */
69 float64x2_t m = vbslq_f64 (MantissaMask, x, v_f64 (0.5));
70 int64x2_t exp_bias = d->exp_bias;
71 uint64x2_t ia12 = vshrq_n_u64 (iax, 52);
72 int64x2_t e = vsubq_s64 (vreinterpretq_s64_u64 (ia12), exp_bias);
73
74 /* Calculate rough approximation for cbrt(m) in [0.5, 1.0], starting point
75 for Newton iterations. */
76 float64x2_t p = v_pairwise_poly_3_f64 (m, vmulq_f64 (m, m), d->poly);
77 float64x2_t one_third = d->one_third;
78 /* Two iterations of Newton's method for iteratively approximating cbrt. */
79 float64x2_t m_by_3 = vmulq_f64 (m, one_third);
80 float64x2_t two_thirds = vaddq_f64 (one_third, one_third);
81 float64x2_t a
82 = vfmaq_f64 (vdivq_f64 (m_by_3, vmulq_f64 (p, p)), two_thirds, p);
83 a = vfmaq_f64 (vdivq_f64 (m_by_3, vmulq_f64 (a, a)), two_thirds, a);
84
85 /* Assemble the result by the following:
86
87 cbrt(x) = cbrt(m) * 2 ^ (e / 3).
88
89 We can get 2 ^ round(e / 3) using ldexp and integer divide, but since e is
90 not necessarily a multiple of 3 we lose some information.
91
92 Let q = 2 ^ round(e / 3), then t = 2 ^ (e / 3) / q.
93
94 Then we know t = 2 ^ (i / 3), where i is the remainder from e / 3, which
95 is an integer in [-2, 2], and can be looked up in the table T. Hence the
96 result is assembled as:
97
98 cbrt(x) = cbrt(m) * t * 2 ^ round(e / 3) * sign. */
99
100 float64x2_t ef = vcvtq_f64_s64 (e);
101 float64x2_t eb3f = vrndnq_f64 (vmulq_f64 (ef, one_third));
102 int64x2_t em3 = vcvtq_s64_f64 (vfmsq_f64 (ef, eb3f, v_f64 (3)));
103 int64x2_t ey = vcvtq_s64_f64 (eb3f);
104
105 float64x2_t my = (float64x2_t){ d->table[em3[0] + 2], d->table[em3[1] + 2] };
106 my = vmulq_f64 (my, a);
107
108 /* Vector version of ldexp. */
109 float64x2_t y = vreinterpretq_f64_s64 (
110 vshlq_n_s64 (vaddq_s64 (ey, vaddq_s64 (exp_bias, v_s64 (1))), 52));
111 y = vmulq_f64 (y, my);
112
113 if (unlikely (v_any_u32h (special)))
114 return special_case (x, vbslq_f64 (d->abs_mask, y, x), special);
115
116 /* Copy sign. */
117 return vbslq_f64 (d->abs_mask, y, x);
118 }
119
120 /* Worse-case ULP error assumes that scalar fallback is GLIBC 2.40 cbrt, which
121 has ULP error of 3.67 at 0x1.7a337e1ba1ec2p-257 [1]. Largest observed error
122 in the vector path is 1.79 ULP.
123 [1] Innocente, V., & Zimmermann, P. (2024). Accuracy of Mathematical
124 Functions in Single, Double, Double Extended, and Quadruple Precision. */
125 TEST_ULP (V_NAME_D1 (cbrt), 3.17)
126 TEST_SIG (V, D, 1, cbrt, -10.0, 10.0)
127 TEST_SYM_INTERVAL (V_NAME_D1 (cbrt), 0, inf, 1000000)
128