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 are positive integers, pairwise co-prime; that is,
for any , . Let be
arbitrary integers. The system of congruences with an unknown *x*:

has a single solution modulo the product .

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:

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

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, itself is a trivial solution. But so
are all , for each integer *k*. One of these will be a solution
to the second congruence as well. Let's call this solution ; it is a
solution for the first two congruences. We can continue the same approach; since
is unique modulo , we'll start looking for a solution
for the third congruence by checking 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:

Is a solution, where and is the multiplicative inverse of modulo . 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). |