Solution of Systems of Linear Diophantine Equations ^{⋆}

Gennadi I. Malaschonok

Tambov State University, 33 Internatsionalnaya, 392622 Tambov, Russia

                      malaschonok@math-iu.tambov.su

Abstract. Two new methods to solve linear systems of Diophantine equations are proposed – modular (CRT) and p-adic (Hensel). Each of them allows to obtain solutions of a system with the size n × m with the complexity O(n^{β}m). For quasi-square systems, the p-adic method allows to get solution with the complexity O(n^{3}), and the modular method

– O(n^{β+1}). Both estimations have the accuracy up to the logarithmic multipliers, β being the power in the estimation of matrix multiplication time.

1 Introduction

The usual method to solve systems of linear Diophantine equations is based on the computation of Smith form of the coefficient matrix [6], [5]. In the recent papers [11], [10] the probabilistic methods for finding the solutions of Diophantine systems were obtained. The main lack of these methods is the probabilistic character of the solution and an absence of a capability to prove the incompatibility of a system. These methods do not allow to distinguish a compatible system, in which the solution is not obtained in the process of solving, from an incompatible system.

The main result of the article is the following.

There are offered deterministic p-adic and modular methods which are independent of the Smith form computation. Both methods reduce the problem to the solving of systems of linear equations in the residue ring.

The p-adic method is based on the linear p-adic Hensel lifting. This is a way to obtain the basic set of system solutions on the plane of all solutions. Then we compute all integer points on this plane.

The second method is based on the Chinese remainder theorem. The result of its application is the diagonal form of a non-singular submatrix of maximal rank. Then we solve the system in modulus which equals the determinant of this submatrix.

The complexity of the solution is estimated for systems in the ring of integers Z and in the ring of polynomials F [x] over a field F.

⋆ in: Computer Algebra in Scientific Computing - CASC’01, Springer, Berlin, 2001

For essentially rectangular systems, when m > 2n, both methods have the same order of the complexity – O(n^{β+1}m) within a logarithmic factor. For quasisquare systems the complexity of a p-adic method is O(n^{3}) and the complexity of a modular method is O(n^{β+1}). Here β is the power in the estimation of complexity of matrix multiplication. Known for today value of β is less than 2.356 [3].

Next section has an auxiliary character, it is devoted to the problem of solving systems of linear equations in a residue ring. Here is presented a method with complexity O(n^{β}) ring operations. For a system of the order n the asymptotic evaluation of this method is less than 6.54M(n), where M(n) is the asymptotic evaluation of Strassen multiplication of matrices of order n.

2 QuasiEuclidean rings

Let R be a commutative ring with identity, α ∈ R^{2} – be a non-zero pair. An unimodular matrix E_{α} ∈ R^{2×2} is called the Diophantine matrix of α if

E_{α}α = ( \begin{matrix}g \\ 0 \end{matrix} ) , g ≠ 0, det E_{α} = 1.

The unit matrix is called the Diophantine matrix of a zero pair.

We must select a special class of the rings, in which an algorithm for computing the Diophantine matrix exists.

Definition 1 . A ring R is called quasiEuclidean, if an algorithm for computing the Diophantine matrix for any α ∈ R^{2} exists.

If α = ( \begin{matrix}p \\ q \end{matrix} ), then ( \begin{matrix}p \\ q \end{matrix} ) = E_{α}^{−1} ( \begin{matrix}g \\ 0 \end{matrix} ). So the ideal (p, q) is a principal ideal

generated by g. Therefore each of the quasiEuclidean rings is a ring of principal ideals, where the algorithm for computing of ideal generator exists.

To the contrary. Let R be the ring of principal ideals, in which an algorithm for calculation of the generator g of an ideal (p, q), p, q ∈ R exists, denote it GI algorithm:

{g, a, b} = GI(p, q), (g) = (p, q), ap + bq = g.

For a non-zero pair α = ( \begin{matrix}p \\ q \end{matrix} ), a Diophantine matrix is E_{α} = ( \begin{matrix}a & b \\ −q/g & p/g \end{matrix} ) .

Therefore a commutative ring is quasiEuclidean if and only if it is a ring of principal ideals, where an algorithm for computing of ideal generator exists.

Note that the numbers a, b and g are not uniquely defined. Therefore the Diophantine matrix is not uniquely defined. The Diophantine matrix can be determined uniquely if an order is introduced into R.

If R is the Euclidean ring, then the extended Euclidean algorithm may be taken as GI algorithm.

Proposal 1 . QuasiEuclidivity of a ring is inherited by its quotient-ring.

A quatient-ring of the ring of principal ideal is a ring of principal ideal. So to prove this proposal we need to show how to construct the Diophantine matrices in quatient-ring.

Let R′ be a quotient-ring of a quasiEuclidean ring R with some order relation, and

ϕ : R → R′

be a canonical homomorphism. Let us introduce an order relation in R, and define a function

ψ : R′ → R,

selecting a representative of a coset for p′ ∈ R′ in a pre-image ϕ^{−1}(p′) by a

standard way: ψ(p′) = min(ϕ^{−1}(p′)). Let ψ(p′) = p, ψ(q′) = q, α = ( \begin{matrix}p \\ q \end{matrix} ), E_{α} –

