The Chinese Remainder Theorem is a tool for solving systems of modular equations.
Suppose we have a system of equations such as,
\[\begin{aligned} x &\equiv a_1 \pmod{n_1} \\ x &\equiv a_2 \pmod{n_2} \\ x &\equiv a_3 \pmod{n_3} \\ &\vdots \\ x &\equiv a_k \pmod{n_k} \\ \end{aligned}\]If for all \(u, v \in \mathbb{N}\) the \(\gcd(n_u, n_v) = 1\), then there's a unique solution for \(x \pmod{\prod_{i=1}^k n_i}\).
Suppose we have the three equations,
\[\begin{aligned} x &\equiv 2 \pmod{3} \\ x &\equiv 2 \pmod{4} \\ x &\equiv 1 \pmod{5} \\ \end{aligned}\]and we want to solve for \(x\). We know that \(x\) is going to be a combination of,
\[ x \equiv x_1 \pmod{3} + x_2 \pmod{4} + x_3 \pmod{5} \]so we start off by adding some initial coefficients that will enable us to cancel out components modulo \(n_i\).
\[ \begin{aligned} x &\equiv 5\cdot 4\cdot x_1 + 5\cdot 3\cdot x_2 + 4 \cdot 3 \cdot x_3 \\ &\equiv 20 \cdot x_1 + 15 \cdot x_2 + 12 \cdot x_3 \\ \end{aligned}\]Notice that since the last term is modulo \(5\), we add a factor of \(5\) to both the first and second terms. This will ensure that when we take the entire equation modulo \(5\), then those two terms will become \(0\), leaving only the third term. This only works because we know all of the \(n\) values are co-prime, meaning none of them are multiples of the others.
If we were to take \(x \pmod{3}\), then the \(x_2\) and \(x_3\) terms would go to \(0\), and \(20 \equiv 2 \pmod{3}\), which fits. Next, we can scale \(15\) by \(2\) since \(30 \equiv 2\pmod{4}\). Finally, if we were to take \(x \pmod{5}\), the first two terms would now go to \(0\), but \(12 \equiv 2 \pmod{5}\), so we need to find a multiple of \(12\) which yields \(1 \pmod{5}\). This is satisfied by scaling \(12\) by \(3\) since \(36 \equiv 1\pmod{5}\).
\[ x \equiv 20 + 30 + 36 = 86 = 26 \pmod{60} \]So, the basic idea is to try and isolate each component by scaling the coefficients such that the other terms will go to zero when we take some modulo \(k\) across all three. Ultimately, we then scale individual terms in order to find something that fits.
Now, this approach is how we can accomplish this by hand, but we incorporated quite a bit of inspection and a bit of trial and error to get the solution. In order to code a solver, we're going to need a step-by-step process which can be applied in the most general way.
It turns out that each componet can be computed with the following formula.
Let \(\mu_i = \prod_{i=1}^k n_i\) and let \(\phi_i\) be the modular inverse of \(\mu_i\) modulo \(n_i\). Therefore,
\[ x_i = a_i \cdot \bigg\lfloor \frac{\mu_i}{n_i} \bigg\rfloor \cdot \phi_i \,,\]and the final solution of \(x\) is,
\[ x = \sum_{i=1}^k x_i \,.\]Implementing this in Julia, we write a function, modsystemsolver, which accepts an array of \(a_i\) values and an array of \(n_i\) values,
where \(q\) is the solution,
\[q = x \mod{\prod_{i=1}^k n_i}\,.\]In our implementation, we initially verify that all of the components of n are relatively prime to each other. We can leverage the gcd function in Julia's standard library which will take an array of integers and compute the greatest common divisor of all of them.
Next, we initialize the variable N with the product of the values of n. Then we iterate over the indices \(1 \leq i \leq |n|\), computing each \(x_i\) value and using the reduce function to accumulate the sum. We finally return the result modulo \(N\).
function modsystemsolver(a, n)
gcd(n) > 1 && throw(
ArgumentError("mod N values are not relatively prime"))
N = prod(n)
result = reduce((acc, i) -> begin
ai, ni = a[i], n[i]
bi = fld(N, ni)
acc + ai * bi * invmod(bi, ni)
end, 1:length(n), init=zero(Int))
mod(result, N)
end