Solving Simultaneous Congruences (Chinese Remainder Theorem)

x \equiv 2 \pmod{3}

The equation above is a congruence. What it says is that x % 3 is 2.

The equals sign with three bars means “is equivalent to”, so more literally what the equation says is “x is equivalent to 2, when we are looking at only the integers mod 3”.

5 is a solution, so is 8, so is 11 and so is -1!

The general solution for x is x = 2 +3*N where N is any integer, aka N \in \mathbb{Z}.

What if you were trying to find an integer value x that satisfied the multiple congruences below? How would you solve it?

x \equiv 2 \pmod{3}
x \equiv 2 \pmod{4}
x \equiv 1 \pmod{5}

You could use brute force, but you’d have to check every value from 0 to 60 (not including 60) since 60 is the least common multiple (In this case, you calculate LCM 3*4*5 = 60, since the numbers have a greatest common divisor of 1, aka they are co-prime).

As the number of equations grew, or the mod values were larger, brute force could grow to be a large problem very quickly.

Well luckily there is a better way called the Chinese Remainder Theorem.

BTW – the answer is 26. Or 26 mod 60 more correctly, or 26 + 60*N where N is any integer.

The Chinese Remainder Theorem

The CRT was first published sometime in the 3rd-5th centuries by Sun Tzu – but not the Sun Tzu that wrote “The Art of War”, that was a different person. It was later proven to be true by Gauss in 1801.

When trying to learn the CRT, I found this amazing video that explains it extremely well. You should give it a look: The Chinese Remainder Theorem made easy.

In that video, it talks about modular inversion. You can read about the details of that in the last post I made: Modular Multiplicative Inverse.

An Example

Let’s work through an example that you can follow along with just to make sure you understand it, since I pointed you elsewhere for explanations 😛

Here is our set of congruences:
x \equiv 5 \pmod{7}
x \equiv 8 \pmod{11}
x \equiv 2 \pmod{3}

First up, we want to break it into 3 sections, combined with addition.
x = \_\_\_\_\_ + \_\_\_\_\_ + \_\_\_\_\_

In the above, the first empty space will be the solution to the first equation, the second empty space will be the solution to the second equation, and the third empty space will be the solution to the third equation.

Each section is going to get a co-efficient which is the mods for every equation multiplied together, except the one we are trying to solve. We still have an unknown though per term that we are going to solve in a moment.

x = 11*3*N_1 + 7*3*N_2 + 7*11*N_3
or
x = 33*N_1 + 21*N_2 + 77*N_3

The reason we make those coefficients is because it isolates our terms so that we can solve each equation individually then combine them to get the combined solution.

This works because you can see that in the first term, no matter what we end up choosing for N_1, it will always be a perfect multiple of 3 and 11, which means that the value will be zero in the other terms / other equations, so won’t affect whatever value we come up for them.

Next up, we need to solve for the N’s.

The first term needs to have an N_1 such that 33*N_1 \pmod 7 = 5.

How do we solve that? We could use brute force, testing every number 0 through 6 (otherwise known as 7-1), but in the last post we showed a better way to do it using the extended euclidean algorithm. So, that link once again is: Modular Multiplicative Inverse.

The answer when solving using inversion comes out to be 15. Note that our answer is in “mod 7” space, so you could use any value 1+7*I where I \in \mathbb{Z}, but we can stick to using 15 to make it easier to follow.

That gives us this:
x = 33*15 + 21*N_2 + 77*N_3
or
x = 495 + 21*N_2 + 77*N_3

The value of N2 and N3 end up being -8 and -2 respectively. That gives us:
x = 33*15 + 21*-8 + 77*-2
or
x = 495 - 168 - 154
or
x = 173

Since this is just one of many solutions, the real answer is:
x = 173 \pmod{231}
or
x = 173 + 231*N where N \in \mathbb{Z}

Let’s check our result and see if we got it right…

173 % 7 is indeed 5.
173 % 11 is 8.
173 % 3 is 2.

Woot, it worked!

Example Code

#include <stdio.h>
#include <algorithm>
#include <array>

// x = a mod m
struct SEquation
{
    int a;
    int m;
};
 
