cnfy-uint 0.2.3

//! Computes a 64-bit reciprocal of a normalized single limb via Newton-Raphson (MG10 Algorithm 3).
use super::consts::RECIP_TABLE;
use super::KnuthD;

impl KnuthD {
    /// Computes `floor((2^128 - 1) / d) - 2^64` for a normalized `d >= 2^63`.
    ///
    /// Uses MG10 Algorithm 3: table lookup for a ~10-bit seed, then three
    /// Newton-Raphson iterations doubling precision to the full 64 bits.
    /// All arithmetic is wrapping.
    pub(super) fn reciprocal_1(d: u64) -> u64 {
        let d0 = d & 1;
        let d9 = (d >> 55) as usize;
        let d40 = (d >> 24).wrapping_add(1);
        let d63 = d.wrapping_add(1) >> 1;

        let v0 = RECIP_TABLE[d9 - 256] as u64;

        // Newton iteration 1: ~21 bits
        let v1 = (v0 << 11)
            .wrapping_sub(v0.wrapping_mul(v0).wrapping_mul(d40) >> 40)
            .wrapping_sub(1);

        // Newton iteration 2: ~42 bits
        let v2 = (v1 << 13)
            .wrapping_add(v1.wrapping_mul((1u64 << 60).wrapping_sub(v1.wrapping_mul(d40))) >> 47);

        // Newton iteration 3: full 64-bit precision
        let e = ((v2 >> 1) & 0u64.wrapping_sub(d0)).wrapping_sub(v2.wrapping_mul(d63));
        let v3 = ((((v2 as u128) * (e as u128)) >> 64) as u64 >> 1).wrapping_add(v2 << 31);

        // Final adjustment
        let adj = (((v3 as u128) * (d as u128)).wrapping_add(d as u128) >> 64) as u64;
        v3.wrapping_sub(adj).wrapping_sub(d)
    }
}

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

    /// Reciprocal of 2^63 is u64::MAX.
    #[test]
    fn power_of_two() {
        assert_eq!(KnuthD::reciprocal_1(1u64 << 63), u64::MAX);
    }

    /// Reciprocal of u64::MAX is 1.
    #[test]
    fn max() {
        assert_eq!(KnuthD::reciprocal_1(u64::MAX), 1);
    }

    /// Matches naive computation for several values.
    #[test]
    fn matches_naive() {
        for &d in &[
            1u64 << 63,
            (1u64 << 63) + 1,
            u64::MAX,
            u64::MAX - 1,
            0xD555_5555_5555_5555,
            0xAAAA_AAAA_AAAA_AAAA,
        ] {
            let got = KnuthD::reciprocal_1(d);
            let naive = ((u128::MAX / (d as u128)) - (1u128 << 64)) as u64;
            assert_eq!(got, naive, "reciprocal mismatch for d={d:#X}");
        }
    }
}