| 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 { (-x, quo) } else { (x, quo) } |
| 95 | } |
| 96 | |
| 97 | #[cfg (test)] |
| 98 | mod 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 | } |
| 107 | |