1 | use core::u64; |
2 | |
3 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
4 | pub 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 | } |
81 | |