pxfm/triangle/
cathetusf.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
31/// Computes the missing leg of a right triangle
32///
33/// Given a hypotenuse `x` and a known leg `y`, returns
34/// `sqrt(x^2 - y^2)` = the length of the other leg.
35///
36/// Domain: requires `|x| >= |y|`. Returns NaN if the input
37/// is outside this range.
38pub fn f_cathetusf(x: f32, y: f32) -> f32 {
39    let x_abs = x.abs();
40    let y_abs = y.abs();
41
42    let x_bits = x_abs.to_bits();
43    let y_bits = y_abs.to_bits();
44
45    let a_u = x_bits.max(y_bits);
46
47    if a_u >= EXP_MASK_F32 {
48        // x or y is inf or nan
49        if f32::from_bits(x_bits).is_nan() || f32::from_bits(y_bits).is_nan() {
50            return f32::NAN;
51        }
52        if f32::from_bits(x_bits).is_infinite() || f32::from_bits(y_bits).is_infinite() {
53            if f32::from_bits(x_bits).is_infinite() && f32::from_bits(y_bits).is_infinite() {
54                // ∞² - ∞² is undefined
55                return f32::NAN;
56            }
57            return f32::INFINITY;
58        }
59        return f32::from_bits(x_bits);
60    }
61    if x_abs < y_abs {
62        // Would yield sqrt(negative), undefined
63        return f32::NAN;
64    }
65    if x_abs == y_abs {
66        // sqrt(c² - c²) = 0
67        return 0.0;
68    }
69
70    let dx = x as f64;
71    let dy = y as f64;
72
73    #[cfg(any(
74        all(
75            any(target_arch = "x86", target_arch = "x86_64"),
76            target_feature = "fma"
77        ),
78        all(target_arch = "aarch64", target_feature = "neon")
79    ))]
80    {
81        use crate::common::f_fmla;
82        // for FMA environment we're using Kahan style summation which is short and reliable.
83        let w = dy * dy; // RN(bc)
84        let e = f_fmla(-dy, dy, w); // RN(w − bc)
85        let f = f_fmla(dx, dx, -w); // RN(ad − w)
86        let r = e + f; // RN(f + e)
87        let cath = r.sqrt(); // sqrt(x^2 - y^2)
88        cath as f32
89    }
90    #[cfg(not(any(
91        all(
92            any(target_arch = "x86", target_arch = "x86_64"),
93            target_feature = "fma"
94        ),
95        all(target_arch = "aarch64", target_feature = "neon")
96    )))]
97    {
98        use crate::double_double::DoubleDouble;
99        let dy2 = DoubleDouble::from_exact_mult(dy, dy);
100        let fdx = DoubleDouble::from_exact_mult(dx, dx);
101        // element must follow condition |x| > |y| so it always follows fasttwosum requirements
102        let f = DoubleDouble::add_f64(fdx, -dy2.hi).to_f64();
103        let r = dy2.lo + f;
104        let cath = r.sqrt();
105        cath as f32
106    }
107}
108
109#[cfg(test)]
110mod tests {
111    use super::*;
112    #[test]
113    fn test_cathetusf_edge() {
114        assert_eq!(f_cathetusf(5., 3.), 4.);
115        assert_eq!(f_cathetusf(5., 4.), 3.);
116        assert_eq!(f_cathetusf(13., 12.), 5.);
117        assert_eq!(f_cathetusf(65., 16.), 63.);
118        assert_eq!(f_cathetusf(25., 24.), 7.);
119        assert!(f_cathetusf(24., 25.).is_nan());
120    }
121
122    #[test]
123    fn test_cathetusf_edge_cases() {
124        assert_eq!(f_cathetusf(0.0, 0.0), 0.0);
125        assert_eq!(f_cathetusf(f32::INFINITY, 0.0), f32::INFINITY);
126        assert_eq!(f_cathetusf(0.0, f32::INFINITY), f32::INFINITY);
127        assert!(f_cathetusf(f32::INFINITY, f32::INFINITY).is_nan());
128        assert_eq!(f_cathetusf(f32::NEG_INFINITY, 0.0), f32::INFINITY);
129        assert_eq!(f_cathetusf(0.0, f32::NEG_INFINITY), f32::INFINITY);
130        assert!(f_cathetusf(f32::NEG_INFINITY, f32::NEG_INFINITY).is_nan());
131        assert!(f_cathetusf(f32::NAN, 1.0).is_nan());
132        assert!(f_cathetusf(1.0, f32::NAN).is_nan());
133    }
134}
pxfm/triangle/cathetusf.rs

pxfm/triangle/
cathetusf.rs
