libm/math/
remquo.rs

1#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
2pub fn remquo(mut x: f64, mut y: f64) -> (f64, i32) {
3    let ux: u64 = x.to_bits();
4    let mut uy: u64 = y.to_bits();
5    let mut ex = ((ux >> 52) & 0x7ff) as i32;
6    let mut ey = ((uy >> 52) & 0x7ff) as i32;
7    let sx = (ux >> 63) != 0;
8    let sy = (uy >> 63) != 0;
9    let mut q: u32;
10    let mut i: u64;
11    let mut uxi: u64 = ux;
12
13    if (uy << 1) == 0 || y.is_nan() || ex == 0x7ff {
14        return ((x * y) / (x * y), 0);
15    }
16    if (ux << 1) == 0 {
17        return (x, 0);
18    }
19
20    /* normalize x and y */
21    if ex == 0 {
22        i = uxi << 12;
23        while (i >> 63) == 0 {
24            ex -= 1;
25            i <<= 1;
26        }
27        uxi <<= -ex + 1;
28    } else {
29        uxi &= (!0) >> 12;
30        uxi |= 1 << 52;
31    }
32    if ey == 0 {
33        i = uy << 12;
34        while (i >> 63) == 0 {
35            ey -= 1;
36            i <<= 1;
37        }
38        uy <<= -ey + 1;
39    } else {
40        uy &= (!0) >> 12;
41        uy |= 1 << 52;
42    }
43
44    q = 0;
45
46    if ex + 1 != ey {
47        if ex < ey {
48            return (x, 0);
49        }
50        /* x mod y */
51        while ex > ey {
52            i = uxi.wrapping_sub(uy);
53            if (i >> 63) == 0 {
54                uxi = i;
55                q += 1;
56            }
57            uxi <<= 1;
58            q <<= 1;
59            ex -= 1;
60        }
61        i = uxi.wrapping_sub(uy);
62        if (i >> 63) == 0 {
63            uxi = i;
64            q += 1;
65        }
66        if uxi == 0 {
67            ex = -60;
68        } else {
69            while (uxi >> 52) == 0 {
70                uxi <<= 1;
71                ex -= 1;
72            }
73        }
74    }
75
76    /* scale result and decide between |x| and |x|-|y| */
77    if ex > 0 {
78        uxi -= 1 << 52;
79        uxi |= (ex as u64) << 52;
80    } else {
81        uxi >>= -ex + 1;
82    }
83    x = f64::from_bits(uxi);
84    if sy {
85        y = -y;
86    }
87    if ex == ey || (ex + 1 == ey && (2.0 * x > y || (2.0 * x == y && (q % 2) != 0))) {
88        x -= y;
89        // TODO: this matches musl behavior, but it is incorrect
90        q = q.wrapping_add(1);
91    }
92    q &= 0x7fffffff;
93    let quo = if sx ^ sy { -(q as i32) } else { q as i32 };
94    if sx { (-x, quo) } else { (x, quo) }
95}
96
97#[cfg(test)]
98mod tests {
99    use super::remquo;
100
101    #[test]
102    fn test_q_overflow() {
103        // 0xc000000000000001, 0x04c0000000000004
104        let _ = remquo(-2.0000000000000004, 8.406091369059082e-286);
105    }
106}
libm/math/remquo.rs

libm/math/
remquo.rs
