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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
//! Constructs a Lehmer step from the current GCD state.
use crate::lehmer_step::LehmerStep;
use super::GcdState;
impl GcdState {
/// Advances the GCD state by one iteration, returning a [`LehmerStep`]
/// matrix when the Lehmer batching produced a non-trivial result.
///
/// Internally extracts the top-bit approximations from `r0` and `r1`,
/// then decides which reduction strategy to use:
///
/// - **`b == 0`**: The approximation of `r1` vanished, so a full-precision
/// Euclidean step ([`full_step`](Self::full_step)) is performed. Returns
/// `None`.
/// - **Identity matrix** (`m01 == 0` after Lehmer steps): The quotient is
/// too large for batching, so a single scalar step
/// ([`scalar_step`](Self::scalar_step)) is performed using the
/// approximate quotient with looped correction. Returns `None`.
/// - **Non-trivial matrix**: The Lehmer loop produced a multi-step cofactor
/// matrix. Returns `Some(lehmer)` so the caller can apply it to both
/// the remainder and coefficient pairs.
///
/// # Examples
///
/// ```
/// use cnfy_uint::gcd_state::GcdState;
/// use cnfy_uint::u256::U256;
/// let P = U256::from_be_limbs([
/// 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF,
/// 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFEFFFFFC2F,
/// ]);
///
/// let mut state = GcdState::new(P, U256::from_be_limbs([0, 0, 0, 7]));
/// // First iteration: P is 256 bits, value is small → b_top == 0 → full_step.
/// assert!(state.lehmer().is_none());
/// ```
#[inline]
pub fn lehmer(&mut self) -> Option<LehmerStep> {
let (a, b) = self.extract_top_bits();
if b == 0 {
self.full_step();
return None;
}
let lehmer = LehmerStep::new(a, b).steps_u64().steps_u128();
if lehmer.m01 == 0 {
self.scalar_step(a / b);
return None;
}
Some(lehmer)
}
}
#[cfg(test)]
mod ai_tests {
use crate::u256::U256;
const P: U256 = U256::from_be_limbs([0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFEFFFFFC2F]);
use super::*;
/// Returns None and performs full_step when b_top is zero.
#[test]
fn full_step_when_b_zero() {
// P is 256 bits, value is tiny → b_top == 0 after extraction.
let mut state = GcdState::new(P, U256::from_be_limbs([0, 0, 0, 7]));
let r0_before = state.r0;
let result = state.lehmer();
assert!(result.is_none());
// full_step advanced the state.
assert_ne!(state.r0, r0_before);
}
/// Returns Some(lehmer) when both values are large and close.
#[test]
fn returns_lehmer_when_both_large() {
let a = U256::from_be_limbs([0x8000_0000_0000_0000, 0, 0, 0]);
let b = U256::from_be_limbs([0x7FFF_FFFF_FFFF_FFFF, 0xFFFF_FFFF_FFFF_FFFF, 0, 0]);
let mut state = GcdState::new(a, b);
let result = state.lehmer();
// Both values are close in magnitude → Lehmer should produce a non-trivial matrix.
assert!(result.is_some());
let lehmer = result.unwrap();
assert_ne!(lehmer.m01, 0);
}
/// Scalar step path when quotient is too large for Lehmer batching.
#[test]
fn scalar_step_on_large_quotient() {
// r0 >> r1: quotient ≈ 2^32, Lehmer can't batch → scalar_step.
let mut state = GcdState::new(
U256::from_be_limbs([0, 0, 0, 1 << 48]),
U256::from_be_limbs([0, 0, 0, 1 << 15]),
);
let r0_before = state.r0;
let result = state.lehmer();
assert!(result.is_none());
assert_ne!(state.r0, r0_before);
}
/// Parity toggles after internal full_step.
#[test]
fn parity_toggles_on_full_step() {
let mut state = GcdState::new(P, U256::from_be_limbs([0, 0, 0, 7]));
assert!(state.even);
state.lehmer();
assert!(!state.even);
}
}