Diophantine matrix for α and E_{α}α = ( \begin{matrix}g \\ 0 \end{matrix} ).

For a non-zero pair α′ = ( \begin{matrix}p′ \\ q′ \end{matrix} ) ∈ R′^{2} we define a Diophantine matrix E_{α′}

as follows: E_{α′} = ϕ(E_{α}). Such definition is correct since

E_{α′}α′ = ( \begin{matrix}g′ \\ 0 \end{matrix} ) , det E_{α′} = 1, g′ ≠ 0.

Let’s show, that g′ = ϕ(g) ≠ 0. Due to the invertibility of Diophantine matrices, p and q are divided by g. Therefore if ϕ(g) = 0, then ϕ(p) = ϕ(q) = 0. But by the condition, at least one of the numbers p′ = ϕ(p), q′ = ϕ(q) is not a zero.

Therefore a set of quasiEuclidean rings contains all Euclidean rings and their quotient-rings.

Let now R be a quasiEuclidean ring, A ∈ R^{n×m}. A matrix A is called upper triangular, if all elements under the main diagonal are zeros: (i, j) = 0 for all n ≥ j > i ≥ 1. (We say nothing about other matrix elements.)

We say that the matrix A is decomposable, if it is a product of an unimodular matrix Ã and an upper triangular matrix Ā: A = ÃĀ.

Proposal 2 . In a quasiEuclidean ring any matrix is decomposable.

Proof. We compute such decomposition for a matrix A = (a_{ij}) of the size

n × m. Let α_{ij} = ( \begin{matrix}a_{jj} \\ a_{ij} \end{matrix} ), i > j, E_{α_{i,j}} be a Diophantine matrix for α_{ij}. Uni-

modular matrix E_{α_{i,j}}^{n} is a matrix that is obtained from the unit matrix I_{n} of the order n after replacement of its four elements (jj), (ji), (ij), (ii) by the four elements (18), (19), (21), (22) of matrix E_{α_{i,j}}, correspondingly.

Let E_{j} = E_{α_{n,j}} ⋯ E_{α_{j+1,j}}, U = E_{n−1} ⋯ E_{1}. Then U is a unimodular matrix and UA is an upper triangular matrix. Such computation leads to the elimination of non-zero elements under the diagonal element in the first column, then in the second column, and so on up to the n − 1 column.

Denote by C_{D} an asymptotic estimation of the complexity of a Diophantine matrix calculation, M(n) – an asymptotic estimation of the complexity of multiplying matrices of the order n, M(n) = αn^{β}. The best known for today value β is less than 2.356 [3].

Lemma 1. Let R be a quasiEuclidean ring A, C ∈ R^{n×n}, A is an upper square

matrix. Then the matrix ( \begin{matrix}A \\ C \end{matrix} ) may be decomposed ( \begin{matrix}A \\ C \end{matrix} ) = Q ( Ā ; 0) with an

asymptotic estimation of complexity which is not superior than

C_{s}(n) = n^{2}C_{D} + αλ_{β}n^{β},  (1)

here λ_{β} = 7(2^{β−2} − 1)^{−1}.

Proof. The decomposition of the 2 × 1 matrix needs the multiplying by its Diophantine matrix. Without loss of generality, we suppose, that n = 2^{p} and decompose the matrices A and C into four square blocks:

A = ( \begin{matrix}a & b \\ 0 & d \end{matrix} ), C = ( \begin{matrix}e & f \\ g & h \end{matrix} ).

Let us have four decompositions.

Q_{1} ( \begin{matrix}a \\ e \end{matrix} ) = ( \begin{matrix}a_{1} \\ 0 \end{matrix} ) , denote ( \begin{matrix}b_{1} \\ f_{1} \end{matrix} ) = Q_{1} ( \begin{matrix}b \\ f \end{matrix} ) .

Q_{2} ( \begin{matrix}a_{1} \\ g \end{matrix} ) = ( \begin{matrix}a_{2} \\ 0 \end{matrix} ) , denote ( \begin{matrix}b_{2} \\ h_{1} \end{matrix} ) = Q_{2} ( \begin{matrix}b_{1} \\ h \end{matrix} ) .

Q_{3} ( \begin{matrix}d \\ f_{1} \end{matrix} ) = ( \begin{matrix}d_{1} \\ 0 \end{matrix} ) .

Q_{4} ( \begin{matrix}d_{1} \\ h_{1} \end{matrix} ) = ( \begin{matrix}d_{2} \\ 0 \end{matrix} ) .

Matrices Q_{i} = ( \begin{matrix}α_{i} & β_{i} \\ γ_{i} & δ_{i} \end{matrix} ), i = 1, …, 4, are unimodular, a, a_{1}, a_{2}, d, d_{1}, d_{2} are

upper triangular matrices. Then we get

