The Solution of Linear Equations with Integer Coefficients Using Modular Arithmetic and Cramer's Rule

by
Kenneth M. Levasseur
Department of Mathematical Sciences
University of Massachusetts, Lowell
Lowell, MA 01854
levasseuk@woods.uml.edu

Published in
Mathematica in Education and Research
Vol.4 No.1
Copyright 1995 TELOS/Springer-Verlag

Introduction

The solution to a system of n linear equation in n unknowns with integer coefficients is generally rational and not necessarily integer.  For this reason, it would not appear that the solutions to corresponding modular systems would be related.  However, by a clever application of Cramer's Rule, the numerators and denominators of the solution to these systems can be obtained.  The method also can make use of the Chinese Remainder Theorem for large systems.  This method was published in Lipson [Lipson, 1981], and with  Mathematica, it becomes feasible.  With Mathematica's capabilities, this method is not really efficient for ordinary systems, but it's a nice way to combine topics from linear algebra and number theory.  For a truly huge system being solved on a computer with parallel processing capabilities the method outlined here might speed up the solution.   This would be done by replacing a large modulus with many relatively small moduli and then combining the results using the Chinese Remainder Theorem.

The Evaluation of Integer Expressions Using Modular Arithmetic

Let U be an "arithmetic expression," by which is meant a finite number of integer calculations (additions, subtractions, multiplications and divisions) on a set of integers that yields a result upon evaluation that is one or more integer values e(U) = u.  We use e here to denote the evaluation function.  For example U could be the expression for the binomial coefficient "four choose two," ${{4\,3}\over {2\,1}}$ and u = e(U) = 6.

If m is a positive integer such that no divisor of U shares a common factor with m, we let U mod m = R be the corresponding modulo m expression. with result e(R) = r, the mod m result.

The homomorphism properties of the mod function assure us that u and r differ by a multiple of  m; that is,  $u = r + ms$ for some integer s.  When u and r represent more than one value, each of the corresponding values are related in this way.

If u is nonnegative and m is sufficiently large, $m > u$, then  u is equal to r.  If the nonnegativity restriction on u is lifted and m is doubled,  $m > 2|u|$, then u can still be determined from r.   	If  $r < {m\over 2}$ then u = r, and if $r ³ {m\over 2}$ then $u = r - m$  (see Figure 1).

Mathematica's Solve function automatically makes this adjustment when Modulus==m is included.  Outside of Solve a fuction such as is SignMod (defined below) can be used.

SignMod[k_,m_]:=k//Mod[#,m]&//If[2#>m,#-m,#]&;

If u represents more than one value, then the conditions on m would be m > max{u} and m > 2max{u}, respectively.

A Warm-up Example

The binomial coefficient "four choose two" can be computed modulo seven.  All of the numbers involved in the calculation are reduced to their mod seven remainders and all calculations are done modulo seven.  Division by two becomes multiplication by the mod seven inverse of two, which is four.

Print["Numerator = ",modnum=Mod[4 3,7]];
Print["Denominator = ",modden=Mod[2 1,7]];
Print["Result = ",modresult=Mod[modnum PowerMod[modden,-1,7],7]];

Any modulus less than seven is too small. If five is selected as the modulus, the result is one, which is congurent to six modulo five.

Print["Numerator = ",modnum=Mod[4 3,5]];
Print["Denominator = ",modden=Mod[2 1,5]];
Print["Result = ",modresult=Mod[modnum PowerMod[modden,-1,5],5]];

It isn't really necessary that the modulus be prime, as long as the it doesn't share any common factor with the divisors in the calculation.  In this example, the modulus can't be even since the only denominator is two.

Print["Numerator = ",modnum=Mod[4 3,27]];
Print["Denominator = ",modden=Mod[2 1,27]];
Print["Result = ",modresult=Mod[modnum PowerMod[modden,-1,27],27]];

Smaller moduli, even if they alone do not give correct results, provide important information that can be combined using the Chinese Remainder Theorem.  We will consider this idea later.

Cramer's Rule

Cramer's Rule is named after Gabriel Cramer  (1704-1752) but this result is thought to have been first derived by Colin Maclaurin (1698-1746).  It states that if Ax = b is a system of n equations in n unknowns over a field where 

				d= Det[A],		C[i] = A with column i replaced with b, and			c[i] = Det[C[i]] 
				
			then     $x[i]={{c[i]}\over d}$   or equivalently   c[i] = d x[i]

Cramer's Rule applies to any field, including the field of integers modulo p (p a prime).

