/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.27/src/triangle/cathetus.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::double_double::DoubleDouble;

/// Computes the missing leg of a right triangle
///
/// Given a hypotenuse `x` and a known leg `y`, returns
/// `sqrt(x^2 - y^2)` = the length of the other leg.
///
/// Domain: requires `|x| >= |y|`. Returns NaN if the input
/// is outside this range.
pub fn f_cathetus(x: f64, y: f64) -> f64 {
    let x_abs = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffffu64);
    let y_abs = f64::from_bits(y.to_bits() & 0x7fff_ffff_ffff_ffffu64);

    let x_bits = x_abs.to_bits();
    let y_bits = y_abs.to_bits();

    let a_u = x_bits.max(y_bits);

    let mut dx = x;
    let mut dy = y;

    const EXP_MASK_F64: u64 = 0x7FF0_0000_0000_0000;
    if a_u >= EXP_MASK_F64 {
        // x or y is inf or nan
        if f64::from_bits(x_bits).is_nan() || f64::from_bits(y_bits).is_nan() {
            return f64::NAN;
        }
        if f64::from_bits(x_bits).is_infinite() || f64::from_bits(y_bits).is_infinite() {
            if f64::from_bits(x_bits).is_infinite() && f64::from_bits(y_bits).is_infinite() {
                // ∞² - ∞² is undefined
                return f64::NAN;
            }
            return f64::INFINITY;
        }
        return f64::from_bits(x_bits);
    }
    if x_abs < y_abs {
        // Would yield sqrt(negative), undefined
        return f64::NAN;
    }
    if x_abs == y_abs {
        // sqrt(c² - c²) = 0
        return 0.0;
    }

    let e_x = x_bits >> 52;
    let e_y = y_bits >> 52;
    let unbiased_e_x = (e_x as i32).wrapping_sub(1023);
    let mut scale = 1f64;

    if e_y == 0 {
        if e_x - e_y > 52 {
            // y is too small to make difference, so result is just |x|
            return x_abs;
        }
        dx *= f64::from_bits(0x6bb0000000000000); // 2^700
        dy *= f64::from_bits(0x6bb0000000000000); // 2^700
        scale = f64::from_bits(0x1430000000000000); // 2^(-700 / 2)
    } else if unbiased_e_x >= 510 {
        dx *= f64::from_bits(0x1430000000000000); // 2^-700
        dy *= f64::from_bits(0x1430000000000000); // 2^-700
        scale = f64::from_bits(0x6bb0000000000000); // 2^(700 / 2)
    } else if unbiased_e_x <= -450 {
        dx *= f64::from_bits(0x6bb0000000000000); // 2^700
        dy *= f64::from_bits(0x6bb0000000000000); // 2^700
        scale = f64::from_bits(0x1430000000000000); // 2^(-700)
    }

    let dy2 = DoubleDouble::from_exact_mult(dy, dy);
    let dx2 = DoubleDouble::from_exact_mult(dx, dx);
    let p = DoubleDouble::sub(dx2, dy2);
    let cath = p.fast_sqrt(); // sqrt(x^2 - y^2)
    cath.to_f64() * scale
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cathethus() {
        assert_eq!(
            f_cathetus(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002248996583584318,
                       0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002842248694776204),
            0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002248996583584318
        );
        assert_eq!(
            f_cathetus(0.00003241747618121237, 0.00003241747618121195),
            5.219099637789996e-12
        );
        assert_eq!(f_cathetus(0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003382112264930946,
                              -0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000005284550413954603),
                   0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003381699384228079);
        assert_eq!(f_cathetus(5., 3.), 4.);
        assert_eq!(f_cathetus(5., 4.), 3.);
        assert_eq!(f_cathetus(13., 12.), 5.);
        assert_eq!(f_cathetus(65., 16.), 63.);
        assert_eq!(f_cathetus(25., 24.), 7.);
        assert!(f_cathetus(24., 25.).is_nan());
    }

    #[test]
    fn test_cathetus_edge_cases() {
        assert_eq!(f_cathetus(0.0, 0.0), 0.0);
        assert_eq!(f_cathetus(f64::INFINITY, 0.0), f64::INFINITY);
        assert_eq!(f_cathetus(0.0, f64::INFINITY), f64::INFINITY);
        assert!(f_cathetus(f64::INFINITY, f64::INFINITY).is_nan());
        assert_eq!(f_cathetus(f64::NEG_INFINITY, 0.0), f64::INFINITY);
        assert_eq!(f_cathetus(0.0, f64::NEG_INFINITY), f64::INFINITY);
        assert!(f_cathetus(f64::NEG_INFINITY, f64::NEG_INFINITY).is_nan());
        assert!(f_cathetus(f64::NAN, 1.0).is_nan());
        assert!(f_cathetus(1.0, f64::NAN).is_nan());
    }
}

Coverage Report

Created: 2026-01-22 07:28

