//! Applies the accumulated Lehmer 2×2 matrix to a pair of U256 values.

use crate::u256::U256;

use super::LehmerStep;

impl LehmerStep {
    /// Applies the accumulated 2×2 cofactor matrix to a pair `(x, y)`,
    /// returning the transformed pair `(x', y')`.
    ///
    /// When `even` is true the transformation is:
    /// - `x' = m00*x - m01*y`
    /// - `y' = m11*y - m10*x`
    ///
    /// When `even` is false the signs are swapped:
    /// - `x' = m01*y - m00*x`
    /// - `y' = m10*x - m11*y`
    ///
    /// Each product is computed as a U320 (`U256 * u64`), subtracted in 320-bit
    /// arithmetic, and truncated to 256 bits via [`U320::low_u256`]. For
    /// remainder updates the result is guaranteed non-negative and fits in
    /// 256 bits (the Jebelean termination condition ensures this). For Bézout
    /// coefficient updates the truncation implements wrapping 256-bit arithmetic,
    /// which is correct for the unsigned coefficient tracking scheme.
    #[inline]
    pub fn apply(&self, x: U256, y: U256) -> (U256, U256) {
        if self.even {
            (
                (x * self.m00 - y * self.m01).low_u256(),
                (y * self.m11 - x * self.m10).low_u256(),
            )
        } else {
            (
                (y * self.m01 - x * self.m00).low_u256(),
                (x * self.m10 - y * self.m11).low_u256(),
            )
        }
    }
}

#[cfg(test)]
mod ai_tests {
    use crate::u256::U256;

    use super::*;

    /// Identity matrix (m00=1, m01=0, m10=0, m11=1, even=true) leaves values unchanged.
    #[test]
    fn identity_even() {
        let step = LehmerStep {
            a: 0,
            b: 0,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        };
        let x = U256::from_be_limbs([0, 0, 0, 100]);
        let y = U256::from_be_limbs([0, 0, 0, 37]);
        let (rx, ry) = step.apply(x, y);
        assert_eq!(rx, U256::from_be_limbs([0, 0, 0, 100]));
        assert_eq!(ry, U256::from_be_limbs([0, 0, 0, 37]));
    }

    /// Identity matrix with odd parity negates both (wrapping).
    #[test]
    fn identity_odd() {
        let step = LehmerStep {
            a: 0,
            b: 0,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: false,
        };
        // odd: x' = m01*y - m00*x = 0*y - 1*x = -x (wrapping)
        //       y' = m10*x - m11*y = 0*x - 1*y = -y (wrapping)
        let x = U256::from_be_limbs([0, 0, 0, 5]);
        let y = U256::from_be_limbs([0, 0, 0, 3]);
        let (rx, ry) = step.apply(x, y);
        // -5 mod 2^256
        assert_eq!(
            rx,
            U256::ZERO
                .overflowing_sub(&U256::from_be_limbs([0, 0, 0, 5]))
                .0
        );
        // -3 mod 2^256
        assert_eq!(
            ry,
            U256::ZERO
                .overflowing_sub(&U256::from_be_limbs([0, 0, 0, 3]))
                .0
        );
    }

    /// Single Euclidean step matrix: swap with q=1.
    #[test]
    fn single_step_even() {
        // After one Euclidean step with q=1:
        // m00=0, m01=1, m10=1, m11=1, even=false
        // This means: x' = 1*y - 0*x = y, y' = 1*x - 1*y = x - y
        let step = LehmerStep {
            a: 0,
            b: 0,
            m00: 0,
            m01: 1,
            m10: 1,
            m11: 1,
            even: false,
        };
        let x = U256::from_be_limbs([0, 0, 0, 10]);
        let y = U256::from_be_limbs([0, 0, 0, 3]);
        let (rx, ry) = step.apply(x, y);
        assert_eq!(rx, U256::from_be_limbs([0, 0, 0, 3]));
        assert_eq!(ry, U256::from_be_limbs([0, 0, 0, 7]));
    }

    /// Multi-step matrix with even parity applied to large values.
    #[test]
    fn multi_step_even() {
        // Matrix from two q=1 steps:
        // Step 1: [[0,1],[1,1]] (odd)
        // Step 2: [[1,1],[1,2]] (even)
        // even: x' = 1*x - 1*y, y' = 2*y - 1*x
        let step = LehmerStep {
            a: 0,
            b: 0,
            m00: 1,
            m01: 1,
            m10: 1,
            m11: 2,
            even: true,
        };
        let x = U256::from_be_limbs([0, 0, 0, 100]);
        let y = U256::from_be_limbs([0, 0, 0, 38]);
        let (rx, ry) = step.apply(x, y);
        // x' = 1*100 - 1*38 = 62
        assert_eq!(rx, U256::from_be_limbs([0, 0, 0, 62]));
        // y' = 2*38 - 1*100 = -24 (wrapping)
        assert_eq!(
            ry,
            U256::ZERO
                .overflowing_sub(&U256::from_be_limbs([0, 0, 0, 24]))
                .0
        );
    }

    /// Apply with large U256 values exercises the U320 intermediate arithmetic.
    #[test]
    fn large_values() {
        // even, m00=3, m01=1, m10=1, m11=1
        // x' = 3*x - 1*y, y' = 1*y - 1*x
        let step = LehmerStep {
            a: 0,
            b: 0,
            m00: 3,
            m01: 1,
            m10: 1,
            m11: 1,
            even: true,
        };
        let x = U256::from_be_limbs([0, 0, 1, 0]);
        let y = U256::from_be_limbs([0, 0, 2, 0]);
        let (rx, ry) = step.apply(x, y);
        // x' = 3 * 2^64 - 2 * 2^64 = 2^64
        assert_eq!(rx, U256::from_be_limbs([0, 0, 1, 0]));
        // y' = 1 * 2*2^64 - 1 * 2^64 = 2^64
        assert_eq!(ry, U256::from_be_limbs([0, 0, 1, 0]));
    }

    /// Verifies consistency: running steps_u64 then apply reproduces a known GCD step.
    #[test]
    fn steps_then_apply() {
        // gcd(r0=100, r1=37): 100 = 2*37 + 26, 37 = 1*26 + 11, 26 = 2*11 + 4, ...
        // Use Lehmer on 64-bit approximations, then apply to full U256.
        let step = LehmerStep {
            a: 100,
            b: 37,
            m00: 1,
            m01: 0,
            m10: 0,
            m11: 1,
            even: true,
        }
        .steps_u128(); // small values, goes straight to u128 phase
        let r0 = U256::from_be_limbs([0, 0, 0, 100]);
        let r1 = U256::from_be_limbs([0, 0, 0, 37]);
        let (r0, r1) = step.apply(r0, r1);
        // After applying, r0 and r1 should be valid GCD remainders.
        // The matrix reduces (100, 37) by multiple Euclidean steps.
        // Verify r0 > r1 >= 0 and gcd is preserved.
        assert!(r0 >= r1);
    }
}