How will we use Cramer's Rule?  Our plan will be to compute the values of d and the c[i]'s, which are integer values, using modular arithmetic.  Only at the end will we formally write the solution as a quotient as prescribed.  We accomplish this here using Mathematica's Solve function by adding the equation Modulus==m.  We compute the determinant separately using the Det function with the option Modulus->m.  The determinant could be obtained as a byproduct of Gaussian elimination.

Example 1.

Solve:   134 x1 + 53 x2 = -90
             -67 x1 + 88 x2 = 78

A1 ={{134,53},{-67, 88}};
b1 = {-90,78};
x1 = {x[1],x[2]};
system1 = A1.x1==b1;
Thread[system1]//MatrixForm

Notice the use of Thread in this situation to display a system of equations.  The equation system1 is a single vector equation and so its MatrixForm is itself.   Thread will distribute its head of its argument over lists in the argument.  This allows us to equate coordinates in the vector equation and the resulting MatrixForm is more readable.  This illustrated by the following example

FullForm[{s,t}=={u,v}]

Thread[{s,t}=={u,v}]

Returning to the system above, let's look at the solution.

sol1=x1/.Solve[system1,x1][[1]]

As we see here the solution is not a pair of integers, but we can view the process of solving this system as computing four numbers, two numerators and two denominators.

A modular solution:	Using hindsight, we can pick a large enough prime so that one modular solution will give us the solution. We need a prime that is twice as large as the largest integer in the rational solution, which appears to be 15343.

m1=Prime[4000]

37813 is large enough

modsol1=x1/.First[Solve[{system1,Modulus==m1},x1]]

mdet=Det[A1,Modulus->m1]

By Cramer's Rule:

c[1] = SignMod[mdet modsol1[[1]],  m1]

c[2] = SignMod[mdet modsol1[[2]],  m1]

{c[1],c[2]}/mdet

Therefore, the modular solution produces the correct rational solution in this case.

Solution of linear equations with rational coefficients using several moduli

A modular solution to a system of linear equations with integer coefficients with large modulus can be obtained using the Chinese Remainder Theorem and several relatively small primes instead of using one large prime.

Needs["NumberTheory`NumberTheoryFunctions`"]

?ChineseRemainderTheorem

The ring of integers modulo a product of primes is not a field, but all field operations will have been done by the time the small modular solutions are consolidated.

Example 2.

We consider the system

$$  3\,x_1 + 4\,x_2 - x_3 = 8  $$

 $$   x_1 + 5\,x_2 + 2\,x_3 = 3   $$

 $$  2\,x_1 + 4\,x_2 + 8\,x_3 = 11 $$

\end

Without knowing the solution to this system, we assume that the absolute values of all integers involved in the solution are less than 500.   If this is true, we need only select a product of primes that exceeds 1000.  We select 7, 11 and 13, whose product is 1001.

A={{3, 4, -1},{1, 5, 2}, {2, 4, 8}};
b={8,3,11};
x2=Array[x,{3}];
mods={7,11,13};
bigm=Apply[Times,mods];

The function ModSol computes a list consisting of the mod j determinant of the matrix  of coefficients followed by the the mod j solution to the system.

ModSol[j_]:={Det[A,Modulus->j],
         x2/.Solve[{A.x2==b, Modulus==j},x2]}//Flatten

Several modular solutions are put together and the transpose of the list is taken so that the different modular values of the determinant are put together as are the different modular values for each coordinate of the solution.

data=ModSol/@mods//Transpose

The computation above is where parallel computing capabilies could be made use of.  The modular systems can all be solved concurrently and so the time it takes to generate the value of data would be essentially constant as the length of the list mods is increased.

We now know from the first entry of the output above that
           det[A] = 2 (mod 7)
           det[A] = -2 (mod 11)
     and det[A] = -5 (mod 13).

The second entry gives information on the value of x[1] mod 1001. 
           x[1] = 2 (mod 7)
           x[1] = 3 (mod 11)
     and x[1] = -1 (mod 13).
The two remaining entries of the list give similar information for x[2] and x[3]

The ChineseRemainderTheorem  function can can be applied to obtain the values of det[A], x[1] , x[2] and x[3] mod 1001.

