Solving Linear Equations Modulo a Composite Number

Alternative Title: Put a Ring on It: Solving Matrix Equations Over $\mathbb{Z}/N\mathbb{Z}$

Recently, I implemented the index calculus algorithm for computing discrete logs (i.e. solving $a=b^x\text{mod}N$), in order to understand it better. The trickiest part, by far, was figuring out how to solve systems of linear equations modulo $N$, where $N$ is an arbitrary integer which may not be prime. Web search turned up a number of resources, all of which left important details out and didn’t easily translate to code. Well, hopefully the next person who finds themself in the same predicament will stumble across this blog post.

The easy case: when $N$ is prime

The standard method for solving a system of linear equations with real-number (or rational) variables is to convert the system to a matrix equation and solve it using Gaussian elimination or LU decomposition. Here’s an animation of the Gaussian method:

The same methods work just as well if the values in your linear equations are integers modulo some prime. The only difference is whenever you would have divided two real numbers $a$ and $b$, instead, you multiply $a$ by the modular inverse of $b$ ($a{b}^{-1}$). Here’s how that looks when the prime modulus is 7:

Both Gaussian elimination and LU decomposition fail when the modulus $N$ is not prime, because then you can’t “divide”; some values have no inverse, so you can’t take $a{b}^{-1}$.

The slightly harder case: when $N$ has no repeated prime factors

For concreteness, say that $N=6=2\cdot 3$.

You might be aware that if you can find the value of an integer variable modulo $p_1$ and modulo $p_2$, you can combine those two solutions using the Chinese Remainder Theorem to get its value modulo $p_1\cdot p_2$. Here’s the formula; $x$ is the integer, $r_1$ is its remainder modulo $p_1$, and $r_2$ is its remainder modulo $p_2$:

$$x\equiv r_1\cdot p_2\cdot (p_2^{-1}\text{mod}{p}_1)+r_2\cdot p_1\cdot (p_1^{-1}\text{mod}{p}_2)$$

More generally, we can combine any number of remainders with the following formula. The total number of remainders is $k$, and $N$ is the product of all the moduli. (This formula only works if all the moduli are relatively prime, meaning that they don’t share any prime factors.)

$$x\equiv \sum _{i=1}^k{r}_i\frac{N}{p_i}\left({\left(\frac{N}{p_i}\right)}^{-1}\text{mod}{p}_i\right)$$

This immediately suggests: can we solve a linear system mod 6 by finding its solutions mod 2 and mod 3, and combining them using the CRT formula?

Answer: Yes! That works perfectly.

The important thing to note is that the matrix must be non-singular (i.e. there must be a unique solution to the matrix equation) both when reduced mod 3 and when reduced mod 2. This is not automatic; even if the original matrix of mod-6 integers is non-singular, it can become singular when you convert it to a matrix of mod-3 or mod-2 integers.

The still harder case: when $N$ is a squared prime $p^2$

Write the values which we are solving for as $x_{i}p+y_i$, where both $x_i$ and $y_i$ are in $\mathbb{Z}/p\mathbb{Z}$ (their values are between 0 and $p-1$). Find the $y_i$ values by solving the equation mod $p$. (Remember, when solving a linear system modulo a prime, we can use our “ordinary” methods, such as Gaussian elimination or LU decomposition.)

Once we know all the $y_i$ values, we can figure out the next step with a bit of algebra:

$$A\left[\begin{array}{c}{x}_{1}p+y_1\\ x_{2}p+y_2\\ x_{3}p+y_3\\ ⋮\end{array}\right]\equiv b$$

$$pA\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\end{array}\right]+A\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\end{array}\right]\equiv b$$

$$pA\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\end{array}\right]\equiv b-A\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\end{array}\right]$$

Because the left side is multiplied by $p$, the right side must be a multiple of $p$ (mod $p^2$). So we can divide both sides by $p$ (using ordinary integer division, not inverting elements in the modular field). That gives us:

$$A\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\end{array}\right]\equiv \left(b-A\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\end{array}\right]\right)/p$$

All the values on the right-hand side are known, so we can calculate a new right-hand-side “$b$” vector, and we are back to $Ax=b$, where $x$ is a “vector” of values in $\mathbb{Z}/p\mathbb{Z}$.^[1] Solve that matrix equation mod $p$ to get the $x_i$ values.

The even-harder-than-that case: when $N$ is a prime power $p^k$, for arbitrary $k$

Apply the idea of the previous section to solve the linear system mod $p^1$, then mod $p^2$, then mod $p^3$, and so on up to $p^k$.

For example, say that $N=p^3$. Then we can write the values which we are solving for as $x_i{p}^2+y_{i}p+z_i$, where the $x$, $y$, and $z$-values are in $\mathbb{Z}/p\mathbb{Z}$. Our linear system is then:

$$A\left[\begin{array}{c}{x}_1{p}^2+y_{1}p+z_1\\ x_2{p}^2+y_{2}p+z_2\\ x_3{p}^2+y_{3}p+z_3\\ ⋮\end{array}\right]\equiv b$$

To get the $z_i$ values, solve the system mod $p$.

To get the $y_i$ values, solve this equation mod $p$:

$$A\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\end{array}\right]\equiv \left(b-A\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ ⋮\end{array}\right]\right)/p$$

To get the $x_i$ values, solve this equation mod $p$:

$$A\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\end{array}\right]\equiv \left(b-A\left[\begin{array}{c}{y}_{1}p+z_1\\ y_{2}p+z_2\\ y_{3}p+z_3\\ ⋮\end{array}\right]\right)/p^2$$

If $k>3$, do the same thing iteratively all the way up to $k$.

A JavaScript implementation might look something like this, assuming you already have some helper routines for basic vector operations, as well as a linear-solve routine:

function linearSolveModPrimePower(A, b, p, k) {
  let x = linearSolveModPrime(A, b, p); // solve mod p
  let pK = p;
  while (--k > 0) {
    // Lift solution to next higher power of p
    let correction = A.multiplyVector(x);
    let _x = linearSolveModPrime(A, b.sub(correction).scalarDiv(pK), p);
    x = x.add(_x.scalarMul(pK));
    pK *= p;
  }
  return x;
}

The hardest case, which isn’t hard any more if you read the previous sections: when $N$ is an arbitrary composite integer $p_1^{k_1}{p}_2^{k_2}{p}_3^{k_3}\dots$

Need I write anything more? For each prime-power factor $p^k$ of $N$, use the method from the previous section to solve the linear system mod $p^k$. Then use the CRT formula to combine the solutions. That’s it! Of course, this means you do need to factor $N$ before you can start solving linear systems mod $N$.

In JavaScript, the final algorithm looks something like:

function linearSolveModComposite(A, b, N) {
  let x = newZeroVector(A.nCols);
  for (const [p, k] of factorizeInteger(N)) {
    partialSolution = solveLinearSystemModPrimePower(A, b, p, k);
    // Use CRT to combine solutions
    let Mi = N / (p ** k);
    let [xi, yi] = extendedGCD(Mi, p ** k); // invert p**k mod Mi using extended GCD algorithm
    x = x.add(partialSolution.scalarMul(Mi * yi));
  }
  return x;
}

Thanks for reading!

May your matrices always be invertible, and your moduli coprime!

[1] It’s not technically a “vector” (as in a member of a vector space), but you know what I mean. ⏎