Ā = ( \begin{matrix}a_{2} & b_{2} \\ 0 & d_{2} \end{matrix} ) and Q = ( \begin{matrix}α_{1}α_{2} & 0 & β_{1}α_{2} & β_{2} \\ α_{4}β_{3}γ_{1} + β_{4}γ_{2}α_{1} & α_{4}α_{3} & α_{4}β_{3}δ_{1} + β_{4}γ_{2}β_{1} & β_{4}δ_{2} \\ δ_{3}γ_{1} & γ_{3} & δ_{3}δ_{1} & 0 \\ γ_{4}β_{3}γ_{1} + δ_{4}γ_{2}α_{1} & γ_{4}α_{3} & γ_{4}β_{3}δ_{1} + δ_{4}γ_{2}β_{1} & δ_{4}δ_{2} \end{matrix} ) .

We obtain matrix Q as a product

( \begin{matrix}1 & 0 & 0 & 0 \\ 0 & α_{4} & 0 & β_{4} \\ 0 & 0 & 1 & 0 \\ 0 & γ_{4} & 0 & δ_{4} \end{matrix} ) ( \begin{matrix}1 & 0 & 0 & 0 \\ 0 & α_{3} & β_{2} & 0 \\ 0 & γ_{2} & δ_{2} & 0 \\ 0 & 0 & 0 & 1 \end{matrix} ) ( \begin{matrix}α_{2} & 0 & 0 & β_{2} \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ γ_{3} & 0 & 0 & δ_{3} \end{matrix} ) ( \begin{matrix}α_{1} & 0 & β_{1} & 0 \\ 0 & 1 & 0 & 0 \\ γ_{1} & 0 & δ_{1} & 0 \\ 0 & 0 & 0 & 1 \end{matrix} ) .

Therefore, these calculations need 28 multiplications of blocks of the order n/2 and four decompositions of blocks of the size n × n/2.

We reduce step by step the calculations to the blocks of size 2 × 1. Finally we obtain the following asymptotic estimation of decomposition

C_{s}(n) = n^{2}C_{D} + 7 ∑_{i=1}^{p} 4^{i}M(\frac{n}{2^{i}}) < n^{2}C_{D} + αλ_{β}n^{β}.

Here the number λ_{β} = 7(2^{β−2} − 1)^{−1} lays between λ_{3} = 7 and λ_{log 6} = 14.

Corollary 1 . Let R be a quasiEuclidean ring, A be an upper square matrix,

A ∈ R^{m×m}, C ∈ R^{(n−m)×m}. Then the matrix ( \begin{matrix}A \\ C \end{matrix} ) may be decomposed with

an asymptotic estimation of complexity which is not superior

C_{s}^{1}(n, m) = (n − m)mC_{D} + αm^{β−1}(λ_{β}(n − m) + 4n − 8m),

for n ≥ 2m and

C_{s}^{2}(n, m) = (n − m)mC_{D} + α(n − m)^{β−2}(λ_{β}m(n − m) + 10mn − 4n^{2} − 4m^{2}),

for m ≤ n ≤ 2m.

Proof.

a. Let n ≥ 2m and m be a divisor of n. Let C_{i}, i = 1, …, (n − m)/m, be an m order block of the matrix C. We count the matrix C blocks from top to bottom. And then we decompose consequently the matrices:

( \begin{matrix}A_{i+1} \\ 0 \end{matrix} ) = Q_{i} ( \begin{matrix}A_{i} \\ C_{i} \end{matrix} ), i = 1, …, (n − 2m)/m, where A_{1} = A.

We need four multiplications of matrices of the size m × m at each step. So the total number of operations is

C_{s}^{1}(n, m) = C_{s}(m)\frac{(n − m)}{m} + 4M(m)\frac{(n − 2m)}{m}.

b. Let m ≤ n ≤ 2m and n − m be a divisor of m. Let us break matrices A and C into blocks of the size (n − m) × (n − m). The matrix C contains one block row. We denote C = (C_{1}, C̄_{1}), here C_{1} is a left block of the matrix C. The matrix A contains k = m/(n − m) block rows and columns.

At first we decompose the matrix that is formed by the upper left diagonal block A_{1} of the matrix A and the left block C_{1} of the matrix C:

( \begin{matrix}A_{1} \\ C_{1} \end{matrix} ) = Q_{1} ( Ā_{1} ; 0 ), Q_{1} is a unimodular matrix.

We compute the product ( \begin{matrix}A′ & 1 \\ C′ & 2 \end{matrix} ) = Q_{1} ( Ā_{1} ; C̄_{1} ), here Ā_{1} is a part of the first

block row of the matrix A situated on the left from the diagonal block A_{1}. Denote C′_{2}=(C_{2}, C̄_{2}), here C_{2} is a left square block. We need to do 4(k − 1) block multiplications.

Then we sequentially decompose matrices that is formed by the i-th diagonal block A_{i} of the matrix A and the left block C_{i} of the matrix C′_{i}:

( \begin{matrix}A_{i} \\ C_{i} \end{matrix} ) = Q_{i} ( Ā_{i} ; 0 ), and compute the product ( \begin{matrix}A′ & i \\ C_{i+1} \end{matrix} ) = Q_{i} ( Ā_{i} ; C̄_{i} ), i = 2, …, k.

Here Q_{i} is a unimodular matrix, Ā_{i} is the part of the i-th block row of the matrix A situated on the left from the diagonal block A_{i}, C′_{i}=(C_{i}, C̄_{i}).

Substitute the blocks (1,1), (1,2), (2,1), (2,2) of the matrix Q_{i} instead of the blocks (i,i), (i,k+1), (k+1,i), (k+1,k+1) of the unit matrix I_{n} correspondingly.

