//! Constant-time fused square-and-reduce for sparse primes.
use super::U256;

impl U256 {
    /// Computes `self^2 mod modulus` for sparse primes in constant time.
    ///
    /// Same algorithm as `square_mod_sparse_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 square_mod_sparse_ct(&self, modulus: &U256, comp: u64) -> U256 {
        let a = self.0;

        // === Phase 1: Squaring (10 products via symmetry, identical to VT) ===

        let p = (a[0] as u128) * (a[1] as u128);
        let c1 = p as u64;
        let carry = (p >> 64) as u64;

        let p = (a[0] as u128) * (a[2] as u128) + (carry as u128);
        let c2 = p as u64;
        let carry = (p >> 64) as u64;

        let p = (a[0] as u128) * (a[3] as u128) + (carry as u128);
        let c3 = p as u64;
        let c4 = (p >> 64) as u64;

        let p = (a[1] as u128) * (a[2] as u128) + (c3 as u128);
        let c3 = p as u64;
        let carry = (p >> 64) as u64;

        let p = (a[1] as u128) * (a[3] as u128) + (c4 as u128) + (carry as u128);
        let c4 = p as u64;
        let c5 = (p >> 64) as u64;

        let p = (a[2] as u128) * (a[3] as u128) + (c5 as u128);
        let c5 = p as u64;
        let c6 = (p >> 64) as u64;

        let c7 = c6 >> 63;
        let c6 = (c6 << 1) | (c5 >> 63);
        let c5 = (c5 << 1) | (c4 >> 63);
        let c4 = (c4 << 1) | (c3 >> 63);
        let c3 = (c3 << 1) | (c2 >> 63);
        let c2 = (c2 << 1) | (c1 >> 63);
        let c1 = c1 << 1;

        let sq = (a[0] as u128) * (a[0] as u128);
        let r0 = sq as u64;
        let sum = (c1 as u128) + (sq >> 64);
        let r1 = sum as u64;
        let carry = (sum >> 64) as u64;

        let sq = (a[1] as u128) * (a[1] as u128);
        let sum = (c2 as u128) + (sq as u64 as u128) + (carry as u128);
        let r2 = sum as u64;
        let sum = (c3 as u128) + (sq >> 64) + (sum >> 64);
        let r3 = sum as u64;
        let carry = (sum >> 64) as u64;

        let sq = (a[2] as u128) * (a[2] as u128);
        let sum = (c4 as u128) + (sq as u64 as u128) + (carry as u128);
        let r4 = sum as u64;
        let sum = (c5 as u128) + (sq >> 64) + (sum >> 64);
        let r5 = sum as u64;
        let carry = (sum >> 64) as u64;

        let sq = (a[3] as u128) * (a[3] as u128);
        let sum = (c6 as u128) + (sq as u64 as u128) + (carry as u128);
        let r6 = sum as u64;
        let sum = (c7 as u128) + (sq >> 64) + (sum >> 64);
        let r7 = sum as u64;

        // === Phase 2: Sparse reduce (always both rounds) ===

        let c = comp as u128;

        let sum = (r0 as u128) + c * (r4 as u128);
        let r0 = sum as u64;
        let carry = sum >> 64;

        let sum = (r1 as u128) + c * (r5 as u128) + carry;
        let r1 = sum as u64;
        let carry = sum >> 64;

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

        let sum = (r3 as u128) + c * (r7 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 as u64 as u128);
        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;

        // === Phase 3: 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,
    ]);
    const COMP: u64 = 0x1000003D1;

    /// 7^2 = 49.
    #[test]
    fn small_square() {
        let a = U256::from_be_limbs([0, 0, 0, 7]);
        assert_eq!(a.square_mod_sparse_ct(&P, COMP), U256::from_be_limbs([0, 0, 0, 49]));
    }

    /// 0^2 mod P = 0.
    #[test]
    fn square_zero() {
        assert_eq!(U256::ZERO.square_mod_sparse_ct(&P, COMP), 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,
        ]);
        assert_eq!(p_minus_1.square_mod_sparse_ct(&P, COMP), U256::ONE);
    }

    /// CT matches VT for diverse inputs.
    #[test]
    fn matches_vt() {
        let cases = [
            U256::ZERO,
            U256::ONE,
            U256::from_be_limbs([0, 0, 0, 42]),
            U256::from_be_limbs([
                0x123456789ABCDEF0, 0xFEDCBA9876543210,
                0x1111222233334444, 0x5555666677778888,
            ]),
            U256::from_be_limbs([u64::MAX, u64::MAX, u64::MAX, u64::MAX]),
            U256::from_be_limbs([u64::MAX, u64::MAX, u64::MAX, 0xFFFFFFFEFFFFFC2E]),
        ];
        for a in &cases {
            assert_eq!(
                a.square_mod_sparse_ct(&P, COMP),
                a.square_mod_sparse_vt(&P, COMP),
                "CT/VT mismatch for {a:?}"
            );
        }
    }
}
