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
//! Computes a 2-limb reciprocal for the divisor's top two normalized limbs (MG10 Algorithm 6).
use super::KnuthD;
impl KnuthD {
/// Computes the MG10 2-limb reciprocal for this divisor's top two
/// normalized limbs `(d1, d0)`.
///
/// Bootstraps from [`reciprocal_1`](Self::reciprocal_1) (MG10 Algorithm 3),
/// then adjusts for `d0`'s contribution (MG10 Algorithm 6). The result
/// approximates `floor((2^128 - 1) / (d1, d0)) - 2^64`, enabling fast
/// quotient estimation in [`div_3x2`](Self::div_3x2).
pub(super) fn reciprocal_2(&self) -> u64 {
let mut v = Self::reciprocal_1(self.d1);
// p = d1 * v + d0 (wrapping, track overflow)
let (mut p, overflow) = self.d1.wrapping_mul(v).overflowing_add(self.d0);
if overflow {
v = v.wrapping_sub(1);
if p >= self.d1 {
v = v.wrapping_sub(1);
p = p.wrapping_sub(self.d1);
}
p = p.wrapping_sub(self.d1);
}
// Fine-tune using d0's contribution: t = v * d0 (128-bit)
let t = (v as u128) * (self.d0 as u128);
let t1 = (t >> 64) as u64;
let t0 = t as u64;
let (p2, overflow2) = p.overflowing_add(t1);
if overflow2 {
v = v.wrapping_sub(1);
let combined = ((p2 as u128) << 64) | (t0 as u128);
if combined >= self.d {
v = v.wrapping_sub(1);
}
}
v
}
}
#[cfg(test)]
mod ai_tests {
use super::*;
/// Matches single-limb reciprocal when d0 = 0.
#[test]
fn d0_zero() {
let kd = KnuthD::new(&[0, 1u64 << 63, 0, 0], 2);
assert_eq!(kd.v, KnuthD::reciprocal_1(1u64 << 63));
}
/// Non-zero d0 adjusts the reciprocal downward.
#[test]
fn with_d0() {
let kd = KnuthD::new(&[u64::MAX, 1u64 << 63, 0, 0], 2);
assert!(kd.v <= KnuthD::reciprocal_1(1u64 << 63));
}
/// reciprocal_2 called directly matches stored v.
#[test]
fn matches_stored() {
let kd = KnuthD::new(&[42, 0xC000_0000_0000_0000, 0, 0], 2);
assert_eq!(kd.reciprocal_2(), kd.v);
}
}