//=================================================================================
unsigned int ExtendedEuclidianAlgorithm (int smaller, int larger, int &s, int &t)
{
    // make sure A <= B before starting
    bool swapped = false;
    if (larger < smaller)
    {
        swapped = true;
        std::swap(smaller, larger);
    }
 
    // set up our storage for the loop.  We only need the last two values so will
    // just use a 2 entry circular buffer for each data item
    std::array<int, 2> remainders = { larger, smaller };
    std::array<int, 2> ss = { 1, 0 };
    std::array<int, 2> ts = { 0, 1 };
    int indexNeg2 = 0;
    int indexNeg1 = 1;
 
    // loop
    while (1)
    {
        // calculate our new quotient and remainder
        int newQuotient = remainders[indexNeg2] / remainders[indexNeg1];
        int newRemainder = remainders[indexNeg2] - newQuotient * remainders[indexNeg1];
 
        // if our remainder is zero we are done.
        if (newRemainder == 0)
        {
            // return our s and t values as well as the quotient as the GCD
            s = ss[indexNeg1];
            t = ts[indexNeg1];
            if (swapped)
                std::swap(s, t);
            return remainders[indexNeg1];
        }
 
        // calculate this round's s and t
        int newS = ss[indexNeg2] - newQuotient * ss[indexNeg1];
        int newT = ts[indexNeg2] - newQuotient * ts[indexNeg1];
 
        // store our values for the next iteration
        remainders[indexNeg2] = newRemainder;
        ss[indexNeg2] = newS;
        ts[indexNeg2] = newT;
 
        // move to the next iteration
        std::swap(indexNeg1, indexNeg2);
    }
}
 
//=================================================================================
void WaitForEnter ()
{
    printf("\nPress Enter to quit");
    fflush(stdin);
    getchar();
}
 
//=================================================================================
int main(int argc, char **argv)
{
    const SEquation equations[] =
    {
        { 2, 3 },
        { 2, 4 },
        { 1, 5 }
    };

    const int c_numEquations = sizeof(equations) / sizeof(equations[0]);

    // print out the equations
    printf("Solving for x:\n");
    for (int i = 0; i < c_numEquations; ++i)
        printf("eq %i:  x = %i mod %i\n", i, equations[i].a, equations[i].m);
    printf("\n");

    // make sure the m's are pairwise co-prime
    for (int i = 0; i < c_numEquations; ++i)
    {
        for (int j = i + 1; j < c_numEquations; ++j)
        {
            int s, t;
            int gcd = ExtendedEuclidianAlgorithm(equations[i].m, equations[j].m, s, t);
            if (gcd != 1)
            {
                printf("%i and %i are not co-prime (index %i and %i)\n", equations[i].m, equations[j].m, i, j);
                WaitForEnter();
                return 0;
            }
        }
    }

    // calculate the coefficients for each term
    std::array < int, c_numEquations > coefficients;
    coefficients.fill(1);
    for (int i = 0; i < c_numEquations; ++i)
    {
        for (int j = 0; j < c_numEquations; ++j)
        {
            if (i != j)
                coefficients[i] *= equations[j].m;
        }
    }

    // now figure out how much to multiply each coefficient by to make it have the specified modulus residue (remainder)
    int result = 0;
    for (int i = 0; i < c_numEquations; ++i)
    {
        int s, t;
        ExtendedEuclidianAlgorithm(coefficients[i], equations[i].m, s, t);
        coefficients[i] *= t * equations[i].a;
    }

    // calculate result and simplify it to the smallest positive integer mod lcm
    // lcm is the product when they are pairwise coprime, as the gcd of any two is 1.
    int lcm = 1;
    for (int i = 0; i < c_numEquations; ++i)
    {
        lcm *= equations[i].m;
        result += coefficients[i];
    }
    result = result % lcm;
    if (result < 0)
        result += lcm;

    // print out the answer
    printf("x = %i mod %i\nor...\n", result, lcm);
    printf("x = %i + %i*N\nWhere N is any positive or negative integer.\n\n", result, lcm);

    // verify that our result is correct
    printf("Verifying Equations:\n");
    for (int i = 0; i < c_numEquations; ++i)
        printf("eq %i:  %i mod %i = %i (%s)\n", i, result, equations[i].m, result % equations[i].m, (result % equations[i].m == equations[i].a) ? "PASS" : "!!FAIL!!");

    WaitForEnter();
    return 0;
}

Here’s the output of the program:

Links

One final note. If the program can’t find an answer, it doesn’t necessarily mean that there is no answer. For instance, you could imagine multiplying everything by 2, which would make the modulus values not be co-prime (they would have 2 as a common denominator). That would make this algorithm fail, even though there was a valid answer.

You can try the method of successive substitution as an alternate method when that happens.

Wikipedia: The Chinese Remainder Theorem

Wikipedia: Method Of Successive Substitution

Also, it looks like Khan Academy has a good bit on modular math!

Khan Academy: What is Modular Arithmetic?

Comments

comments

About Demofox

I'm a game and engine programmer at Blizzard Entertainment and have been making games since 1990 (starting out with QBasic and TI-85 games) My shipped titles include: * Heroes of the Storm * StarCraft II: Heart of the Swarm & Legacy of the void * Insanely Twisted Shadow Planet (PC) * Gotham City Impostors (PC, 360, PS3) * Line Rider (PC, Wii, DS) I also like hiking, making music, learning cool new stuff and attempting the impossible.