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;
30
31#[inline(always)]
32pub(crate) fn halley_refine_d(x: f64, a: f64) -> f64 {
33    let tx = x * x * x;
34    x * f_fmla(2., a, tx) / f_fmla(2., tx, a)
35}
36
37#[inline(always)]
38const fn halley_refine(x: f32, a: f32) -> f32 {
39    let tx = x * x * x;
40    x * (tx + 2f32 * a) / (2f32 * tx + a)
41}
42
43/// Cbrt for given value for const context.
44/// This is simplified version just to make a good approximation on const context.
45#[inline]
46pub const fn cbrtf(x: f32) -> f32 {
47    let u = x.to_bits();
48    let au = u.wrapping_shl(1);
49    if au < (1u32 << 24) || au >= (0xffu32 << 24) {
50        if au >= (0xffu32 << 24) {
51            return x + x; /* inf, nan */
52        }
53        if au == 0 {
54            return x; /* +-0 */
55        }
56    }
57
58    const B1: u32 = 709958130;
59    let mut t: f32;
60    let mut ui: u32 = x.to_bits();
61    let mut hx: u32 = ui & 0x7fffffff;
62
63    hx = (hx / 3).wrapping_add(B1);
64    ui &= 0x80000000;
65    ui |= hx;
66
67    t = f32::from_bits(ui);
68    t = halley_refine(t, x);
69    halley_refine(t, x)
70}
71
72/// Computes cube root
73///
74/// Peak ULP on 64 bit = 0.49999577
75#[inline]
76pub fn f_cbrtf(x: f32) -> f32 {
77    let u = x.to_bits();
78    let au = u.wrapping_shl(1);
79    if au < (1u32 << 24) || au >= (0xffu32 << 24) {
80        if au >= (0xffu32 << 24) {
81            return x + x; /* inf, nan */
82        }
83        if au == 0 {
84            return x; /* +-0 */
85        }
86    }
87
88    let mut ui: u32 = x.to_bits();
89    let mut hx: u32 = ui & 0x7fffffff;
90
91    if hx < 0x00800000 {
92        /* zero or subnormal? */
93        if hx == 0 {
94            return x; /* cbrt(+-0) is itself */
95        }
96        const TWO_EXP_24: f32 = f32::from_bits(0x4b800000);
97        ui = (x * TWO_EXP_24).to_bits();
98        hx = ui & 0x7fffffff;
99        const B2: u32 = 642849266;
100        hx = (hx / 3).wrapping_add(B2);
101    } else {
102        const B1: u32 = 709958130;
103        hx = (hx / 3).wrapping_add(B1);
104    }
105    ui &= 0x80000000;
106    ui |= hx;
107
108    let mut t = f32::from_bits(ui) as f64;
109    let dx = x as f64;
110    t = halley_refine_d(t, dx);
111    halley_refine_d(t, dx) as f32
112}
113
114#[cfg(test)]
115mod tests {
116    use super::*;
117
118    #[test]
119    fn test_fcbrtf() {
120        assert_eq!(f_cbrtf(0.0), 0.0);
121        assert_eq!(f_cbrtf(-27.0), -3.0);
122        assert_eq!(f_cbrtf(27.0), 3.0);
123        assert_eq!(f_cbrtf(64.0), 4.0);
124        assert_eq!(f_cbrtf(-64.0), -4.0);
125        assert_eq!(f_cbrtf(f32::NEG_INFINITY), f32::NEG_INFINITY);
126        assert_eq!(f_cbrtf(f32::INFINITY), f32::INFINITY);
127        assert!(f_cbrtf(f32::NAN).is_nan());
128    }
129}
pxfm/cube_roots/cbrtf.rs

pxfm/cube_roots/
cbrtf.rs
