Field element multiplication with big numbers in Rust

2021/07/26 | by Alex Ren | cryptography

Introduction

We saw in the addition part 1 multiple ways to represent our field elements which is key to optimization.

The final version was 4 x u52 and 1 x 48 to store our 256 bits leaving the extra bits for the carry.

pub struct Fe {
    d: [u64; 5]
}

The long multiplication

Algo with 4 x u64

To make it easier we can start working with 4 u64 digits as we did for the addition, so $a \times b$ :

$\begin{aligned} [t_{7}] & [t_{6}] & [t_{5}] & [t_{4}] & [t_{3}] & [t_{2}] & [t_{1}] & [t_{0}] \\ a_{3} & a_{2} & a_{1} & a_{0} \\ \times & b_{3} & b_{2} & b_{1} & b_{0} \\ b_{0} \cdot a_{3} & b_{0} \cdot a_{2} & b_{0} \cdot a_{1} & b_{0} \cdot a_{0} \\ b_{1} \cdot a_{3} & b_{1} \cdot a_{2} & b_{1} \cdot a_{1} & b_{1} \cdot a_{0} & 0 \\ b_{2} \cdot a_{3} & b_{2} \cdot a_{2} & b_{2} \cdot a_{1} & b_{2} \cdot a_{0} & 0 & 0 \\ b_{3} \cdot a_{3} & b_{3} \cdot a_{2} & b_{3} \cdot a_{1} & b_{3} \cdot a_{0} & 0 & 0 & 0 \\ + c_{6} & + c_{5} & + c_{4} & + c_{3} & + c_{2} & + c_{1} & + c_{0} & 0 \end{aligned}$

I labeled each column with $t_{i}$ with $i := 0 \to 7$ so it’s easier. And finally we also define $c_{i}$ with $i := 0 \to 6$ as the carries for each column.

Since we represent $a$ and $b$ with 4 digits, there are now 8 terms in our result:

Now we have an issue here, because most of these terms can exceed 128 bits while we only have a u128. For example $t_{1}$ $(???)$ is $64 \times 64 + 64 \times 64 + 64 \to 129$ , meaning we would need 129 bits.

Fortunately if we go back to our u52 strategy, we don’t have this issue anymore because at worse $t_{3}$ $(???)$ is $52 \times 52 + 52 \times 52 + 52 \times 52 + 52 \times 52 + 52 \to 106$ i.e. requires 106 bits

Algo with 4 x u52 + u48

Let’s get back to our u52 baby

Now our beautiful multiplication looks more like:

$\begin{aligned} [t_{9}] & [t_{8}] & [t_{7}] & [t_{6}] & [t_{5}] & [t_{4}] & [t_{3}] & [t_{2}] & [t_{1}] & [t_{0}] \\ a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ \times & b_{4} & b_{3} & b_{2} & b_{1} & b_{0} \\ b_{0} \cdot a_{4} & b_{0} \cdot a_{3} & b_{0} \cdot a_{2} & b_{0} \cdot a_{1} & b_{0} \cdot a_{0} \\ b_{1} \cdot a_{4} & b_{1} \cdot a_{3} & b_{1} \cdot a_{2} & b_{1} \cdot a_{1} & b_{1} \cdot a_{0} & 0 \\ b_{2} \cdot a_{4} & b_{2} \cdot a_{3} & b_{2} \cdot a_{2} & b_{2} \cdot a_{1} & b_{2} \cdot a_{0} & 0 & 0 \\ b_{3} \cdot a_{4} & b_{3} \cdot a_{3} & b_{3} \cdot a_{2} & b_{3} \cdot a_{1} & b_{3} \cdot a_{0} & 0 & 0 & 0 \\ b_{4} \cdot a_{4} & b_{4} \cdot a_{3} & b_{4} \cdot a_{2} & b_{4} \cdot a_{1} & b_{4} \cdot a_{0} & 0 & 0 & 0 & 0 \\ + c_{8} & + c_{7} & + c_{6} & + c_{5} & + c_{4} & + c_{3} & + c_{2} & + c_{1} & + c_{0} \end{aligned}$