As the result we get the matrix Q′_{i}. To obtain the unimodular matrix Q there is required to multiply the unimodular matrices Q′_{i}:

Q = Q′_{k} ⋯ Q′_{2}Q′_{1}.

One such matrix multiplication needs four block multiplications and the i-th step needs 4(k − i) + 4 block multiplications, i = 2, …, k.

So the total number of operations is (∑_{i=1}^{k} 4(k − i + 1) − 4)M(n − m)+ kC_{s}(n − m)+ < kC_{s}(n − m) + 2(k^{2} + k − 2)M(n − m) = m(n − m)C_{D} + M(n − m)(λ_{β}m/(n − m) +2(5mn − 2n^{2} − 2m^{2})/(n − m)^{2}) = C_{s}^{2}(n, m).

Remark. It is easy to show that we can change the upper triangular matrix A to a lower triangular matrix in the condition of Lemma 1. We also can change the order of multipliers in the decomposition formula, the evaluations in the lemma and in the corollary will be the same.

Theorem 1. In a quasiEuclidean ring any quadratic matrix A of the order n may be decomposed with an asymptotic estimation of complexity which is not superior

C_{t}(n) = μ_{β}M(n) + \frac{n^{2} − n}{2}C_{D},  (2)

where μ_{β} = \frac{9·2^{β−2}+19}{(2^{β}−2)(2^{β−2}−1)}.

Proof. The decomposition of a second order matrix needs the multiplying of its first column by the Diophantine matrix.

We suppose that n = 2^{p} and decompose the matrix A into four square blocks:

A = ( \begin{matrix}A & B \\ C & D \end{matrix} ) .

1. Let the block A be decomposable:

A = Ã^{−1}Ā.

Denote A_{1} = ÃA, B_{1} = ÃB.

2. Let the block ( \begin{matrix}A_{1} \\ C \end{matrix} ) be decomposable:

( \begin{matrix}A_{1} \\ C \end{matrix} ) = Q^{−1} ( \begin{matrix}A_{2} \\ 0 \end{matrix} ) ,  (3)

Q – unimodular, A_{2} – upper triangular. Denote ( \begin{matrix}B_{2} \\ D_{1} \end{matrix} ) = Q ( \begin{matrix}B_{1} \\ D \end{matrix} ).

3. Let the block D_{1} be decomposable: D_{1} = D̃_{1}^{−1}D̄_{1}. Denote D_{2} = D̄_{1}.

Then the matrix A is decomposable:

A = Ã^{−1}Ā, Ã = diag(I, D̃_{1}) · Q · diag(Ã, I), Ā = ( \begin{matrix}A_{2} & B_{2} \\ 0 & D_{2} \end{matrix} ) .

Such decomposition requires 9 block multiplications, two block decompositions and one decomposition of the rectangular n × n/2 block.

We denote by C_{t}(n) the complexity of decomposition of a matrix of the order n. Then we obtain according to (1)

C_{t}(n) = 2C_{t}(n/2) + 9M(n/2) + C_{s}(n).

According to the lemma 1 we receive

C_{t}(n) = 2C_{t}(n/2) + (n/2)^{2}C_{D} + α(9 + λ_{β})(n/2)^{β}.

We shall break matrices A and D into the four square blocks, and continue such procedure up to the second order blocks.

The asymptotic estimation of the complexity is

∑_{i=1}^{p} 2^{i−1}((\frac{n}{2^{i}} )^{2}C_{D} + α(9 + λ_{β})(n/2^{i})^{β}) < \frac{n^{2} − n}{2}C_{D} + αμ_{β}n^{β} = C_{t}(n)

where the number μ_{β} = \frac{1}{2^{β}−2}(9 + \frac{7}{2^{β−2}−1}) lays between μ_{3} = 8/3 and μ_{log 6} = 23/4.

Corollary 2 In a quasiEuclidean ring any rectangular matrix A of the size n×m is decomposable. Such decomposition is possible to execute with an asymptotic estimation of complexity which is not higher, than C_{t}^{0}(n, m) = C_{D}(n^{2} − n)/2 + αn^{β}(μ_{β} + (m − n)/n), for n < m, C_{t}^{1}(n, m) = C_{D}(nm − (m^{2} + m)/2) + αμ_{β}m^{β}+ α(n − m)^{β−1}(λ_{β}m(n − m) + 10mn − 4n^{2} − 4m^{2}), for m ≤ n ≤ 2m and C_{t}^{2}(n, m) = C_{D}(nm − (m^{2} + m)/2)+ αm^{β−1}(μ_{β}m + λ_{β}(n − m) + 4n − 8m), for n ≥ 2m.

Proof. The proof follows from the estimations for C_{t}(n) and C_{s}(n) in (1) and (2).

For n < m it is enough to decompose the left square block and then to multiply the left block from m − n columns of a matrix A by Ã.

For n ≥ m it is enough to decompose the upper square block and to use the corollary of Lemma 1.

Theorem 2. In a quasiEuclidean ring any system of linear equations Ay = c with the upper triangular coefficient matrix A ∈ R^{n×n} may be solved with an asymptotic estimation of complexity

