1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
//! Divides a 192-bit value by a 128-bit divisor using a precomputed reciprocal (MG10 Algorithm 5).
use super::KnuthD;
impl KnuthD {
/// Divides a 192-bit value `(u2, u1, u0)` by this divisor using the
/// precomputed reciprocal.
///
/// Precondition: `(u2, u1) < (d1, d0)` (ensures quotient fits in 64 bits).
/// Returns `(quotient, remainder)` where remainder is 128 bits.
/// MG10 Algorithm 5.
#[inline]
pub(crate) fn div_3x2(&self, u2: u64, u1: u64, u0: u64) -> (u64, u128) {
// q = u2 * v + (u2, u1)
let q = ((u2 as u128) * (self.v as u128)).wrapping_add(((u2 as u128) << 64) | (u1 as u128));
let q_lo = q as u64;
let q_hi = (q >> 64) as u64;
// Partial remainder: r1 = u1 - q_hi * d1 (wrapping)
let r1 = u1.wrapping_sub(q_hi.wrapping_mul(self.d1));
// t = d0 * q_hi (128-bit)
let t = (self.d0 as u128) * (q_hi as u128);
// r = (r1, u0) - t - d (wrapping 128-bit)
let mut r = (((r1 as u128) << 64) | (u0 as u128))
.wrapping_sub(t)
.wrapping_sub(self.d);
let mut q_final = q_hi.wrapping_add(1);
// First correction: detect overflow via high bits
if (r >> 64) as u64 >= q_lo {
q_final = q_final.wrapping_sub(1);
r = r.wrapping_add(self.d);
}
// Second correction (rare): remainder still >= divisor
if r >= self.d {
q_final = q_final.wrapping_add(1);
r = r.wrapping_sub(self.d);
}
(q_final, r)
}
}
#[cfg(test)]
mod ai_tests {
use super::*;
/// Simple division: (0, 2^62, 0) / 2^127 = quotient 0, remainder 2^126.
#[test]
fn simple() {
let kd = KnuthD::new(&[0, 1u64 << 63, 0, 0], 2);
let (q, r) = kd.div_3x2(0, 1 << 62, 0);
assert_eq!(q, 0);
assert_eq!(r, 1u128 << 126);
}
/// Division with non-zero quotient.
#[test]
fn nonzero_quotient() {
// Divisor = 2^63 (normalized, d0=0)
let kd = KnuthD::new(&[0, 1u64 << 63, 0, 0], 2);
// (0, 2^63, 0) / 2^127 = quotient 1, remainder 0
let (q, r) = kd.div_3x2(0, 1u64 << 63, 0);
// u21 = 2^63, d = 2^127, but 2^63 < 2^127, so div_3x2 is valid
// Result: 2^127 / 2^127 = 1 remainder 0
assert_eq!(q, 1);
assert_eq!(r, 0);
}
/// Division where quotient needs correction.
#[test]
fn with_correction() {
// Use a divisor where corrections are likely: d1=MAX, d0=MAX
let kd = KnuthD::new(&[u64::MAX, u64::MAX, 0, 0], 2);
// Small dividend: (0, 0, 42) / ~2^128 = quotient 0, remainder 42
let (q, r) = kd.div_3x2(0, 0, 42);
assert_eq!(q, 0);
assert_eq!(r, 42);
}
}