1 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
2 | pub 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 { |
95 | (-x, quo) |
96 | } else { |
97 | (x, quo) |
98 | } |
99 | } |
100 | |
101 | #[cfg (test)] |
102 | mod tests { |
103 | use super::remquo; |
104 | |
105 | #[test ] |
106 | fn test_q_overflow() { |
107 | // 0xc000000000000001, 0x04c0000000000004 |
108 | let _ = remquo(-2.0000000000000004, 8.406091369059082e-286); |
109 | } |
110 | } |
111 | |