iron_shapes/
math.rs

1// Copyright (c) 2018-2020 Thomas Kramer.
2// SPDX-FileCopyrightText: 2018-2022 Thomas Kramer
3//
4// SPDX-License-Identifier: AGPL-3.0-or-later
5
6//! Math helper functions.
7
8use num_traits::PrimInt;
9
10/// Implement the fast approximate computation of '1/sqrt(x)' for a type.
11pub trait FastInvSqrt {
12    /// Fast approximate computation of 1/sqrt(x).
13    ///
14    /// See: <https://en.wikipedia.org/wiki/Fast_inverse_square_root>
15    /// And: <http://www.lomont.org/papers/2003/InvSqrt.pdf>
16    fn fast_invsqrt(x: Self) -> Self;
17}
18
19impl FastInvSqrt for f32 {
20    #[inline]
21    fn fast_invsqrt(x: Self) -> Self {
22        fast_invsqrt_32(x)
23    }
24}
25
26impl FastInvSqrt for f64 {
27    #[inline]
28    fn fast_invsqrt(x: Self) -> Self {
29        fast_invsqrt_64(x)
30    }
31}
32
33/// Fast approximate computation of 1/sqrt(x).
34/// The error should be below 0.2%.
35///
36/// See: <https://en.wikipedia.org/wiki/Fast_inverse_square_root>
37/// And: <http://www.lomont.org/papers/2003/InvSqrt.pdf>
38pub fn fast_invsqrt_32(x: f32) -> f32 {
39    let approx = 0x5f375a86u32;
40    let xhalf = 0.5 * x;
41    let i = x.to_bits();
42    let i = approx - (i >> 1);
43    let y = f32::from_bits(i);
44    y * (1.5 - xhalf * y * y)
45}
46
47/// Fast approximate computation of 1/sqrt(x).
48/// The error should be below 0.2%.
49///
50/// See: <https://en.wikipedia.org/wiki/Fast_inverse_square_root>
51/// And: <http://www.lomont.org/papers/2003/InvSqrt.pdf>
52pub fn fast_invsqrt_64(x: f64) -> f64 {
53    let approx = 0x5fe6eb50c7aa19f9u64;
54    let xhalf = 0.5 * x;
55    let i = x.to_bits();
56    let i = approx - (i >> 1);
57    let y = f64::from_bits(i);
58    y * (1.5 - xhalf * y * y)
59}
60
61#[test]
62fn test_fast_invsqrt_32() {
63    for i in 1..1000000 {
64        let x = i as f32;
65        let inv_sqrt_approx = fast_invsqrt_32(x);
66        let inv_sqrt = 1. / x.sqrt();
67        let abs_diff = (inv_sqrt - inv_sqrt_approx).abs();
68        assert!(abs_diff / inv_sqrt < 0.002, "Error should be below 0.2%.")
69    }
70}
71
72#[test]
73fn test_fast_invsqrt_64() {
74    for i in 1..1000000 {
75        let x = i as f64;
76        let inv_sqrt_approx = fast_invsqrt_64(x);
77        let inv_sqrt = 1. / x.sqrt();
78        let abs_diff = (inv_sqrt - inv_sqrt_approx).abs();
79        assert!(abs_diff / inv_sqrt < 0.002, "Error should be below 0.2%.")
80    }
81}
82
83/// Compute square root of integers using Newtons method.
84/// Returns the biggest integer which is smaller or equal to the actual square root of `n`.
85/// Similar to `(i as f64).sqrt().floor() as T` but without type conversions.
86///
87/// See: <https://en.wikipedia.org/wiki/Integer_square_root>
88///
89/// # Panics
90/// Panics when given a negative number.
91///
92/// # Example
93/// ```
94/// use iron_shapes::math::int_sqrt_floor;
95/// assert_eq!(int_sqrt_floor(16), 4);
96/// assert_eq!(int_sqrt_floor(17), 4);
97/// assert_eq!(int_sqrt_floor(24), 4);
98/// assert_eq!(int_sqrt_floor(25), 5);
99/// ```
100pub fn int_sqrt_floor<T: PrimInt>(n: T) -> T {
101    assert!(
102        n >= T::zero(),
103        "Cannot compute the square root of a negative number."
104    );
105    let one = T::one();
106    let two = one + one;
107    let mut x = n;
108    let mut y = (x + one) / two;
109    while y < x {
110        x = y;
111        y = (x + (n / x)) / two;
112    }
113    x
114}
115
116#[test]
117pub fn test_int_sqrt_floor_i32() {
118    for i in 0..100000i32 {
119        let sqrt = int_sqrt_floor(i);
120        let sqrt_reference = (i as f64).sqrt().floor() as i32;
121        assert_eq!(sqrt, sqrt_reference);
122
123        let i_lower = sqrt * sqrt;
124        let i_upper = (sqrt + 1) * (sqrt + 1);
125        assert!(i_lower <= i && i < i_upper);
126    }
127}
iron_shapes/math.rs

iron_shapes/
math.rs
