Reduced Radix Schoolbook Multiplication modulo p

mul.rs

fn mul(self, rhs: &'b FieldElement5x52) -> FieldElement5x52 {
    const LOW_52_BIT_MASK: u128 = (1 << 52) - 1;
    const LOW_48_BIT_MASK: u128  = (1 << 48) - 1;
    const R: u128 = 0x1000003D10; // R = 2^256 mod p

    #[inline(always)]
    fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) }

    let a: &[u64; 5] = &self.0;
    let b: &[u64; 5] = &rhs.0;

    // Schoolbook multiplication
    let mut l0: u128 = m(a[0], b[0]);
    let mut l1: u128 = m(a[0], b[1]) + m(a[1], b[0]);
    let mut l2: u128 = m(a[0], b[2]) + m(a[1], b[1]) + m(a[2], b[0]);
    let mut l3: u128 = m(a[0], b[3]) + m(a[1], b[2]) + m(a[2], b[1]) + m(a[3], b[0]);
    let mut mid: u128 = m(a[0], b[4]) + m(a[1], b[3]) + m(a[2], b[2]) + m(a[3], b[1]) + m(a[4], b[0]);
    let mut h0: u128 = m(a[1], b[4]) + m(a[2], b[3]) + m(a[3], b[2]) + m(a[4], b[1]);
    let mut h1: u128 = m(a[2], b[4]) + m(a[3], b[3]) + m(a[4], b[2]);
    let mut h2: u128 = m(a[3], b[4]) + m(a[4], b[3]);
    let mut h3: u128 = m(a[4], b[4]);

    // Setup for arithmetic
    let mut out = [0u64; 5];
    let mut carry: u128 = 0;

    // The idea is we multiply the high bits with R and add it to the low bits.
    // Then we set the carries for each with the the high bits beyond the limb.

    // c*2^156
    l3 += (h3 & LOW_52_BIT_MASK) * R;
    out[3] = (l3 & LOW_52_BIT_MASK) as u64;
    h3 >>= 52;
    l3 >>= 52;

    // c*2^208
    mid += ((l3 as u64) as u128) + (((h3 as u64) as u128) * R);
    out[4] = ((mid & LOW_52_BIT_MASK) & LOW_48_BIT_MASK) as u64;
    carry = (mid & LOW_52_BIT_MASK) >> 48;
    mid >>= 52;

    // c*2^0
    h0 += ((mid as u64) as u128);
    l0 += (((((h0 & LOW_52_BIT_MASK) << 4) | carry) as u64) as u128) * (((R as u64) >> 4) as u128);
    out[0] = (l0 & LOW_52_BIT_MASK) as u64;
    h0 >>= 52;
    l0 >>= 52;

    // c*2^52
    h1 += ((h0 as u64) as u128);
    l1 += (((l0 as u64) as u128)) + ((h1 & LOW_52_BIT_MASK) * R);
    out[1] = (l1 & LOW_52_BIT_MASK) as u64;
    h1 >>= 52;
    l1 >>= 52;

    // c*2^104
    h2 += ((h1 as u64) as u128);
    l2 += ((l1 as u64) as u128) + ((h2 & LOW_52_BIT_MASK) * R);
    out[2] = (l2 & LOW_52_BIT_MASK) as u64;
    h2 >>= 52;
    l2 >>= 52;

    // c*2^156
    l3 = ((l2 as u64) as u128) + (((h2 as u64) as u128) * R) + out[3] as u128;
    out[3] = (l3 & LOW_52_BIT_MASK) as u64;
    l3 >>= 52;

    // c*2^208
    out[4] = (((l3 as u64) as u128) + out[4] as u128) as u64;

    FieldElement5x52(out)
}