const_soft_float/soft_f32/
sqrt.rs

1/* origin: Rust libm https://github.com/rust-lang/libm/blob/4c8a973741c014b11ce7f1477693a3e5d4ef9609/src/math/sqrtf.rs */
2/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
3/*
4 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
5 */
6
7use crate::soft_f32::SoftF32;
8use core::cmp::Ordering;
9
10pub(crate) const fn sqrtf(x: SoftF32) -> SoftF32 {
11    const TINY: SoftF32 = SoftF32(1.0e-30);
12
13    let mut z: SoftF32;
14    let sign: i32 = 0x80000000_u32 as i32;
15    let mut ix: i32;
16    let mut s: i32;
17    let mut q: i32;
18    let mut m: i32;
19    let mut t: i32;
20    let mut i: i32;
21    let mut r: u32;
22
23    ix = x.to_bits() as i32;
24
25    /* take care of Inf and NaN */
26    if (ix as u32 & 0x7f800000) == 0x7f800000 {
27        return x.mul(x).add(x); /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
28    }
29
30    /* take care of zero */
31    if ix <= 0 {
32        if (ix & !sign) == 0 {
33            return x; /* sqrt(+-0) = +-0 */
34        }
35        if ix < 0 {
36            return (x.sub(x)).div(x.sub(x)); /* sqrt(-ve) = sNaN */
37        }
38    }
39
40    /* normalize x */
41    m = ix >> 23;
42    if m == 0 {
43        /* subnormal x */
44        i = 0;
45        while ix & 0x00800000 == 0 {
46            ix <<= 1;
47            i = i + 1;
48        }
49        m -= i - 1;
50    }
51    m -= 127; /* unbias exponent */
52    ix = (ix & 0x007fffff) | 0x00800000;
53    if m & 1 == 1 {
54        /* odd m, double x to make it even */
55        ix += ix;
56    }
57    m >>= 1; /* m = [m/2] */
58
59    /* generate sqrt(x) bit by bit */
60    ix += ix;
61    q = 0;
62    s = 0;
63    r = 0x01000000; /* r = moving bit from right to left */
64
65    while r != 0 {
66        t = s + r as i32;
67        if t <= ix {
68            s = t + r as i32;
69            ix -= t;
70            q += r as i32;
71        }
72        ix += ix;
73        r >>= 1;
74    }
75
76    /* use floating add to find out rounding direction */
77    if ix != 0 {
78        z = SoftF32(1.0).sub(TINY); /* raise inexact flag */
79        if ge(z, 1.0) {
80            z = SoftF32(1.0).add(TINY);
81            if gt(z, 1.0) {
82                q += 2;
83            } else {
84                q += q & 1;
85            }
86        }
87    }
88
89    ix = (q >> 1) + 0x3f000000;
90    ix += m << 23;
91    SoftF32::from_bits(ix as u32)
92}
93
94const fn gt(l: SoftF32, r: f32) -> bool {
95    if let Some(ord) = l.cmp(SoftF32(r)) {
96        match ord {
97            Ordering::Greater => true,
98            _ => false,
99        }
100    } else {
101        panic!("Failed to compare values");
102    }
103}
104
105const fn ge(l: SoftF32, r: f32) -> bool {
106    if let Some(ord) = l.cmp(SoftF32(r)) {
107        match ord {
108            Ordering::Less => false,
109            _ => true,
110        }
111    } else {
112        panic!("Failed to compare values");
113    }
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119    use core::f32::*;
120
121    #[test]
122    fn sanity_check() {
123        assert_eq!(sqrtf(SoftF32(100.0)).0, 10.0);
124        assert_eq!(sqrtf(SoftF32(4.0)).0, 2.0);
125    }
126
127    /// The spec: https://en.cppreference.com/w/cpp/numeric/math/sqrt
128    #[test]
129    fn spec_tests() {
130        // Not Asserted: FE_INVALID exception is raised if argument is negative.
131        assert!(sqrtf(SoftF32(-1.0)).0.is_nan());
132        assert!(sqrtf(SoftF32(NAN)).0.is_nan());
133        for f in [0.0, -0.0, INFINITY].iter().copied() {
134            assert_eq!(sqrtf(SoftF32(f)).0, f);
135        }
136    }
137}
const_soft_float/soft_f32/sqrt.rs

const_soft_float/soft_f32/
sqrt.rs