Line	Count	Source
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		*/
29		use crate::double_double::DoubleDouble;
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.
38	0	pub fn f_cathetus(x: f64, y: f64) -> f64 {
39	0	let x_abs = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffffu64);
40	0	let y_abs = f64::from_bits(y.to_bits() & 0x7fff_ffff_ffff_ffffu64);
41
42	0	let x_bits = x_abs.to_bits();
43	0	let y_bits = y_abs.to_bits();
44
45	0	let a_u = x_bits.max(y_bits);
46
47	0	let mut dx = x;
48	0	let mut dy = y;
49
50		const EXP_MASK_F64: u64 = 0x7FF0_0000_0000_0000;
51	0	if a_u >= EXP_MASK_F64 {
52		// x or y is inf or nan
53	0	if f64::from_bits(x_bits).is_nan() \|\| f64::from_bits(y_bits).is_nan() {
54	0	return f64::NAN;
55	0	}
56	0	if f64::from_bits(x_bits).is_infinite() \|\| f64::from_bits(y_bits).is_infinite() {
57	0	if f64::from_bits(x_bits).is_infinite() && f64::from_bits(y_bits).is_infinite() {
58		// ∞² - ∞² is undefined
59	0	return f64::NAN;
60	0	}
61	0	return f64::INFINITY;
62	0	}
63	0	return f64::from_bits(x_bits);
64	0	}
65	0	if x_abs < y_abs {
66		// Would yield sqrt(negative), undefined
67	0	return f64::NAN;
68	0	}
69	0	if x_abs == y_abs {
70		// sqrt(c² - c²) = 0
71	0	return 0.0;
72	0	}
73
74	0	let e_x = x_bits >> 52;
75	0	let e_y = y_bits >> 52;
76	0	let unbiased_e_x = (e_x as i32).wrapping_sub(1023);
77	0	let mut scale = 1f64;
78
79	0	if e_y == 0 {
80	0	if e_x - e_y > 52 {
81		// y is too small to make difference, so result is just \|x\|
82	0	return x_abs;
83	0	}
84	0	dx *= f64::from_bits(0x6bb0000000000000); // 2^700
85	0	dy *= f64::from_bits(0x6bb0000000000000); // 2^700
86	0	scale = f64::from_bits(0x1430000000000000); // 2^(-700 / 2)
87	0	} else if unbiased_e_x >= 510 {
88	0	dx *= f64::from_bits(0x1430000000000000); // 2^-700
89	0	dy *= f64::from_bits(0x1430000000000000); // 2^-700
90	0	scale = f64::from_bits(0x6bb0000000000000); // 2^(700 / 2)
91	0	} else if unbiased_e_x <= -450 {
92	0	dx *= f64::from_bits(0x6bb0000000000000); // 2^700
93	0	dy *= f64::from_bits(0x6bb0000000000000); // 2^700
94	0	scale = f64::from_bits(0x1430000000000000); // 2^(-700)
95	0	}
96
97	0	let dy2 = DoubleDouble::from_exact_mult(dy, dy);
98	0	let dx2 = DoubleDouble::from_exact_mult(dx, dx);
99	0	let p = DoubleDouble::sub(dx2, dy2);
100	0	let cath = p.fast_sqrt(); // sqrt(x^2 - y^2)
101	0	cath.to_f64() * scale
102	0	}
103
104		#[cfg(test)]
105		mod tests {
106		use super::*;
107
108		#[test]
109		fn test_cathethus() {
110		assert_eq!(
111		f_cathetus(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002248996583584318,
112		0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002842248694776204),
113		0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002248996583584318
114		);
115		assert_eq!(
116		f_cathetus(0.00003241747618121237, 0.00003241747618121195),
117		5.219099637789996e-12
118		);
119		assert_eq!(f_cathetus(0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003382112264930946,
120		-0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000005284550413954603),
121		0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003381699384228079);
122		assert_eq!(f_cathetus(5., 3.), 4.);
123		assert_eq!(f_cathetus(5., 4.), 3.);
124		assert_eq!(f_cathetus(13., 12.), 5.);
125		assert_eq!(f_cathetus(65., 16.), 63.);
126		assert_eq!(f_cathetus(25., 24.), 7.);
127		assert!(f_cathetus(24., 25.).is_nan());
128		}
129
130		#[test]
131		fn test_cathetus_edge_cases() {
132		assert_eq!(f_cathetus(0.0, 0.0), 0.0);
133		assert_eq!(f_cathetus(f64::INFINITY, 0.0), f64::INFINITY);
134		assert_eq!(f_cathetus(0.0, f64::INFINITY), f64::INFINITY);
135		assert!(f_cathetus(f64::INFINITY, f64::INFINITY).is_nan());
136		assert_eq!(f_cathetus(f64::NEG_INFINITY, 0.0), f64::INFINITY);
137		assert_eq!(f_cathetus(0.0, f64::NEG_INFINITY), f64::INFINITY);
138		assert!(f_cathetus(f64::NEG_INFINITY, f64::NEG_INFINITY).is_nan());
139		assert!(f_cathetus(f64::NAN, 1.0).is_nan());
140		assert!(f_cathetus(1.0, f64::NAN).is_nan());
141		}
142		}