$\begin{aligned} t_{0} = b_{0} \cdot a_{0} \\ t_{1} = b_{0} \cdot a_{1} + b_{1} \cdot a_{0} + c_{0} \\ t_{2} = b_{0} \cdot a_{2} + b_{1} \cdot a_{1} + b_{2} \cdot a_{0} + c_{1} \\ t_{3} = b_{0} \cdot a_{3} + b_{1} \cdot a_{2} + b_{2} \cdot a_{1} + b_{3} \cdot a_{0} + c_{2} \\ t_{4} = b_{0} \cdot a_{4} + b_{1} \cdot a_{3} + b_{2} \cdot a_{2} + b_{3} \cdot a_{1} + b_{4} \cdot a_{0} + c_{3} \\ t_{5} = b_{1} \cdot a_{4} + b_{2} \cdot a_{3} + b_{3} \cdot a_{2} + b_{4} \cdot a_{1} + c_{4} \\ t_{6} = b_{2} \cdot a_{4} + b_{3} \cdot a_{3} + b_{4} \cdot a_{2} + c_{5} \\ t_{7} = b_{3} \cdot a_{4} + b_{4} \cdot a_{3} + c_{6} \\ t_{8} = b_{4} \cdot a_{4} + c_{7} \\ t_{9} = c_{8} \end{aligned}$

We can also see our result R as $R = R_{1} \cdot 2^{256} + R_{0}$ with:

$R_{0}$ composed of the digits $t_{0}, t_{1}, t_{2}, t_{3}, t_{4}$
$R_{1}$ composed of the digits $t_{5}, t_{6}, t_{7}, t_{8}, t_{9}$

Now because of the size of $R_{1}$ , it’s actually more interesting to reduce our result modulo $P$ directly in the multiplication as opposed to our previous add function in which we could store the carry in the extra bits.

The final purpose is to be able to make multiplications of numbers close to $P$ so the carry $R_{1}$ won’t fit in our 64-bits reserved space.

Similarly to the addition, we have 3 cases:

$R_{1} = 0$ and $R_{0} < P$

$R = R_{0}$

$R_{1} = 0$ , $P <= R_{0} < 2^{256}$

$R = R_{0} - P$

$R_{1} > 0$

$R = R_{0} + R_{1} \cdot (2^{256} - P)$

For more clarity, let’s define $P_{0} = 2^{256} - P = 2^{32} + 977$ , we can now have:

$\begin{aligned} t_{0} = b_{0} \cdot a_{0} + t_{5} \cdot P_{0} \\ t_{1} = b_{0} \cdot a_{1} + b_{1} \cdot a_{0} + c_{0} + t_{6} \cdot P_{0} \\ t_{2} = b_{0} \cdot a_{2} + b_{1} \cdot a_{1} + b_{2} \cdot a_{0} + c_{1} + t_{7} \cdot P_{0} \\ t_{3} = b_{0} \cdot a_{3} + b_{1} \cdot a_{2} + b_{2} \cdot a_{1} + b_{3} \cdot a_{0} + c_{2} + t_{8} \cdot P_{0} \\ t_{4} = b_{0} \cdot a_{4} + b_{1} \cdot a_{3} + b_{2} \cdot a_{2} + b_{3} \cdot a_{1} + b_{4} \cdot a_{0} + c_{3} + t_{9} \cdot P_{0} \end{aligned}$

Cool, now we have new carries to propagate:

$\begin{aligned} t_{0} = b_{0} \cdot a_{0} + t_{5} \cdot P_{0} + c_{4}^{'} \cdot P_{0} \\ t_{1} = b_{0} \cdot a_{1} + b_{1} \cdot a_{0} + c_{0} + t_{6} \cdot P_{0} + c_{0}^{'} + c_{0}^{″} \\ t_{2} = b_{0} \cdot a_{2} + b_{1} \cdot a_{1} + b_{2} \cdot a_{0} + c_{1} + t_{7} \cdot P_{0} + c_{1}^{'} + c_{1}^{″} \\ t_{3} = b_{0} \cdot a_{3} + b_{1} \cdot a_{2} + b_{2} \cdot a_{1} + b_{3} \cdot a_{0} + c_{2} + t_{8} \cdot P_{0} + c_{2}^{'} + c_{2}^{″} \\ t_{4} = b_{0} \cdot a_{4} + b_{1} \cdot a_{3} + b_{2} \cdot a_{2} + b_{3} \cdot a_{1} + b_{4} \cdot a_{0} + c_{3} + t_{9} \cdot P_{0} + c_{3}^{'} + c_{3}^{″} \end{aligned}$

It’s actually way easier to see it this way:

$\begin{aligned} A \times B = & R_{1} \cdot 2^{256} + R_{0} \\ (A \times B) \mod P = & R_{0} + R_{1} \cdot (2^{256} - P) \\ = & R_{0} + R_{1} \cdot P_{0} \end{aligned}$

NOW IN RUST LADIES AND GENTLEMEN

    pub fn mul(&self, b: &Self) -> Self {
        const M52: u128 = 0x000fffffffffffffu128; // 2^52 - 1
        const M48: u128 = 0x0000ffffffffffffu128; // 2^48 - 1
        const P0: u128 = 0x1000003d1u128; // 2^32 + 977

        let (a0, a1, a2, a3, a4) = (self.d[0], self.d[1], self.d[2], self.d[3], self.d[4]);
        let (b0, b1, b2, b3, b4) = (b.d[0], b.d[1], b.d[2], b.d[3], b.d[4]);

        let (
            mut t0, mut t1, mut t2,
            mut t3, mut t4, mut t5,
            mut t6, mut t7, mut t8
        ): (
            u64, u64, u64,
            u64, u64, u64,
            u64, u64, u64
        );
        let t9: u64;

        let mut t: u128;
        let mut c: u128;

        t = a0 as u128 * b0 as u128;
        t0 = (t & M52) as u64;
        t >>= 52;
        t += a0 as u128 * b1 as u128 +
            a1 as u128 * b0 as u128;
        t1 = (t & M52) as u64;
        t >>= 52;
        t += a0 as u128 * b2 as u128 +
            a1 as u128 * b1 as u128 +
            a2 as u128 * b0 as u128;
        t2 = (t & M52) as u64;
        t >>= 52;
        t += a0 as u128 * b3 as u128 +
            a1 as u128 * b2 as u128 +
            a2 as u128 * b1 as u128 +
            a3 as u128 * b0 as u128;
        t3 = (t & M52) as u64;
        t >>= 52;
        t += a0 as u128 * b4 as u128 +
            a1 as u128 * b3 as u128 +
            a2 as u128 * b2 as u128 +
            a3 as u128 * b1 as u128 +
            a4 as u128 * b0 as u128;
        t4 = (t & M52) as u64;
        t >>= 52;
        c = t4 as u128 >> 48;
        t4 &= M48 as u64;

        t += a1 as u128 * b4 as u128 +
            a2 as u128 * b3 as u128 +
            a3 as u128 * b2 as u128 +
            a4 as u128 * b1 as u128;
        t5 = (t & M52) as u64;
        t >>= 52;
        t5 = t5 << 4 | c as u64;
        c = t5 as u128 >> 52;
        t5 &= M52 as u64;

        t += a2 as u128 * b4 as u128 +
            a3 as u128 * b3 as u128 +
            a4 as u128 * b2 as u128;
        t6 = (t & M52) as u64;
        t >>= 52;
        t6 = t6 << 4 | c as u64;
        c = t6 as u128 >> 52;
        t6 &= M52 as u64;
        
        t += a3 as u128 * b4 as u128 +
            a4 as u128 * b3 as u128;
        t7 = (t & M52) as u64;
        t >>= 52;
        t7 = t7 << 4 | c as u64;
        c = t7 as u128 >> 52;
        t7 &= M52 as u64;

        t += a4 as u128 * b4 as u128;
        t8 = (t & M52) as u64;
        t >>= 52;
        t8 = t8 << 4 | c as u64;
        c = t8 as u128 >> 52;
        t8 &= M52 as u64;

        t9 = (t << 4 | c) as u64;

        // 1st reduction R = R1 + R0 * P0
        t = t0 as u128 + t5 as u128 * P0;
        t0 = (t & M52) as u64;
        t >>= 52;

        t += t1 as u128 + t6 as u128 * P0;
        t1 = (t & M52) as u64;
        t >>= 52;

        t += t2 as u128 + t7 as u128 * P0;
        t2 = (t & M52) as u64;
        t >>= 52;

        t += t3 as u128 + t8 as u128 * P0;
        t3 = (t & M52) as u64;
        t >>= 52;

        t += t4 as u128 + t9 as u128 * P0;
        t4 = (t & M52) as u64;
        t >>= 52;
        c = (t4 >> 48) as u128;
        t4 &= M48 as u64;

        c = c | t << 4;

        // 2nd pass
        t = t0 as u128 + c * P0;
        t0 = (t & M52) as u64;
        t >>= 52;

        t += t1 as u128;
        t1 = (t & M52) as u64;
        t >>= 52;

        t += t2 as u128;
        t2 = (t & M52) as u64;
        t >>= 52;

        t += t3 as u128;
        t3 = (t & M52) as u64;
        t >>= 52;

        t += t4 as u128;
        t4 = (t & M48) as u64;

        Self { d: [t0, t1, t2, t3, t4] }
    }

Finally let’s write a quick test:

    #[test]
    fn it_multiply_field_elements() {
        // A=0xfffffffffffffffffffffffffffffffffffffffffffffffffffffbfefffffc2f = p - 2^42
        // B=0xfffffffffffffffffffffffffffffffffffffffffffffffffffff7fefffffc2f = p - 2^43
        let a = Fe::new(
            0xffffffffffffffffu64,
            0xffffffffffffffffu64,
            0xffffffffffffffffu64,
            0xfffffbfefffffc2fu64,
        );
        let b = Fe::new(
            0xffffffffffffffffu64,
            0xffffffffffffffffu64,
            0xffffffffffffffffu64,
            0xfffff7fefffffc2fu64,
        );
        let r = a.mul(&b);
        // r = ((p - 2^42) * (p - 2^43)) % p
        let expected = Fe::new(
            0x0000000000000000u64,
            0x0000000000000000u64,
            0x0000000000200000u64,
            0x0000000000000000u64,
        );
        assert_eq!(r, expected);
    }

and we now have a multiplication.

JUST LIKE THAT

Verbose, but it works.

Optimizations

Carry-saver

Okay, now let’s just think about it, we actually don’t have to propate the carry after the 2nd pass. Same as the addition, we can just use our carry-save storage.

A very simple way to prove it is to work with the worst case scenario: $(P - 1) (P - 1)$

After the 1st pass, it gives us a result $R 0 = {1000003 d 0 fffffffffffffffffffffffffffffffffffffffffffffffefffff 85 dfff 16 f 60}_{16}$ Or $R 0 = 2^{256} \times {1000003 d 0}_{16} + {fffffffffffffffffffffffffffffffffffffffffffffffefffff 85 dfff 16 f 60}_{16}$

Now applying the 2nd pass: $\begin{aligned} R 1 & = {fffffffffffffffffffffffffffffffffffffffffffffffefffff 85 dfff 16 f 60}_{16} + {1000003 d 0}_{16} \times (2^{32} + 977) \end{aligned}$

${1000003 d 0}_{16} \times (2^{32} + 977)$ is 65 bits so our first digit needs 66 bits. We only have 12 extra bits on the 1st digit $d_{0}$ , so we need to propagate the carry on the 2nd digit $d_{1}$

And we don’t have to go further for now, we can still work with few operations before we overflow our 320 bits storage.

Let’s remove the last carry propagation on t3 and t4:

        // ...
        // 2nd pass
        t = t0 as u128 + c * P0;
        t0 = (t & M52) as u64;
        t >>= 52;

        t += t1 as u128;
        t1 = (t & M52) as u64;

        Self { d: [t0, t1, t2, t3, t4] }

Combine operations

Another optimization would be to combine the operations on $t_{0}$ , we have:

        // 1st pass
        t = t0 as u128 + t5 as u128 * P0;
        t0 = (t & M52) as u64;
        // ...
        // 2nd pass
        t = t0 as u128 + c * P0;
        t0 = (t & M52) as u64;

Could be nice to combine those just so we save some instructions, better to do $t_{0} + (a + b) P_{0}$ than $t_{0} + a P_{0} + b P_{0}$

We can use the fact that we don’t need to propagate the carry over $t_{3}$ and $t_{4}$ on the 2nd pass. Therefore, we can work with them first.

Final version

    pub fn mul(&mut self, b: &Self) -> Self {
        const M52: u128 = 0x000fffffffffffffu128; // 2^52 - 1
        const M48: u64 = 0x0000ffffffffffffu64; // 2^48 - 1
        const P0: u128 = 0x1000003d1u128; // 2^32 + 977
        const P1: u128 = 0x1000003d10u128; // 2^32 + 977 << 4

        let (a0, a1, a2, a3, a4) = (
            self.d[0] as u128, self.d[1] as u128, self.d[2] as u128,
            self.d[3] as u128, self.d[4] as u128
        );
        let (b0, b1, b2, b3, b4) = (
            b.d[0] as u128, b.d[1] as u128, b.d[2] as u128,
            b.d[3] as u128, b.d[4] as u128
        );
        let mut tx: u128;
        let mut cx: u128;
        let (t0, t1, t2, mut t3, mut t4, mut t5): (u64, u64, u64, u64, u64, u128);
        let c4: u64;

        // t3
        tx = a0 * b3 + a1 * b2 + a2 * b1 + a3 * b0;
        // t8
        cx = a4 * b4;
        // t3 + t8 * P1
        tx += (cx & M52) * P1;
        cx >>= 52;
        t3 = (tx & M52) as u64;
        tx >>= 52;

        // t4
        tx += a0 * b4 + a1 * b3 + a2 * b2 + a3 * b1 + a4 * b0;
        // (c3 + t4) + (c8 + t9) * P1
        tx += cx * P1;
        t4 = (tx & M52) as u64;
        tx >>= 52;
        c4 = t4 >> 48;
        t4 &= M48;

        // t5
        cx = tx + a1 * b4 + a2 * b3 + a3 * b2 + a4 * b1;
        // t0
        tx = a0 * b0;
        t5 = cx & M52;
        cx >>= 52;
        t5 = (t5 << 4) | c4 as u128;
        // c9 + t0 + (c4 + t5) * P0
        tx += t5 * P0;
        t0 = (tx & M52) as u64;
        tx >>= 52;

        // t1
        tx += a0 * b1 + a1 * b0;
        // t6
        cx += a2 * b4 + a3 * b3 + a4 * b2;
        // c0 + t1 + (c5 + t6) * P1
        tx += (cx & M52) * P1;
        cx >>= 52;
        t1 = (tx & M52) as u64;
        tx >>= 52;

        // t2
        tx += a0 * b2 + a1 * b1 + a2 * b0;
        // t7
        cx += a3 * b4 + a4 * b3;
        // t2 + t7 * P1
        tx += (cx & M52) * P1;
        cx >>= 52;
        t2 = (tx & M52) as u64;
        tx >>= 52;

        // t23
        tx += cx * P1 + t3 as u128;
        t3 = (tx & M52) as u64;
        tx >>= 52;
        // t24
        tx += t4 as u128;
        t4 = tx as u64;

        Self { d: [t0, t1, t2, t3, t4] }
    }

#cryptography #rust #math