Computing the Chinese Remainder Theorem

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.