Tags Math , Go

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 i\neq 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_1\times n_2\times \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

  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]

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