C_{n} = C_{D}(n^{2} − n)/2 + αϕ_{β}n^{β}, ϕ_{β} = \frac{1}{2^{β} − 2}(5 + \frac{7}{2^{β−2} − 1}).  (4)

Proof. Solving the system Ay = c is reduced to solving two systems of twice less size and one matrix decomposing.

Let ( \begin{matrix}A & B \\ 0 & D \end{matrix} ) ( \begin{matrix}y_{1} \\ y_{0} \end{matrix} ) = ( \begin{matrix}f \\ h \end{matrix} ) be the block form of the system Ay = c. All the

matrix A blocks are square and A, D are the upper triangular blocks.

(1). Let y_{0} = V s_{0} + v be a solution of the system Dy_{0} = h, V ∈ R^{n/2×n/2}, v, s_{0} ∈ R^{n/2}, s_{0} is a vector of parameters. Substituting it into the system Ay_{1} + By_{0} = f we obtain Ay_{1} + Gs_{0} = g, here G = BV, g = f − Bv.

(2). Decompose the matrix (A, G): (A, G)P = (Ā, 0). P is a unimodular matrix, Ā is an upper triangular matrix. The complexity of such decomposition is C_{s}(n/2) according to the lemma 1 and the remark.

(3). Let z_{1} = W t_{1} + w be a solution of the system Āz_{1} = g, W ∈ R^{n/2×n/2}, w, t_{1} ∈ R^{n/2}, t_{1} is a vector of parameters.

Then

( \begin{matrix}z_{1} \\ z_{0} \end{matrix} ) = diag(W, I) ( \begin{matrix}t_{1} \\ t_{0} \end{matrix} ) + ( \begin{matrix}w \\ 0 \end{matrix} )

is the solution of the system

(Ā, 0) ( \begin{matrix}z_{1} \\ z_{0} \end{matrix} ) = g.

Let P = ( \begin{matrix}P_{a} & P_{b} \\ P_{c} & P_{d} \end{matrix} ) be the block form of the matrix P. Since ( \begin{matrix}y_{1} \\ s_{0} \end{matrix} ) = P ( \begin{matrix}z_{1} \\ z_{0} \end{matrix} )

and y_{0} = V s_{0} + v, the solution of the initial system Ay = c is obtained in the form

y = ( \begin{matrix}P_{a}W & P_{b} \\ V P_{c}W & V P_{d} \end{matrix} ) ( \begin{matrix}t_{1} \\ t_{0} \end{matrix} ) + ( \begin{matrix}P_{a}w \\ V P_{c}w + v \end{matrix} ) .

In the whole we need five block multiplications. Denote by C_{n} the complexity of solving a system with square coefficient matrix of the order n. Then

C_{n} = 2C_{n/2} + 5M(n/2) + C_{s}(n/2) = 2C_{n/2} + α(5 + λ_{β})n^{n/2} + (n/2)^{2}C_{D}.

Hence we obtain the asymptotic evaluation

C_{n} = C_{D}(n^{2} − n)/2 + αϕ_{β}n^{β}, where ϕ_{β} = \frac{1}{2^{β} − 2}(5 + \frac{7}{2^{β−2} − 1}).

Theorem 3. In a quasiEuclidean ring any system of linear equations Ax = c, A ∈ R^{n×m} may be solved with an asymptotic estimation of the complexity C_{nm}:

C_{nm} = (nm − l)C_{D} + αψ_{β}l^{β},

ψ_{β} = μ_{β} + ϕ_{β} + λ_{β} \frac{n−m}{m} + \frac{4n}{m} − 8 for 2m ≤ n, ψ_{β} = μ_{β} + ϕ_{β} + λ_{β} \frac{m−n}{n} + \frac{6m}{n} − 9 for 2n ≤ m;

C_{nm} = (nm − l)C_{D} + αψ_{β}l^{β} + α|n − m|^{β−2}(λ_{β}l|m − n| + 10nm − 4m^{2} − 4n^{2}),

ψ_{β} = μ_{β} + ϕ_{β} for m ≤ n ≤ 2m, ψ_{β} = μ_{β} + ϕ_{β} + \frac{2m}{n} − 1 for n < m ≤ 2n, here l = min(m, n).

Proof. The decomposition (A, c) = P (A′, c′) with the upper triangular matrix (A′, c′) may be obtained with the complexity C_{t}^{∗}(n, m) according to the corollary 2.

Therefore it is enough to solve the system

A′x = c′  (5)

with upper triangular coefficient matrix.

a. Let m > n, A′ = (A, B), A be a square block with an upper triangular form.

According to the corollary 1 and the remark we can obtain the unimodular matrix Q such that (A, B)Q = (Ā, 0), with the complexity C_{s}^{∗}(n, m). Here Ā is an upper triangular matrix.

The solutions of the system (5) and of the system Āy = c′ are connected

according to the equation: x = Q ( \begin{matrix}y \\ t_{0} \end{matrix} ), where t_{0} is a vector of parameters.

Let Ut_{1} + u be a solution of the system Āy = c′, U ∈ R^{n×n}, t_{1}, u, y ∈ R^{n}, t_{1} be a vector of parameters, t_{1} ∈ R^{n}.

