Chinese Remainder Theorem

A lot of smart contracts use the SafeMath library. It prevents contracts from having incorrect results, but it does so by failing transactions instead of making them correct. Let’s instead try to do the math correct. In this series I will derive some advanced techniques. The Chinese Remainder Theorem is important in these, so I will introduce it first.

contract ChineseRemainder {
    function chineseRemainder(uint256 x0, uint256 x1)
    public pure returns (uint256 r0, uint256 r1) {
        assembly {
            r0 := x0
            r1 := sub(sub(x1, x0), lt(x1, x0))
        }
    }
}

Special Chinese Remainder reconstruction in Solidity

Sufficiently advanced mathematics is indistinguishable from magic.

— Arthur C. Clarke’s 3rd law (edited)

Warning. This is a very mathematical post. Follow up posts will build on the results derived here, but will not be that math heavy.

Notation. I use square brackets with a subscript for the modulo operation. I do this mostly because the conventional notation is cumbersome for nested expressions. The notation is inspired by Iverson’s notation for the floor and ceiling functions, \floor{a}, and \ceil{b} respectively. In typeset math, the floor, ceil and mod-m operations look like:

\begin{aligned} \floor{a} && \ceil{b} && \mod{c}_m \end{aligned}

Chinese Remainder Theorem

Given two coprime numbers m_0 and m_1, any number x less than m_0 · m_1 can be uniquely written as a pair (x_0, x_1):

\begin{aligned} x_0 &= \mod{x}_{m_0} & x_1 &= \mod{x}_{m_1} \end{aligned}

From this pair, the original number can be reconstructed as:

x = \mod{x_0 ⋅ \mod{m_1^{-1}}_{m_0} ⋅ m_1 + x_1 ⋅ \mod{m_0^{-1}}_{m_1} ⋅ m_0}_{m_0 ⋅ m_1}

This is known as the Chinese Remainder Theorem. In fact, it is a special case with only two moduli, the full theorem can have any number of ms as basis. For completeness, here is the full theorem as I expressed it almost a decade ago:

\begin{aligned} M &= \prod_{[1,n]}^i m_i & \mod{x}_M &= \mod{\sum_{[1,n]}^i \frac{M}{m_i} ⋅ \mod{x ⋅ \frac{m_i}{M} }_{m_i} }_M \end{aligned}

The general Chinese Remainder Theorem.

A special modular basis

The trick I discovered is to use this theorem with a basis m_0 and m_1 such that m_1 = m_0 - 1 and thus m_0 = m_1 + 1. This allows for very easy computation, for example they have nice inverses (you can convince yourself of their correctness by computing m_1 · m_1 and m_0 · 1, expanding the result, and taking the modulus.):

\begin{aligned} \mod{m_1^{-1}}_{m_0} &= m_1 & \mod{m_0^{-1}}_{m_1} &= 1 \end{aligned}

This simplifies the before mentioned reconstruction formula to:

x = \mod{x_0 ⋅ m_1^2 + x_1 ⋅ m_0}_{m_0 ⋅ m_1}

If we substitute m_1 = m_0 - 1 and do a bit of algebra it reduces to:

x = \mod{x_0 + \p{x_1 - x_0} ⋅ m_0}_{m_0^2 - m_0}

The modulo operation now has a large modulus compared to the number inside the brackets. The only time it can ‘wrap around’ is when x_1 < x_0 and then only once. When this happens, we need to add an single extra modulus term, m_0^2 - m_0, to get a positive result again. To get rid of the modulo operation, we introduce a variable c that is 1 when we need to add the modulus and zero otherwise:

x = x_0 + \p{x_1 - x_0} ⋅ m_0 + c ⋅\p{m_0^2 - m_0}

This is the reconstruction formula for any basis where m_1 = m_0 - 1. From here on, I use the specific basis m_0 = 2^{256} and thus m_1 = 2^{256} - 1. This is, in a way, the largest and best basis possible in the EVM. It also turns out to lead to even more simplification in the next step.

Output

We now have the full number x in theory, but this number is to big to represent directly in code. To represent it, we want to split it in the least significant and most significant bits of x. I will call these r_0 and r_1 respectively. In a 256 bit machine like the EVM they are defined as:

\begin{aligned} r_0 &= \mod{x}_{2^{256}} & r_1 &= \floor{\frac{x}{2^{256}}} \end{aligned}

If we substitute the reconstruction formula from above, fill in the modular basis and realize that there is an implicit modulo m_0, almost everything cancels out. (do verify if you are not convinced):

\begin{aligned} r_0 &= x_0 & r_1 &= x_1 - x_0 - c \end{aligned}

As mentioned earlier, c is one if x_1 < x_0 and zero otherwise. That is all there is to it; trivially simple!

Conclusion

So, to summarize, if we know a number modulo 2^{256} and 2^{256} - 1, we can convert that to a 512 bit number with the above two equations. This works iff the number is less than 2^{256} · (2^{256} - 1). In solidity the reconstruction algorithm is as follows:

contract ChineseRemainder {
    function chineseRemainder(uint256 x0, uint256 x1)
    public pure returns (uint256 r0, uint256 r1) {
        r0 = x0;
        r1 = x1 - x0 - (x1 < x0 ? 1 : 0);
    }
}

The Solidity compiler, as of version 0.4.18, does not produce very optimal code here. It introduces a conditional jump that is entirely avoidable. It can be hand-optimized as:

contract ChineseRemainder {
    function chineseRemainder(uint256 x0, uint256 x1)
    public pure returns (uint256 r0, uint256 r1) {
        assembly {
            r0 := x0
            r1 := sub(sub(x1, x0), lt(x1, x0))
        }
    }
}

Two subtractions and one comparisons. This is practically for free in terms of gas!

In the next post, I will use this theorem to build a full 512 bit multiplication in Solidity. It will take just two more instructions, can you guess which?