compiler_builtins/math/libm_math/
cbrt.rs

1/* SPDX-License-Identifier: MIT */
2/* origin: core-math/src/binary64/cbrt/cbrt.c
3 * Copyright (c) 2021-2022 Alexei Sibidanov.
4 * Ported to Rust in 2025 by Trevor Gross.
5 */
6
7use super::Float;
8use super::support::{FpResult, Round, cold_path};
9
10/// Compute the cube root of the argument.
11#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
12pub fn cbrt(x: f64) -> f64 {
13    cbrt_round(x, Round::Nearest).val
14}
15
16pub fn cbrt_round(x: f64, round: Round) -> FpResult<f64> {
17    const ESCALE: [f64; 3] = [
18        1.0,
19        hf64!("0x1.428a2f98d728bp+0"), /* 2^(1/3) */
20        hf64!("0x1.965fea53d6e3dp+0"), /* 2^(2/3) */
21    ];
22
23    /* the polynomial c0+c1*x+c2*x^2+c3*x^3 approximates x^(1/3) on [1,2]
24    with maximal error < 9.2e-5 (attained at x=2) */
25    const C: [f64; 4] = [
26        hf64!("0x1.1b0babccfef9cp-1"),
27        hf64!("0x1.2c9a3e94d1da5p-1"),
28        hf64!("-0x1.4dc30b1a1ddbap-3"),
29        hf64!("0x1.7a8d3e4ec9b07p-6"),
30    ];
31
32    let u0: f64 = hf64!("0x1.5555555555555p-2");
33    let u1: f64 = hf64!("0x1.c71c71c71c71cp-3");
34
35    let rsc = [1.0, -1.0, 0.5, -0.5, 0.25, -0.25];
36
37    let off = [hf64!("0x1p-53"), 0.0, 0.0, 0.0];
38
39    /* rm=0 for rounding to nearest, and other values for directed roundings */
40    let hx: u64 = x.to_bits();
41    let mut mant: u64 = hx & f64::SIG_MASK;
42    let sign: u64 = hx >> 63;
43
44    let mut e: u32 = (hx >> f64::SIG_BITS) as u32 & f64::EXP_SAT;
45
46    if ((e + 1) & f64::EXP_SAT) < 2 {
47        cold_path();
48
49        let ix: u64 = hx & !f64::SIGN_MASK;
50
51        /* 0, inf, nan: we return x + x instead of simply x,
52        to that for x a signaling NaN, it correctly triggers
53        the invalid exception. */
54        if e == f64::EXP_SAT || ix == 0 {
55            return FpResult::ok(x + x);
56        }
57
58        let nz = ix.leading_zeros() - 11; /* subnormal */
59        mant <<= nz;
60        mant &= f64::SIG_MASK;
61        e = e.wrapping_sub(nz - 1);
62    }
63
64    e = e.wrapping_add(3072);
65    let cvt1: u64 = mant | (0x3ffu64 << 52);
66    let mut cvt5: u64 = cvt1;
67
68    let et: u32 = e / 3;
69    let it: u32 = e % 3;
70
71    /* 2^(3k+it) <= x < 2^(3k+it+1), with 0 <= it <= 3 */
72    cvt5 += u64::from(it) << f64::SIG_BITS;
73    cvt5 |= sign << 63;
74    let zz: f64 = f64::from_bits(cvt5);
75
76    /* cbrt(x) = cbrt(zz)*2^(et-1365) where 1 <= zz < 8 */
77    let mut isc: u64 = ESCALE[it as usize].to_bits(); // todo: index
78    isc |= sign << 63;
79    let cvt2: u64 = isc;
80    let z: f64 = f64::from_bits(cvt1);
81
82    /* cbrt(zz) = cbrt(z)*isc, where isc encodes 1, 2^(1/3) or 2^(2/3),
83    and 1 <= z < 2 */
84    let r: f64 = 1.0 / z;
85    let rr: f64 = r * rsc[((it as usize) << 1) | sign as usize];
86    let z2: f64 = z * z;
87    let c0: f64 = C[0] + z * C[1];
88    let c2: f64 = C[2] + z * C[3];
89    let mut y: f64 = c0 + z2 * c2;
90    let mut y2: f64 = y * y;
91
92    /* y is an approximation of z^(1/3) */
93    let mut h: f64 = y2 * (y * r) - 1.0;
94
95    /* h determines the error between y and z^(1/3) */
96    y -= (h * y) * (u0 - u1 * h);
97
98    /* The correction y -= (h*y)*(u0 - u1*h) corresponds to a cubic variant
99    of Newton's method, with the function f(y) = 1-z/y^3. */
100    y *= f64::from_bits(cvt2);
101
102    /* Now y is an approximation of zz^(1/3),
103     * and rr an approximation of 1/zz. We now perform another iteration of
104     * Newton-Raphson, this time with a linear approximation only. */
105    y2 = y * y;
106    let mut y2l: f64 = y.fma(y, -y2);
107
108    /* y2 + y2l = y^2 exactly */
109    let mut y3: f64 = y2 * y;
110    let mut y3l: f64 = y.fma(y2, -y3) + y * y2l;
111
112    /* y3 + y3l approximates y^3 with about 106 bits of accuracy */
113    h = ((y3 - zz) + y3l) * rr;
114    let mut dy: f64 = h * (y * u0);
115
116    /* the approximation of zz^(1/3) is y - dy */
117    let mut y1: f64 = y - dy;
118    dy = (y - y1) - dy;
119
120    /* the approximation of zz^(1/3) is now y1 + dy, where |dy| < 1/2 ulp(y)
121     * (for rounding to nearest) */
122    let mut ady: f64 = dy.abs();
123
124    /* For directed roundings, ady0 is tiny when dy is tiny, or ady0 is near
125     * from ulp(1);
126     * for rounding to nearest, ady0 is tiny when dy is near from 1/2 ulp(1),
127     * or from 3/2 ulp(1). */
128    let mut ady0: f64 = (ady - off[round as usize]).abs();
129    let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();
130
131    if ady0 < hf64!("0x1p-75") || ady1 < hf64!("0x1p-75") {
132        cold_path();
133
134        y2 = y1 * y1;
135        y2l = y1.fma(y1, -y2);
136        y3 = y2 * y1;
137        y3l = y1.fma(y2, -y3) + y1 * y2l;
138        h = ((y3 - zz) + y3l) * rr;
139        dy = h * (y1 * u0);
140        y = y1 - dy;
141        dy = (y1 - y) - dy;
142        y1 = y;
143        ady = dy.abs();
144        ady0 = (ady - off[round as usize]).abs();
145        ady1 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();
146
147        if ady0 < hf64!("0x1p-98") || ady1 < hf64!("0x1p-98") {
148            cold_path();
149            let azz: f64 = zz.abs();
150
151            // ~ 0x1.79d15d0e8d59b80000000000000ffc3dp+0
152            if azz == hf64!("0x1.9b78223aa307cp+1") {
153                y1 = hf64!("0x1.79d15d0e8d59cp+0").copysign(zz);
154            }
155
156            // ~ 0x1.de87aa837820e80000000000001c0f08p+0
157            if azz == hf64!("0x1.a202bfc89ddffp+2") {
158                y1 = hf64!("0x1.de87aa837820fp+0").copysign(zz);
159            }
160
161            if round != Round::Nearest {
162                let wlist = [
163                    (hf64!("0x1.3a9ccd7f022dbp+0"), hf64!("0x1.1236160ba9b93p+0")), // ~ 0x1.1236160ba9b930000000000001e7e8fap+0
164                    (hf64!("0x1.7845d2faac6fep+0"), hf64!("0x1.23115e657e49cp+0")), // ~ 0x1.23115e657e49c0000000000001d7a799p+0
165                    (hf64!("0x1.d1ef81cbbbe71p+0"), hf64!("0x1.388fb44cdcf5ap+0")), // ~ 0x1.388fb44cdcf5a0000000000002202c55p+0
166                    (hf64!("0x1.0a2014f62987cp+1"), hf64!("0x1.46bcbf47dc1e8p+0")), // ~ 0x1.46bcbf47dc1e8000000000000303aa2dp+0
167                    (hf64!("0x1.fe18a044a5501p+1"), hf64!("0x1.95decfec9c904p+0")), // ~ 0x1.95decfec9c9040000000000000159e8ep+0
168                    (hf64!("0x1.a6bb8c803147bp+2"), hf64!("0x1.e05335a6401dep+0")), // ~ 0x1.e05335a6401de00000000000027ca017p+0
169                    (hf64!("0x1.ac8538a031cbdp+2"), hf64!("0x1.e281d87098de8p+0")), // ~ 0x1.e281d87098de80000000000000ee9314p+0
170                ];
171
172                for (a, b) in wlist {
173                    if azz == a {
174                        let tmp = if round as u64 + sign == 2 {
175                            hf64!("0x1p-52")
176                        } else {
177                            0.0
178                        };
179                        y1 = (b + tmp).copysign(zz);
180                    }
181                }
182            }
183        }
184    }
185
186    let mut cvt3: u64 = y1.to_bits();
187    cvt3 = cvt3.wrapping_add(((et.wrapping_sub(342).wrapping_sub(1023)) as u64) << 52);
188    let m0: u64 = cvt3 << 30;
189    let m1 = m0 >> 63;
190
191    if (m0 ^ m1) <= (1u64 << 30) {
192        cold_path();
193
194        let mut cvt4: u64 = y1.to_bits();
195        cvt4 = (cvt4 + (164 << 15)) & 0xffffffffffff0000u64;
196
197        if ((f64::from_bits(cvt4) - y1) - dy).abs() < hf64!("0x1p-60") || (zz).abs() == 1.0 {
198            cvt3 = (cvt3 + (1u64 << 15)) & 0xffffffffffff0000u64;
199        }
200    }
201
202    FpResult::ok(f64::from_bits(cvt3))
203}
204
205#[cfg(test)]
206mod tests {
207    use super::*;
208
209    #[test]
210    fn spot_checks() {
211        if !cfg!(x86_no_sse) {
212            // Exposes a rounding mode problem. Ignored on i586 because of inaccurate FMA.
213            assert_biteq!(
214                cbrt(f64::from_bits(0xf7f792b28f600000)),
215                f64::from_bits(0xd29ce68655d962f3)
216            );
217        }
218    }
219}
compiler_builtins/math/libm_math/cbrt.rs

compiler_builtins/math/libm_math/
cbrt.rs