Then the solution of initial system will be x = Qdiag(U, I)t + Q ( \begin{matrix}u \\ 0 \end{matrix} ), where

t ∈ R^{m} is a vector of parameters. The computing of matrix product Qdiag(U, I) needs αn^{β−1}m operations. According to the theorem 2 the solution of the system Āy = c′ may be obtained with the complexity C_{n}.

In the whole we obtain the evaluation of the complexity of the system Ax = c solving

C_{nm} = C_{t}^{0}(n, m) + C_{s}^{i}(m, n) + C_{n} + M(n)(m/n),

i = 1, 2, where i = 1 in the case m ≥ 2n and i = 2 in the case 2n ≥ m > n.

b. Let n ≥ m, (A′, c′) = ( \begin{matrix}A & c_{A} \\ 0 & c_{B} \end{matrix} ), A be a square block with upper triangular

form. The system is incompatible if the vector c_{B} is non-zero.

According to the theorem 2 we can obtain the solution of the system Ax = c_{A} with the complexity C_{m}.

So the total complexity is

C_{nm} = C_{t}^{i}(n, m) + C_{m} i = 1, 2,  (6)

where i = 1 in the case 2m ≥ n ≥ m and i = 2 in the case n ≥ 2m.

Therefore the theorem assertion follows.

In particular for the system with square coefficient matrix we obtain

C_{nn} = n^{2}C_{D} + αψ_{β}n^{β}, ψ_{β} = \frac{7 · 2^{β−1}}{(2^{β} − 2)(2^{β−1} − 2)}.

We have ψ_{3} < 4, 67, ψ_{log7} < 6, 54, ψ_{log6} = 10, 5.

3 Hensel lift

Let R be a commutative domain with identity, K be a quotient field of R, A ∈ R^{n×m}, c ∈ R^{n}. Let

Ax = c  (7)

be a system of linear equations in R.

Suppose, that n < m, rankA = n, k = m − n + 1 is a number of the linearly independent solutions of the system (7). Let matrix A = (A_{0}, A_{1}) have a square block A_{0}, det A_{0} ≠ 0, and the last row of a the matrix A_{0}^{−1}(A_{1}, c) not contain nulls. The system in general may be reduced to such case (it may be done for example as in [10]).

Let e_{j} be the j column of the unit matrix I_{k} = (e_{1}, …, e_{k}) of the order k. We

seek the solution of the system (7) in the form ( \begin{matrix}x_{j} \\ ξ_{j}e_{j} \end{matrix} ), x_{j} ∈ K^{n−1}, ξ_{j} ∈ K\ 0,

(j = 1, …, k). Then we obtain a determined system of the order n.

It is known, that such system may be solved by p-adic lift [4]. If R is a ring of integers, the solution may be obtained with

D_{Z} = O(n^{3}(log n + log ∥A∥ + log p)^{2})

bits operations. If R is a ring of polynomials over a field F, then solving the system needs

D_{F [x]} = O(n^{3}(deg A + deg p)^{2})

operations in the field F [11], [4].

To solve k such systems we need

kD_{Z}, and kD_{F [x]}  (8)

operations correspondingly.

As ξ_{j} ≠ 0, j = 1, …, k, the solutions of these k systems are linearly independent. Reduce the components of each solution to their common denominator (it is a divider of some minor of the matrix A), and write the solutions in the

following way: x_{j} = ( \begin{matrix}x_{j} \\ α_{j}e_{j} \end{matrix} ) χ_{j}^{−1}, j = 1, …, k, x_{j} ∈ R^{n−1}, α_{j}, χ_{j} ∈ R\0.

The probabilistic approach to finding Diophantine solutions consists in the following.

Denote ( \begin{matrix}x_{j} \\ α_{j}e_{j} \end{matrix} ) by w_{j}. It is easy to show that if \{w_{j}χ_{j}^{−1}|i = 1, 2, …, k\} is

a basis set of solutions of the system (7), then

\{\frac{∑_{j=1}^{k} w_{j}η_{j}}{∑_{j=1}^{k} χ_{j}η_{j}} | η_{j} ∈ R, ∑_{j=1}^{k} χ_{j}η_{j} ≠ 0\}

is the set of all solutions of this system [10].

So the system has Diophantine solutions if the ideal generated by the denominators of a basis set of solutions is unit. In this case we can find such η_{j} ∈ R that ∑_{j=1}^{k} χ_{j}η_{j} = 1 and ∑_{j=1}^{k} w_{j}η_{j} is a Diophantine solution. Moreover if x_{i} = w_{i}χ_{i}^{−1}, i = 1, 2, …, h is a basis set of solutions and χ_{1} = 1, then the set w_{1}, w_{i} − w_{1}(χ_{i} − 1), i = 2, 3, …, h, is an integer basis [10].

If the ideal generated by the denominators of a basis set of solutions is not unit then we can rearrange columns in the matrix A and repeat the whole computation process. It might be necessary to repeat such procedure several times.

An expected number of such computations for the ring Z is O(k^{−1}(log n + log log ∥(A, c)∥)) [10], [11].

The deterministic approach to finding Diophantine solutions consists in the following.