bigdata=Map[ChineseRemainderTheorem[#,mods]&,data]

Now we proceed as with a single modulus.

sol2=Map[SignMod[First[bigdata]#,bigm]&,
        Rest[bigdata]]/First[bigdata]

We compare our solution against the solution computed over the rationals and find that they are equal.

x2/.Solve[A.x==b,x2]

A single function to solve systems

The following function summarizes the steps outlined above.

ModularLinearSolve::usage = 
  "ModularLinearSolve[A,B,littlemods] computes the mod m
 solution of A.x==b for x, where m is the product of
  integers in littlemods."
ModularLinearSolve[A_,B_,littlemods_]:=Block[
 {bigm,data,bigdata,k,ModSol,x=Array[s,{Length[B]}]},
    ModSol[j_]:={Det[A,Modulus->j],x/.Solve[{A.x==B,
            Modulus==j},x]}//Flatten;
    data=Map[ModSol, littlemods]//Transpose;
    bigm=Apply[Times,littlemods];
    bigdata=Map[ChineseRemainderTheorem[#,littlemods]&,data];
    Map[SignMod[First[bigdata]#,bigm]&,
        Rest[bigdata]]/First[bigdata]]

ModularLinearSolve[{{2,-1,0},{-1,2,-1},{0,-1,2}},
                   {-2, 5, 3},
                   {13,17,19}]

One problem occurs is when a modulus is selected that  shares a common factor greater than one with the determinant of the system in question.  When this occurs it is possible to salvage some information from the modular calculations, but this has not been built into the function above.   In the example below, the determinant of the system happens to be 47.

ModularLinearSolve[{{7,-2},{-1,7}},{8,4},{47}]

A simpler way to overcome this problem is to expand the function so that for each value in littlemods for which the modular determinant is zero, the data is ignored.  For a reasonably large set of primes, it is highly unlikely that many of the determinants would be zero.

Least squares approximation by powers of x on the interval [0,1].

The least squares approximation of a function on the interval [0,1] with the basis {1, x, . . ., xn-1} is notoriously unstable due to the fact that the successive powers of x are so close to one another.   In general, the best least squares approximation to a function f(x) on [a, b] by linear combinations of functions g1, g2, ..., gn is c1g1 + c2g2 + . . . cngn where the "normal equations"

$$\sum_{j = 1}^{n}\,\left( \int_{a}^{b} \,g_j(t)\,g_k(t) dt \right)c_j 
         =\int_{a}^{b}\,f(t)\,g_k(t) dt $$
 for $1\le k \le n.$
\end

are satisfied.  The matrix of coefficients for this system when on [0, 1] with powers of x is the Hilbert matrix.
For example, to approximate a function with powers of x up to x9 the matrix is the ten by ten Hilbert matrix, H10.  Since H10 has such a miniscule determinant (approx 2 x 10-53) that it is impossible to solve the normal equations using floating point calculations.  We can solve the normal equations over a prime modulus without much difficulty.
  We will approximate f(x) = x10 with lesser powers of x using the least-squares criterion.

Needs["LinearAlgebra`MatrixManipulation`"]

H=HilbertMatrix[10,10];
B=Last[HilbertMatrix[11,10]];

Next we determine the denominators of numbers in each of the normal equations to produce a system of integer equations.

d=Map[(LCM@@Range[#,#+10])&,Range[1,10]]

A formula for the determinant of the Hilbert matrix is well-know (see for example Davis, page 288).  It is  

$$({{{{\left( 1!\,2!\,...\left( n - 1 \right) ! \right) }^4}}\over   {1!\,2!\,3!\,...\left( 2\,n - 1 \right) !}})$$

When n = 10, the determinant is

(Product[k!,{k,1,9}])^4/(Product[k!,{k,1,19}])

Since we will be multiplying the system though by d, the determinant of the system will be

% Apply[Times,d]//N

By experimenting with products of consecutive primes we find that the product of the primes 12 through 18 should be sufficient to solve (d*H).x == d*B.

Times@@Map[Prime,Range[12,18]]//N

c = ModularLinearSolve[d*H,d*B, Map[Prime,Range[12,18]]]

These are the coefficients of the polynomial of best approximation.

y=c.Table[x^k,{k,0,9}];
Short[y]

Naturally, the error is quite small.

Plot[y - x^10, {x,0,1}]

Conclusion

It is generally perceived that Cramer's Rule is purely of theoretical interest (see Axler, page 140).   In a computational introduction to linear algebra it is often given little attention.  This method shows how Cramer's Rule can be applied in a reasonably efficient method for solving systems of equations as opposed to its direct application requiring computation of several determinants.  In combination with the Chinese Remainder Theorem, this method could make use of parallel computation to solve large systems with many relatively small moduli.

References

Sheldon Axler, "Down With Determinants!", The American Mathematical Monthly, 102(1995) 139-154.

Philip J. Davis, Interpolation and Approximation, Dover, 1975.

John D. Lipson, Elements of Algebra and Algebraic Computing, Addison-Welsley, 1981.

About the Author

Ken Levasseur earned a Ph. D. in 1980 from the University of Rhode Island in approximation theory.  Since that time he has been at the University of Massachusetts Lowell.  His main research interests are currently discrete mathematics and applied algebra.

Kenneth M. Levasseur
Department of Mathematical Sciences
University of Massachusetts  Lowell
Lowell, MA  01854
Levasseuk@woods.uml.edu