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:

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:

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):

// 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).