The number LCM(α_{1}, …, α_{k}) is a divider of the minor of a matrix (A, c), constructed from the first n − 1 columns and the last column, according to Crammer’s rule. Denote δ = LCM(α_{1}, …, α_{k}).

Multiply the numerator and the denominator of the solution x_{j} by δα_{j}^{−1} and denote x_{j}δα_{j}^{−1} = v_{j}, χ_{j}δα_{j}^{−1} = δ_{j} for each j. We obtain the basic set of the

solutions of the system (7) in the form x_{j} = ( \begin{matrix}v_{j} \\ δe_{j} \end{matrix} ) δ_{j}^{−1}.

Denote the vector (δ_{1}, …, δ_{k}) by u, u ∈ R^{k}, and the matrix (v_{1}, …, v_{k}) by V, V ∈ R^{(n−1)×k}.

Theorem 4. The system (7) and

\{ \begin{matrix}V y & = & 0 mod δ \\ uy & = & δ \end{matrix}  (9)

are consistent simultaneously in the ring R.

If Y is the set of the solutions of the system (9), then

X = \{ ( \begin{matrix}y \\ V yδ^{−1} \end{matrix} ) | y ∈ Y \} is the set of the solutions of the system (7).

Proof. The set of the solutions of the linear system (7) is a plane of the dimension m−n in the space K^{m}. It is defined by m−n+1 linearly independent solutions. Therefore, for any solution z ∈ K^{m} of the system (7) the numbers σ_{i} ∈ K, i = 1, …, k, exist, such that

∑_{i=1}^{k} x_{i}σ_{i} = z and ∑_{i=1}^{k} σ_{i} = 1.

Denote σ_{i}δ_{i}^{−1} = τ_{i} ∈ K and z = (y, w), y ∈ R^{k}, w ∈ R^{m−k}, t = (τ_{1}, …, τ_{k}). Then

∑_{i=1}^{k} δ_{i}τ_{i}x_{i} = z and ut = 1,

and these equations may be written as the system:

( \begin{matrix}δI_{k} & −I_{k} & 0 \\ V & 0 & −I_{n−1} \\ u & 0 & 0 \end{matrix} ) ( \begin{matrix}t \\ y \\ w \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} )

Multiplying each part of this system by the matrix

( \begin{matrix}I_{k} & 0 & 0 \\ −V & δI_{n−1} & 0 \\ −u & 0 & δ \end{matrix} ) ,

invertible in K, we obtain a system

( \begin{matrix}δI_{k} & −I_{k} & 0 \\ 0 & V & −δI_{n−1} \\ 0 & u & 0 \end{matrix} ) ( \begin{matrix}t \\ y \\ w \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ δ \end{matrix} ) .  (10)

Let the system (7) have the Diophantine solution z = (y, w) ∈ R^{m}. Then the system (10) is consistent. Therefore the system (9) is also consistent.

To the contrary, if the system (9) is consistent and y ∈ R^{k} is its solution, then (δ^{−1}y, y, δ^{−1}V y) is the solution of a system (10). Therefore z = (y, δ^{−1}V y) is the Diophantine solution of the system (7).

Now we estimate the complexity of the solution of the system (9) in the ring of integers and in the ring of polynomials. The matrix V has the size (n − 1) × k, and the module δ is a divider of a minor of the order n. Therefore, for R = Z

log δ ≤ n(log ∥A∥ + 0.5 log n),  (11)

where ∥A∥ is an absolute value of the greatest element of the matrix A. For R = F [z]

deg δ ≤ n deg ∥A∥,  (12)

where deg ∥A∥ is the greatest degree of elements of the matrix A.

Let us use the standard algorithm of multiplying integers and polynomials. Then the multiplying of N bit numbers requires N^{2} bit operations, the multiplying of two polynomials of a degree N requires N^{2} multiplications in the field F.

By Theorem 3 we obtain the following asymptotic estimations of operations for Z and F [x], correspondingly:

η_{β}knM(n)(log ∥A∥ + 0.5 log n)^{2}, η_{β}knM(n) deg^{2} ∥A∥,  (13)

for k ≥ n, and for k < n we obtain

δk^{β−1}(log ∥A∥ + (log n)/2)^{2}(n^{3}(4 + λ_{β}) + n^{2}k(ν_{β} + μ_{β} − λ_{β} − 4)),  (14)

δk^{β−1} deg^{2} ∥A∥(n^{3}(4 + λ_{β}) + n^{2}k(ν_{β} + μ_{β} − λ_{β} − 4)).  (15)

In general, these estimations do not depend on the complexity (8) of k solutions of the system (7).

However, for quasi-square systems it is not so, the contribution of the value of (8) is essential. In this case for Z and F [x] we have the following estimations

O(n^{3}(log n + log ∥A∥ + log p)^{2}), O(n^{3}(deg ∥A∥ + deg p)^{2}).  (16)

Fast methods of multiplying integers and polynomials correspondingly improve these estimations. For example, the Karasuba method demands O(N^{log 3}) operations and the FFT method demands O(N log N log log N) operations.

4 Modular Method

Let us solve the system (7), supposing that A = (A_{0}, A_{1}), block A_{0} is a square block, det A_{0} = δ ≠ 0. A recursive method [9] allows to transform the system (7) to the diagonal form

(δI_{n}, A′_{1})x = c′  (17)

