pxfm/triangle/
hypot3f.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/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::EXP_MASK_F32;
30
31pub fn f_hypot3f(x: f32, y: f32, z: f32) -> f32 {
32    let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffffu32);
33    let y_abs = f32::from_bits(y.to_bits() & 0x7fff_ffffu32);
34    let z_abs = f32::from_bits(z.to_bits() & 0x7fff_ffffu32);
35
36    let a_bits = x_abs.to_bits().max(y_abs.to_bits()).max(z_abs.to_bits());
37    let b_bits = x_abs.to_bits().min(y_abs.to_bits()).min(z_abs.to_bits());
38
39    let a_u = a_bits;
40    let b_u = b_bits;
41
42    if a_u >= EXP_MASK_F32 {
43        // x or y is inf or nan
44        if x_abs.is_nan() || y_abs.is_nan() || z_abs.is_nan() {
45            return f32::NAN;
46        }
47        if x_abs.is_infinite() || y_abs.is_infinite() || z_abs.is_infinite() {
48            return f32::INFINITY;
49        }
50        return f32::from_bits(a_bits);
51    }
52
53    if a_u.wrapping_sub(b_u) >= ((23u32 + 2) << 23) {
54        return x_abs + y_abs + z_abs;
55    }
56
57    #[cfg(any(
58        all(
59            any(target_arch = "x86", target_arch = "x86_64"),
60            target_feature = "fma"
61        ),
62        all(target_arch = "aarch64", target_feature = "neon")
63    ))]
64    {
65        let ad = x as f64;
66        let bd = y as f64;
67        let cd = z as f64;
68        use crate::common::f_fmla;
69        // for FMA environment we're using Kahan style summation which is short and reliable.
70        let w = bd * bd; // RN(bc)
71        let e = f_fmla(-bd, bd, w); // RN(w − bc)
72        let f = f_fmla(ad, ad, w); // RN(ad + w)
73        let f0 = f_fmla(cd, cd, f); // RN(cd + f)
74        let r = e + f0; // RN(f + e)
75        let hyp = r.sqrt(); // sqrt(x^2 + y^2)
76        hyp as f32
77    }
78    #[cfg(not(any(
79        all(
80            any(target_arch = "x86", target_arch = "x86_64"),
81            target_feature = "fma"
82        ),
83        all(target_arch = "aarch64", target_feature = "neon")
84    )))]
85    {
86        let ad = x_abs as f64;
87        let bd = y_abs as f64;
88        let cd = z_abs as f64;
89        use crate::double_double::DoubleDouble;
90        let da = DoubleDouble::from_exact_mult(bd, bd);
91        let db = DoubleDouble::from_exact_mult(ad, ad);
92        let dc = DoubleDouble::from_exact_mult(cd, cd);
93        let f = DoubleDouble::add(DoubleDouble::add(da, db), dc);
94        let cath = f.to_f64().sqrt();
95        cath as f32
96    }
97}
98
99#[cfg(test)]
100mod tests {
101
102    use super::*;
103
104    #[test]
105    fn test_hypot3f() {
106        assert_eq!(f_hypot3f(3.0, 4.0, 12.0), 13.0);
107        assert_eq!(f_hypot3f(6.0, 8.0, 24.0), 26.0);
108        assert_eq!(f_hypot3f(5.0, 12.0, 84.0), 85.0);
109        assert_eq!(f_hypot3f(9.0, 12.0, 20.0), 25.0);
110        assert_eq!(f_hypot3f(1e20, 3.0, 4.0), 1e20);
111        assert_eq!(f_hypot3f(1e-20, 1e-20, 1.0), 1.);
112        assert_eq!(
113            f_hypot3f(f32::MIN_POSITIVE, f32::MIN_POSITIVE, 0.0),
114            1.6624e-38
115        );
116        assert_eq!(f_hypot3f(f32::MAX, f32::MAX, 0.), f32::INFINITY);
117        assert_eq!(f_hypot3f(f32::MAX, 0., 0.), 3.4028235e38);
118        assert_eq!(f_hypot3f(f32::INFINITY, 0., 0.), f32::INFINITY);
119        assert_eq!(f_hypot3f(0., f32::INFINITY, 0.), f32::INFINITY);
120        assert_eq!(f_hypot3f(0., 0., f32::INFINITY), f32::INFINITY);
121        assert!(f_hypot3f(f32::NAN, 0., 0.).is_nan());
122        assert!(f_hypot3f(0., f32::NAN, 0.).is_nan());
123        assert!(f_hypot3f(0., 0., f32::NAN).is_nan());
124    }
125}
pxfm/triangle/hypot3f.rs

pxfm/triangle/
hypot3f.rs
