pxfm/cube_roots/
cbrtf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 4/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::f_fmla;
30use crate::cube_roots::cbrt::{CbrtBackend, GenericCbrtBackend};
31
32#[inline(always)]
33pub(crate) fn halley_refine_d(x: f64, a: f64) -> f64 {
34    let tx = x * x * x;
35    x * f_fmla(2., a, tx) / f_fmla(2., tx, a)
36}
37
38#[inline(always)]
39#[allow(unused)]
40pub(crate) fn halley_refine_d_fma(x: f64, a: f64) -> f64 {
41    let tx = x * x * x;
42    x * f64::mul_add(2., a, tx) / f64::mul_add(2., tx, a)
43}
44
45#[inline(always)]
46const fn halley_refine(x: f32, a: f32) -> f32 {
47    let tx = x * x * x;
48    x * (tx + 2f32 * a) / (2f32 * tx + a)
49}
50
51/// Cbrt for given value for const context.
52/// This is simplified version just to make a good approximation on const context.
53#[inline]
54pub const fn cbrtf(x: f32) -> f32 {
55    let u = x.to_bits();
56    let au = u.wrapping_shl(1);
57    if au < (1u32 << 24) || au >= (0xffu32 << 24) {
58        if au >= (0xffu32 << 24) {
59            return x + x; /* inf, nan */
60        }
61        if au == 0 {
62            return x; /* +-0 */
63        }
64    }
65
66    const B1: u32 = 709958130;
67    let mut t: f32;
68    let mut ui: u32 = x.to_bits();
69    let mut hx: u32 = ui & 0x7fffffff;
70
71    hx = (hx / 3).wrapping_add(B1);
72    ui &= 0x80000000;
73    ui |= hx;
74
75    t = f32::from_bits(ui);
76    t = halley_refine(t, x);
77    halley_refine(t, x)
78}
79
80#[inline(always)]
81fn cbrtf_gen_impl<B: CbrtBackend>(x: f32, backend: B) -> f32 {
82    let u = x.to_bits();
83    let au = u.wrapping_shl(1);
84    if au < (1u32 << 24) || au >= (0xffu32 << 24) {
85        if au >= (0xffu32 << 24) {
86            return x + x; /* inf, nan */
87        }
88        if au == 0 {
89            return x; /* +-0 */
90        }
91    }
92
93    let mut ui: u32 = x.to_bits();
94    let mut hx: u32 = ui & 0x7fffffff;
95
96    if hx < 0x00800000 {
97        /* zero or subnormal? */
98        if hx == 0 {
99            return x; /* cbrt(+-0) is itself */
100        }
101        const TWO_EXP_24: f32 = f32::from_bits(0x4b800000);
102        ui = (x * TWO_EXP_24).to_bits();
103        hx = ui & 0x7fffffff;
104        const B2: u32 = 642849266;
105        hx = (hx / 3).wrapping_add(B2);
106    } else {
107        const B1: u32 = 709958130;
108        hx = (hx / 3).wrapping_add(B1);
109    }
110    ui &= 0x80000000;
111    ui |= hx;
112
113    let mut t = f32::from_bits(ui) as f64;
114    let dx = x as f64;
115    t = backend.halley(t, dx);
116    backend.halley(t, dx) as f32
117}
118
119#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
120#[target_feature(enable = "avx", enable = "fma")]
121unsafe fn cbrtf_fma_impl(x: f32) -> f32 {
122    use crate::cube_roots::cbrt::FmaCbrtBackend;
123    cbrtf_gen_impl(x, FmaCbrtBackend {})
124}
125
126/// Computes cube root
127///
128/// Peak ULP on 64 bit = 0.49999577
129#[inline]
130pub fn f_cbrtf(x: f32) -> f32 {
131    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
132    {
133        cbrtf_gen_impl(x, GenericCbrtBackend {})
134    }
135    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
136    {
137        use std::sync::OnceLock;
138        static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
139        let q = EXECUTOR.get_or_init(|| {
140            if std::arch::is_x86_feature_detected!("avx")
141                && std::arch::is_x86_feature_detected!("fma")
142            {
143                cbrtf_fma_impl
144            } else {
145                fn def_cbrt(x: f32) -> f32 {
146                    cbrtf_gen_impl(x, GenericCbrtBackend {})
147                }
148                def_cbrt
149            }
150        });
151        unsafe { q(x) }
152    }
153}
154
155#[cfg(test)]
156mod tests {
157    use super::*;
158
159    #[test]
160    fn test_fcbrtf() {
161        assert_eq!(f_cbrtf(0.0), 0.0);
162        assert_eq!(f_cbrtf(-27.0), -3.0);
163        assert_eq!(f_cbrtf(27.0), 3.0);
164        assert_eq!(f_cbrtf(64.0), 4.0);
165        assert_eq!(f_cbrtf(-64.0), -4.0);
166        assert_eq!(f_cbrtf(f32::NEG_INFINITY), f32::NEG_INFINITY);
167        assert_eq!(f_cbrtf(f32::INFINITY), f32::INFINITY);
168        assert!(f_cbrtf(f32::NAN).is_nan());
169    }
170}
pxfm/cube_roots/cbrtf.rs

pxfm/cube_roots/
cbrtf.rs