with O(n^{β−1}m) operations in the ring R. The methods described in [7] and [8] can do this transformation with O(n^{2}m) operations.

All coefficients of the system (17) are minors of the order n of matrix (A, c). Therefore estimations (11) and (12) are suitable to the ring of integers and ring of polynomials, correspondingly.

Calculation needs

O(mn^{β+1}(log n + log ∥A∥)^{2}) and O(mn^{β+1} deg^{2} ∥A∥),

operations in the rings Z and F [x], correspondingly.

The Chinese remainder theorem allows to reduce amount of operations up to

O(mn^{β}(log n + log ∥A∥)^{2}) and O(mn^{β} deg^{2} |A |).  (18)

The solving of the system (17) is reduced to the solving of a system

A′_{1}y = c′ mod δ.  (19)

in the quotient-rings R/δR.

Let Y be a set of the solutions of the system (19), then

X = \{ ( \begin{matrix}(c′ − A′_{1}y)δ^{−1} \\ y \end{matrix} ) | y ∈ Y \}

is the set of the solutions of the system (17).

Let us estimate the complexity of the solution of the system (19). The matrix A′ has the size n × (m − n), and the module δ is a minor of the order n.

By Theorem 3 we obtain the following estimations for Z and F [x], correspondingly:

O(n^{β+1}(m − n)(log ∥A∥ + log n)^{2}), O(n^{β+1}(m − n) deg^{2} ∥A∥). (20)

for m − n ≥ n, and for m − n ≤ n

O(n^{2}m(m − n)^{β−1}(log ∥A∥ + log n)^{2}), O(n^{2}m(m − n)^{β−1} deg^{2} ∥A∥). (21)

In general (20) and (21) are the estimations for the modular method. However, for quasi-square systems the value m−n can be neglected and the estimation (18) will give the basic contribution. In this case we obtain the following estimations for the number of operations:

O(n^{β+1}(log n + log ∥A∥)^{2}) and O(n^{β+1} deg^{2} ∥A∥)  (22)

for Z and F [x], correspondingly.

5 Resume

The evaluations (13)–(16) and (20)–(22) obtained in the last two sections allow to formulate the next theorem about the complexity (number of bit operations) of solving a linear Diophantine system of the size n × m.

Theorem 5. For essentially rectangular systems, when m > 2n, both methods have the identical complexity O(mn^{β+1}). For quasi-square systems the complexity of the p-adic method is O(n^{3}) and the complexity of the modular method is O(n^{β+1}).

All estimations in the theorem are given within a logarithmic factors. In the case of the polynomial ring F [x] the analogous theorem is true for the evaluation of the number of operations in the field F.

The comparison of the p-adic and the modular methods shows, that for quasi-square systems the p-adic method has the advantage, it is O(n^{β−2}) times as fast.

The other advantage of the p-adic method for systems of the arbitrary size is the possibility to use the probability approach. The general solution of a system can be obtained with a high probability with O((m − n + 1)n^{3}) operations [10]. An expected amount of attempts will be O((log n + log log ∥A∥)/(m − n + 1)).

These estimations can still be improved, if the fast algorithms of multiplying would be used.

References

1. Abbott, J., Bronstein, M. and Mulders, T.: Fast Deterministic Computation of Determinants of Dense Matrices. In: Proceedings of ISSAC’99. ACM Press, Vancouver (1999) 197–204

2. Bareiss, E. H.: Computational solutions of matrix problems over an integral domain. J. Inst. Maths Applics 10 (1972) 68–104

3. Coppersmith, D. and Winograd, S.: Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation 9 (1990) 251–280

4. Dixon, J.: Exact solution of linear equations using p-adic expansions. Numer. Math. 40 (1982) 137-141

5. Gregory, R. T., Krishnamurthy E. V.: Methods and Applications of Error-Free Computation. Springer-Verlag, Berlin (1984)

6. Hurt, M. F., Waid, C. A.: A generalized inverse which give all the integral solutions to a system of linear equations. SIAM J. Appl. Math. 10, (1970) 547–550 7. Malaschonok, G. I.: Solution of a system of linear equations in an integral do-

main. USSR Journal of Computational Mathematics and Mathematical Physics 23 (1983) 1497–1500

8. Malaschonok, G. I.: Algorithms for the solution of systems of linear equations in commutative rings. In: Mora, T. and Traverso, C. (eds.): Effective Methods in Algebraic Geometry. Progress in Mathematics 94, Birkhauser (1991) 289–298

9. Malaschonok, G. I.: Recursive Method for the Solution of Systems of Linear Equations. In: Sydow, A. (ed.): Computational Mathematics (Proceedings of the 15th IMACS World Congress, Vol. I, Berlin, August 1997). Wissenschaft & Technik Verlag, Berlin (1997) 475–480

10. Malaschonok, G. I.: Effective Matrix Methods in Commutative Domains, In: Krob, D., Mikhalev, A. A., Mikhalev, A. V. (eds.): Formal Power Series and Algebraic Combinatorics. Springer (2000) 506–517

11. Mulders, T., Storjohann, A.: Diophantine Linear System Solving. In: Proceedings of ISSAC’99. ACM Press, Vancouver (1999) 181–188