//! Constant-time 512-bit modular reduction for sparse primes.
use super::U512;
use crate::u256::U256;

impl U512 {
    /// Reduces this 512-bit value modulo a sparse prime in constant time.
    ///
    /// Same algorithm as `sparse_reduce_vt` but with all data-dependent
    /// branches removed. Round 2 always executes and the final reduction
    /// uses branchless conditional subtractions via `ct_select`.
    #[inline]
    pub fn sparse_reduce_ct(&self, modulus: &U256, comp: u64) -> U256 {
        let c = comp as u128;

        // Round 1: lo + hi * comp
        let sum = (self.0[0] as u128) + c * (self.0[4] as u128);
        let r0 = sum as u64;
        let carry = sum >> 64;

        let sum = (self.0[1] as u128) + c * (self.0[5] as u128) + carry;
        let r1 = sum as u64;
        let carry = sum >> 64;

        let sum = (self.0[2] as u128) + c * (self.0[6] as u128) + carry;
        let r2 = sum as u64;
        let carry = sum >> 64;

        let sum = (self.0[3] as u128) + c * (self.0[7] as u128) + carry;
        let r3 = sum as u64;
        let overflow = (sum >> 64) as u64;

        // Round 2: always executed (no-op when overflow=0)
        let correction = (overflow as u128) * c;

        let sum = (r0 as u128) + (correction & 0xFFFF_FFFF_FFFF_FFFF);
        let r0 = sum as u64;
        let carry = sum >> 64;

        let sum = (r1 as u128) + (correction >> 64) + carry;
        let r1 = sum as u64;
        let carry = sum >> 64;

        let sum = (r2 as u128) + carry;
        let r2 = sum as u64;
        let carry = sum >> 64;

        let sum = (r3 as u128) + carry;
        let r3 = sum as u64;
        let overflow2 = (sum >> 64) != 0;

        // Branchless final reduction
        let result = U256([r0, r1, r2, r3]);

        let (r1_sub, borrow1) = result.overflowing_sub(modulus);
        let after1 = r1_sub.ct_select(&result, overflow2 | !borrow1);

        let (r2_sub, borrow2) = after1.overflowing_sub(modulus);
        r2_sub.ct_select(&after1, !borrow2)
    }
}

#[cfg(test)]
mod ai_tests {
    use super::*;
    const P: U256 = U256::from_be_limbs([0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFEFFFFFC2F]);

    /// 7^2 = 49.
    #[test]
    fn small_product() {
        let a = U256::from_be_limbs([0, 0, 0, 7]);
        let product = a * a;
        assert_eq!(product.sparse_reduce_ct(&P, 0x1000003D1), U256::from_be_limbs([0, 0, 0, 49]));
    }

    /// Zero product reduces to zero.
    #[test]
    fn zero_product() {
        let product = U512::from_be_limbs([0; 8]);
        assert_eq!(product.sparse_reduce_ct(&P, 0x1000003D1), U256::ZERO);
    }

    /// (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!(product.sparse_reduce_ct(&P, 0x1000003D1), U256::ONE);
    }

    /// CT matches VT for diverse products.
    #[test]
    fn matches_vt() {
        let cases = [
            (U256::from_be_limbs([0, 0, 0, 1]), U256::from_be_limbs([0, 0, 0, 1])),
            (
                U256::from_be_limbs([
                    0x123456789ABCDEF0, 0xFEDCBA9876543210,
                    0x1111222233334444, 0x5555666677778888,
                ]),
                U256::from_be_limbs([
                    0xAAAABBBBCCCCDDDD, 0xEEEEFFFF00001111,
                    0x2222333344445555, 0x6666777788889999,
                ]),
            ),
            (
                U256::from_be_limbs([u64::MAX; 4]),
                U256::from_be_limbs([u64::MAX; 4]),
            ),
        ];
        for (a, b) in &cases {
            let product = *a * *b;
            assert_eq!(
                product.sparse_reduce_ct(&P, 0x1000003D1),
                product.sparse_reduce_vt(&P, 0x1000003D1),
                "CT/VT mismatch for {a:?} * {b:?}"
            );
        }
    }
}
