//! Phase 2 Lehmer mini-Euclidean loop using u128 overflow detection.

use super::LehmerStep;

impl LehmerStep {
    /// Runs the phase 2 Lehmer mini-Euclidean loop with u128 overflow checks.
    ///
    /// Continues Euclidean reduction on `self.a` and `self.b` after the u64
    /// phase has brought `b` below `LEHMER_LIMIT` (2^32). Cofactors may now
    /// exceed 2^32, so new cofactor candidates `nm10` and `nm11` are computed
    /// in u128 to detect overflow before truncating back to u64.
    ///
    /// Typically runs 0–2 steps before the Jebelean condition fires
    /// (`nm10 > b` or `nm11 > r`), indicating the 64-bit approximation can
    /// no longer faithfully represent the full-precision remainders.
    ///
    /// The loop terminates when `b == 0` (GCD found on the approximations)
    /// or when the Jebelean condition triggers.
    #[inline]
    pub fn steps_u128(mut self) -> Self {
        while self.b != 0 {
            let q = self.a / self.b;
            let r = self.a - q * self.b;

            let nm10 = (self.m00 as u128) + (q as u128) * (self.m10 as u128);
            let nm11 = (self.m01 as u128) + (q as u128) * (self.m11 as u128);

            if (nm10 > (self.b as u128)) | (nm11 > (r as u128)) {
                break;
            }

            self.a = self.b;
            self.b = r;
            self.m00 = self.m10;
            self.m01 = self.m11;
            self.m10 = nm10 as u64;
            self.m11 = nm11 as u64;
            self.even = !self.even;
        }
        self
    }
}

#[cfg(test)]
mod ai_tests {
    use super::*;

    /// No steps when b is already zero.
    #[test]
    fn no_steps_when_b_zero() {
        let step = LehmerStep {
            a: 42,
            b: 0,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        let step = step.steps_u128();
        assert_eq!(step.a, 42);
        assert_eq!(step.b, 0);
        assert_eq!(step.m00, 1);
        assert_eq!(step.m01, 0);
        assert_eq!(step.m10, 0);
        assert_eq!(step.m11, 1);
        assert!(step.even);
    }

    /// Single step with small values below LEHMER_LIMIT.
    #[test]
    fn single_step_small() {
        let step = LehmerStep {
            a: 7,
            b: 3,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // q=2, r=1. nm10=1+2*0=1, nm11=0+2*1=2.
        // Jebelean: nm10(1) <= b(3), nm11(2) > r(1) → break after 0 steps.
        let step = step.steps_u128();
        assert_eq!(step.a, 7);
        assert_eq!(step.b, 3);
        assert!(step.even);
    }

    /// Multiple steps with Fibonacci-like small values.
    #[test]
    fn fibonacci_like() {
        let step = LehmerStep {
            a: 89,
            b: 55,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // Fibonacci: q=1 each step. Should take several steps.
        let step = step.steps_u128();
        // Verify matrix determinant invariant: |det| = 1.
        let det = (step.m00 as i128) * (step.m11 as i128) - (step.m01 as i128) * (step.m10 as i128);
        assert!(det == 1 || det == -1, "determinant should be ±1, got {det}");
        if step.even {
            assert_eq!(det, 1);
        } else {
            assert_eq!(det, -1);
        }
    }

    /// Cofactors can grow beyond u32 range in phase 2.
    #[test]
    fn large_cofactors() {
        // Start with cofactors already near 2^31 (as if phase 1 ran first).
        let step = LehmerStep {
            a: 100,
            b: 31,
            m00: 1 << 31,
            m01: (1 << 31) - 1,
            m10: (1 << 31) + 1,
            m11: 1 << 31,
            even: true,
        };
        // q=3, r=7. nm10 = 2^31 + 3*(2^31+1) = 4*2^31 + 3 ≈ 2^33.
        // nm10(~2^33) > b(31) → Jebelean fires immediately.
        let step = step.steps_u128();
        assert_eq!(step.a, 100);
        assert_eq!(step.b, 31);
        assert!(step.even);
    }

    /// Parity toggles correctly across u128 steps.
    #[test]
    fn parity_toggle() {
        let step = LehmerStep {
            a: 13,
            b: 8,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // Step 1: q=1, r=5. nm10=1, nm11=1. 1<=8, 1<=5 → proceed. even=false.
        // Step 2: a=8, b=5. q=1, r=3. nm10=1+0=1, nm11=0+1=1. Wait—
        //   m00=0, m10=1 → nm10=0+1*1=1. m01=1, m11=1 → nm11=1+1*1=2. 1<=5, 2<=3 → proceed. even=true.
        // Step 3: a=5, b=3. q=1, r=2. m00=1, m10=1 → nm10=1+1=2. m01=1, m11=2 → nm11=1+2=3. 2<=3, 3>2 → break.
        let step = step.steps_u128();
        assert!(step.even);
        assert_eq!(step.a, 5);
        assert_eq!(step.b, 3);
    }

    /// Exact GCD terminates when b reaches zero.
    #[test]
    fn terminates_at_zero() {
        let step = LehmerStep {
            a: 6,
            b: 3,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // q=2, r=0. nm10=1+2*0=1, nm11=0+2*1=2. 1<=3, 2<=0? No → break.
        // Actually nm11(2) > r(0) so Jebelean fires before b reaches 0.
        let step = step.steps_u128();
        assert_eq!(step.a, 6);
        assert_eq!(step.b, 3);
    }
}
