//! Computes the remainder of a [`U512`] divided by a [`U256`] using Knuth Algorithm D.
use crate::u256::U256;
use crate::u512::U512;

use super::KnuthD;

impl KnuthD {
    /// Computes the remainder of dividing a [`U512`] dividend by this divisor.
    ///
    /// Normalizes the dividend by the same left-shift applied to the divisor
    /// during construction, runs the Knuth Algorithm D quotient loop using
    /// MG10 reciprocal-based estimation, then unnormalizes the remainder.
    pub(crate) fn rem(&self, dividend: &U512) -> U256 {
        let u_le = dividend.to_le_limbs();
        let m = 8 - self.n;

        // Normalize dividend by the same shift used for the divisor.
        let mut un = [0u64; 9];
        self.normalize_dividend(&u_le, &mut un);

        // Run the quotient loop (quotient digits are discarded).
        let mut q_arr = [0u64; 8];
        self.knuth_loop(&mut un[..m + self.n + 1], m, &mut q_arr[..m + 1]);

        // Unnormalize the remainder (already LE) and wrap as U256.
        let rem_le = self.unnormalize(&un);
        U256(rem_le)
    }
}

#[cfg(test)]
mod ai_tests {
    use crate::u256::U256;
    use crate::u512::U512;
    const P: U256 = U256::from_be_limbs([0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFEFFFFFC2F]);

    use super::*;

    /// Helper to compute remainder using KnuthD.
    fn knuth_rem(dividend: &U512, divisor: &U256) -> U256 {
        let n = divisor.sig_limbs();
        let v_le = divisor.to_le_limbs();
        KnuthD::new(&v_le, n).rem(dividend)
    }

    /// Small dividend mod P.
    #[test]
    fn small_dividend() {
        let dividend = U512::from_be_limbs([0, 0, 0, 0, 0, 0, 0, 42]);
        assert_eq!(knuth_rem(&dividend, &P), U256::from_be_limbs([0, 0, 0, 42]));
    }

    /// (P-1)^2 mod P = 1.
    #[test]
    fn p_minus_one_squared() {
        let p_minus_1 = U256::from_be_limbs([
            0xFFFFFFFFFFFFFFFF,
            0xFFFFFFFFFFFFFFFF,
            0xFFFFFFFFFFFFFFFF,
            0xFFFFFFFEFFFFFC2E,
        ]);
        let product = p_minus_1 * p_minus_1;
        assert_eq!(knuth_rem(&product, &P), U256::ONE);
    }

    /// k * P mod P = 0.
    #[test]
    fn multiple_of_p() {
        let small = U256::from_be_limbs([0, 0, 0, 42]);
        let product = small * P;
        assert_eq!(knuth_rem(&product, &P), U256::ZERO);
    }

    /// 2-limb divisor.
    #[test]
    fn two_limb_divisor() {
        let dividend = U512::from_be_limbs([0, 0, 0, 0, 0, 0, 1, 0]);
        let divisor = U256::from_be_limbs([0, 0, 1, 1]);
        assert_eq!(
            knuth_rem(&dividend, &divisor),
            U256::from_be_limbs([0, 0, 1, 0])
        );
    }

    /// (MAX)^2 mod P matches mul_mod.
    #[test]
    fn max_squared() {
        let max = U256::from_be_limbs([u64::MAX; 4]);
        let product = max * max;
        let reduced = U256::from_be_limbs([0, 0, 0, 0x1000003D0]);
        let expected = reduced.mul_mod(&reduced, &P);
        assert_eq!(knuth_rem(&product, &P), expected);
    }
}
