Computing the Chinese Remainder Theorem

Last year, I wrote a post about the Chinese Remainder Theorem (CRT),
focusing on the math. Here, I want to talk about implementing solvers for
the CRT.

CRT reminder

Assume n_1,dots,n_k are positive integers, pairwise co-prime; that is,
for any ineq j, (n_i,n_j)=1. Let a_1,dots,a_k be
arbitrary integers. The system of congruences with an unknown x:

[begin{align*}
x &equiv a_1 pmod{n_1} \
&vdots \
x &equiv a_k pmod{n_k}
end{align*}]

has a single solution modulo the product
N=n_1times n_2times cdots times n_k.

See the post linked above for a proof of this theorem, along with all the
required number theory prerequisites.

Naive solution by searching exhaustively

Suppose you have an actual programming problem that maps to the CRT. How
would you go about solving it?

To make things more concrete, say the problem is:

[begin{align*}
x &equiv 0 pmod{3} \
x &equiv 3 pmod{4} \
x &equiv 4 pmod{5} \
end{align*}]

The naive solution is search from 1 to 3*4*5-1 (since by the CRT we expect a
unique solution modulo N, which is the product of all n). When we find a
number that satisfies all the congruences – it’s the solution! In this
particular case, 39 is a solution. Since the solution is only unique modulo
N, we can keep adding 60 to our solution to get additional solutions.

Coding this is trivial:

func crtSearch(a, n []int64) int64 {
var N int64 = 1
for _, nk := range n {
N *= nk
}

search:
for i := int64(0); i < N; i++ {
// Does i satisfy all the congruences?
for k := 0; k < len(n); k++ {
if i%n[k] != a[k] {
continue search
}
}
return i
}

return 1
}

This approach works well for small problems, but fails miserably for anything
even moderately large. The problem is obvious: the complexity of this algorithm
is O(kN) where k is the number of congruences. In number theoretic
algorithms, it’s common to talk about the problem size as the bit size of
the numbers involved; in this formulation, this algorithm is exponential (since
N itself is 2 to the power of its size in bits).

For a concrete example, let’s examine a somewhat larger problem:

[begin{align*}
x &equiv 2292 &pmod{77003} \
x &equiv 3010 &pmod{61223} \
x &equiv 500 &pmod{60161} \
x &equiv 399 &pmod{25873} \
end{align*}]

Here the naive algorithm will have to run for up to 4N iterations, where N
is a rather sizable 19 digit number. That’s just not going to cut it [1]; we
need a better approach.

Search by sieving

The best way to understand this algorithm is to sit down with a piece of paper
and a pencil and try to work through a CRT problem by hand. A key insight that
may help is that there is a unique solution for every subset of the CRT problem
as well.

For example, let’s take the first problem (remainders 0, 3, 4 and moduli 3, 4,
5); looking only at the first two congruences, we can find a unique solution
(modulo 12) – in this case 3. By the CRT itself, we know this solution is
unique, and any number in the arithmetic progression 3+12k is also a solution.

Therefore, we can build a solution by induction, starting with the first
congruence and moving to the next each time.

For just the first congruence, a_1 itself is a trivial solution. But so
are all a_1+kn_1, for each integer k. One of these will be a solution
to the second congruence as well. Let’s call this solution x_2; it is a
solution for the first two congruences. We can continue the same approach; since
x_2 is unique modulo n_1n_2, we’ll start looking for a solution
for the third congruence by checking x_2+kn_1n_2 for each integer k.
And so on, until we find a solution to all the congruences.

If you’re having trouble following this explanation with a piece of paper, try
reading the Wikipedia description for “Search by sieving”.

Here is the code that implements this:

func crtSieve(a, n []int64) int64 {
var N int64 = 1
for _, nk := range n {
N *= nk
}

base := a[0]
incr := n[0]

nextBase:
// This loop goes over the congruences one by one; base is a solution
// to the congruences seen so far.
for i := 1; i < len(a); i++ {
// Find a solution that works for the new congruence a[i] as well.
for candidate := base; candidate < N; candidate += incr {
if candidate%n[i] == a[i] {
base = candidate
incr *= n[i]
continue nextBase
}
}
// Inner loop exited without finding candidate
return 1
}
return base
}

By the way, for this approach to be maximally effective we should sort the
congruences by decreasing modulo (solve the congruence with the largest n
first).

Applied to the larger problem presented at the end of the previous section, this
algorithm solves it in half a millisecond on my machine. Not a bad improvement
vs. the ~forever time it takes using the naive approach!

That said, this approach is still exponential! For really large problems (think
public cryptography-level numbers) we’ll need something better.

Using the proof construction

The proof of the CRT includes a
construction of the solution that we could implement in code. A quick reminder:

[x=a_1 N_1 N’_1+a_2 N_2 N’_2+cdots +a_k N_k N’_k]

Is a solution, where N_k=frac{N}{n_k} and N’_k is the
multiplicative inverse of N_k modulo n_k. Finding a
multiplicative modular inverse can be done efficiently with the (extended)
Euclidean algorithm, so we should be good as long as we can handle potentially
enormous numbers.

In Go, this is easy with the math/big package [2]. Here’s an
implementation; unlike the previous ones, it uses big.Int instead of
int64, so it can handle integers of arbitrary size (machine memory
permitting):

func crtConstruct(a, n []*big.Int) *big.Int {
// Compute N: product(n[…])
N := new(big.Int).Set(n[0])
for _, nk := range n[1:] {
N.Mul(N, nk)
}

// x is the accumulated answer.
x := new(big.Int)

for k, nk := range n {
// Nk = N/nk
Nk := new(big.Int).Div(N, nk)

// N’k (Nkp) is the multiplicative inverse of Nk modulo nk.
Nkp := new(big.Int)
if Nkp.ModInverse(Nk, nk) == nil {
return big.NewInt(1)
}

// x += ak*Nk*Nkp
x.Add(x, Nkp.Mul(a[k], Nkp.Mul(Nkp, Nk)))
}
return x.Mod(x, N)
}

Looking at the formula at the beginning of this section, following this code
should be straightforward. The complexity of this algorithm is quadratic, and
it’s much faster than the sieve method on large inputs.

To make the comparison fair, I’ve implemented a version of crtSieve that
uses big.Int as well (since the version shown above is limited by the size
of int64) and ran it vs. crtConstruct on a large-ish CRT problem where
the solution has 144 bits (a 44-digit number in decimal). crtConstruct
was ~20 times faster, in my measurements.

You can see all the code for this post, along with some tests and simple
benchmarks, in this repository.

[1]
The solution, in case you were wondering, is 4412381708627286819.

[2]
I have to admit Go really shines here. Since it comes with an
industrial-strength crypto package which itself relies on arbitrary-sized
integers, math/big has a whole bunch of goodies that are useful for
number-theoretical computations. For example the ModInverse method is
built-in (and if you need something more general, there’s also a GCD
method for computing the full extended Euclidean algorithm).

Flatlogic Admin Templates banner

Leave a Reply

Your email address will not be published. Required fields are marked *