Solving Linear Equations Modulo a Composite Number
Alternative Title: Put a Ring on It: Solving Matrix Equations Over
Recently, I implemented the index calculus algorithm for computing discrete logs (i.e. solving ), in order to understand it better. The trickiest part, by far, was figuring out how to solve systems of linear equations modulo , where 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 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 and , instead, you multiply by the modular inverse of (). Here’s how that looks when the prime modulus is 7:
Both Gaussian elimination and LU decomposition fail when the modulus is not prime, because then you can’t “divide”; some values have no inverse, so you can’t take .
The slightly harder case: when has no repeated prime factors
For concreteness, say that .
You might be aware that if you can find the value of an integer variable modulo and modulo , you can combine those two solutions using the Chinese Remainder Theorem to get its value modulo . Here’s the formula; is the integer, is its remainder modulo , and is its remainder modulo :
More generally, we can combine any number of remainders with the following formula. The total number of remainders is , and 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.)
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-12 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 is a squared prime
Write the values which we are solving for as , where both and are in (their values are between 0 and ). Find the values by solving the equation mod . (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 values, we can figure out the next step with a bit of algebra:
Because the left side is multiplied by , the right side must be a multiple of (mod ). So we can divide both sides by (using ordinary integer division, not inverting elements in the modular field). That gives us:
All the values on the right-hand side are known, so we can calculate a new right-hand-side “” vector, and we are back to , where is a “vector” of values in .[1] Solve that matrix equation mod to get the values.
The even-harder-than-that case: when is a prime power , for arbitrary
Apply the idea of the previous section to solve the linear system mod , then mod , then mod , and so on up to .
For example, say that . Then we can write the values which we are solving for as , where the , , and -values are in . Our linear system is then:
To get the values, solve the system mod .
To get the values, solve this equation mod :
To get the values, solve this equation mod :
If , do the same thing iteratively all the way up to .
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 is an arbitrary composite integer
Need I write anything more? For each prime-power factor of , use the method from the previous section to solve the linear system mod . Then use the CRT formula to combine the solutions. That’s it! Of course, this means you do need to factor before you can start solving linear systems mod .
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
// Earlier in this article, the CRT formula for two moduli was shown; this is the version for any number
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. ⏎