//! Phase 1 Lehmer mini-Euclidean loop using pure u64 arithmetic.

use super::LehmerStep;

impl LehmerStep {
    /// Runs the phase 1 Lehmer mini-Euclidean loop on 64-bit approximations.
    ///
    /// Performs Euclidean reduction steps on `self.a` and `self.b` (where
    /// `a >= b`) while `b >= 2^32` ([`Self::LEHMER_LIMIT`]). Each step accumulates
    /// the quotient into the 2×2 cofactor matrix (`m00`, `m01`, `m10`, `m11`)
    /// and toggles the swap parity (`even`).
    ///
    /// The loop terminates when `b` drops below the threshold or when the
    /// Jebelean condition fires (`nm10 > b` or `nm11 > r`), indicating the
    /// 64-bit approximation may no longer faithfully represent the full-precision
    /// remainders.
    ///
    /// All arithmetic stays within u64: when `b >= 2^32` the quotients and
    /// cofactors are bounded by `2^31`, so products `q * m` fit in `2^62`.
    #[inline]
    pub fn steps_u64(mut self) -> Self {
        while self.b >= Self::LEHMER_LIMIT {
            let q = self.a / self.b;
            let r = self.a - q * self.b;
            let nm10 = self.m00 + q * self.m10;
            let nm11 = self.m01 + q * self.m11;

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

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

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

    /// Identity matrix when b starts below Self::LEHMER_LIMIT.
    #[test]
    fn no_steps_below_limit() {
        let step = LehmerStep {
            a: 100,
            b: 7,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // b < 2^32 so loop never executes.
        let step = step.steps_u64();
        assert_eq!(step.a, 100);
        assert_eq!(step.b, 7);
        assert_eq!(step.m00, 1);
        assert_eq!(step.m01, 0);
        assert_eq!(step.m10, 0);
        assert_eq!(step.m11, 1);
        assert!(step.even);
    }

    /// Single Euclidean step with large values.
    #[test]
    fn single_step() {
        let step = LehmerStep {
            a: 10_000_000_000,
            b: 7_000_000_000,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // q=1, r=3_000_000_000
        // nm10 = 1 + 1*0 = 1, nm11 = 0 + 1*1 = 1
        // Jebelean: nm10(1) <= b(7B), nm11(1) <= r(3B) → proceed
        // But r=3B < Self::LEHMER_LIMIT(~4.3B) → loop exits after one step.
        let step = step.steps_u64();
        assert_eq!(step.a, 7_000_000_000);
        assert_eq!(step.b, 3_000_000_000);
        assert_eq!(step.m00, 0);
        assert_eq!(step.m01, 1);
        assert_eq!(step.m10, 1);
        assert_eq!(step.m11, 1);
        assert!(!step.even);
    }

    /// Multiple Euclidean steps accumulate the matrix correctly.
    #[test]
    fn multiple_steps() {
        let step = LehmerStep {
            a: 1 << 62,
            b: (1 << 62) - 1,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        let a_orig = step.a;
        let b_orig = step.b;
        // q=1 each time, classic Fibonacci-like reduction.
        let step = step.steps_u64();
        // Verify matrix invariant: m00*m11 - m01*m10 = ±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}");
        // Parity matches determinant sign.
        if step.even {
            assert_eq!(det, 1);
        } else {
            assert_eq!(det, -1);
        }
        // Remainders decreased.
        assert!(step.a < a_orig);
        assert!(step.b < b_orig);
    }

    /// Jebelean condition terminates the loop early.
    #[test]
    fn jebelean_termination() {
        let step = LehmerStep {
            a: (1 << 63) - 1,
            b: 1 << 32,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        let a_orig = step.a;
        // Large quotient: q ≈ 2^31 - 1.
        let step = step.steps_u64();
        // At most 1 step taken (large quotient exhausts approximation).
        assert!(step.a <= a_orig);
    }

    /// Even parity toggles correctly across steps.
    #[test]
    fn parity_toggle() {
        let step = LehmerStep {
            a: 13_000_000_000,
            b: 8_000_000_000,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        // Step 1: q=1, r=5B. nm10=1, nm11=1. Both <= b,r. Proceed. even=false.
        // Step 2: a=8B, b=5B. q=1, r=3B. nm10=1, nm11=2. Both ok. even=true.
        // b=3B < Self::LEHMER_LIMIT → stop.
        let step = step.steps_u64();
        // 2 steps taken → even toggles true→false→true.
        assert!(step.even);
    }
}
