libm/math/
fmod.rs

1use core::u64;
2
3#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
4pub fn fmod(x: f64, y: f64) -> f64 {
5    let mut uxi = x.to_bits();
6    let mut uyi = y.to_bits();
7    let mut ex = (uxi >> 52 & 0x7ff) as i64;
8    let mut ey = (uyi >> 52 & 0x7ff) as i64;
9    let sx = uxi >> 63;
10    let mut i;
11
12    if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff {
13        return (x * y) / (x * y);
14    }
15    if uxi << 1 <= uyi << 1 {
16        if uxi << 1 == uyi << 1 {
17            return 0.0 * x;
18        }
19        return x;
20    }
21
22    /* normalize x and y */
23    if ex == 0 {
24        i = uxi << 12;
25        while i >> 63 == 0 {
26            ex -= 1;
27            i <<= 1;
28        }
29        uxi <<= -ex + 1;
30    } else {
31        uxi &= u64::MAX >> 12;
32        uxi |= 1 << 52;
33    }
34    if ey == 0 {
35        i = uyi << 12;
36        while i >> 63 == 0 {
37            ey -= 1;
38            i <<= 1;
39        }
40        uyi <<= -ey + 1;
41    } else {
42        uyi &= u64::MAX >> 12;
43        uyi |= 1 << 52;
44    }
45
46    /* x mod y */
47    while ex > ey {
48        i = uxi.wrapping_sub(uyi);
49        if i >> 63 == 0 {
50            if i == 0 {
51                return 0.0 * x;
52            }
53            uxi = i;
54        }
55        uxi <<= 1;
56        ex -= 1;
57    }
58    i = uxi.wrapping_sub(uyi);
59    if i >> 63 == 0 {
60        if i == 0 {
61            return 0.0 * x;
62        }
63        uxi = i;
64    }
65    while uxi >> 52 == 0 {
66        uxi <<= 1;
67        ex -= 1;
68    }
69
70    /* scale result */
71    if ex > 0 {
72        uxi -= 1 << 52;
73        uxi |= (ex as u64) << 52;
74    } else {
75        uxi >>= -ex + 1;
76    }
77    uxi |= (sx as u64) << 63;
78
79    f64::from_bits(uxi)
80}
libm/math/fmod.rs

libm/math/
fmod.rs
