A simple randomized algorithm is given for finding an integer solution to a system of linear Diophantine equations. Given as input a system which admits an integer solution, the algorithm can be used to find such a solution with probability at least 1/2. The running time (number of bit operations) is essentially cubic in the dimension of the system. The analogous result is presented for linear systems over the ring of polynomials with coefficients from a field.
Solving a system of linear Diophantine equations is a classical mathematical problem: given an integer matrix -4 and vector b, the goal is to find an integer vector x that satisfies Ax = b. We present a simple randomized algorithm for solving this problem. We also show how to adapt, the algorithm to solve linear systems over the ring of univariate polynomials with coefficients from a field. The algorithm is easy to implement, memory efficient, and faster than previous methods. The algorithm is probabilistic in the following sense. If a solution vector is returned it is guaranteed to be correct;. On the other hand, the algorithm does not currently certify inconsistency and may return NIL for a system which does admit a solution; the chance of this happening is bounded by a user specified error tolerance f > 0.
For the moment we restrict the discussion to integer systems; the case for polynomial systems is analogous. In this paper we analyse our algorithms under the assumption of standard, quadratic integer arithmetic. Let 0(M(t)) bit operations be sufficient to multiply two t-bit integers. Then standard arithmetic has M(t) = t~. FFT-based methods allow M(t) = t log «log log t.
Let n be the larger of the row or column dimension of a linear Diophantine system Ax = b. Let ft be a bound on the magnitudes of entries in A and b. If Ax = b admits an integer solution for x, then it admits an integer solution with entries having length 0(n(log?(. + log ft)) bits — this bound is tight in the worst case. The algorithm we present here uses an expected iiumbo:r of 0~ (n* (log ft)2) bit operations to return
1'ermission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not. made or distributed for profit or commercial advantage, and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. ISSAC '99, Vancouver. British Columbia, Canada. © 1Э99 ACM 1-581 13-073-2 / :i'J I 07 * 5.00
a solution with this size bound. For input systems with logd small compared to n, this improves on the previously-fastest algorithm which uses О' (пл log/.? + nM(n\og/3)) bit operations in the worst case . Moreover, the algorithm we present here, like that in . requires additional storage for only О [n2 (log n + log ft)) bits. This space complexity, which is linear in the size of the output, vector and essentially linear in the size of a dense input, system, is an important, feature of the algorithm.
The global technique of our Diophantine solver is similar to the randomized algorithm in . The key idea is to compute a small number of rational solutions of perturbations of the input system. In , the perturbed systems have the form UALy = Ub. where U and L are randomly chosen Toeplitz matrices. In our algorithm, the perturbed systems have the form APy = b. where P is a random dense matrix. Note that if у is a rational solution to APy = b, then Py is a rational solution to the original system Ax — b. The idea is to compute a sequence of these perturbed rational solutions which can then hopefully be combined to obtain a Diophantine solution. This suggests an iterative algorithm which is supposed to "converge" to a Diophantine solution. A worked example demonstrating this technique is given in Section 2. The challenge is twofold. First, show how to efficiently solve the perturbed rational systems. Second, specify how to choose the entries in the perturbation matrices, and prove that fast convergence can be expected with these choices.
In . the perturbed rational systems are solved using the Wiedemann algorithm together with a homomor-phic imaging scheme: this leads to an algorithm which admits an easy coarse grain parallelization and exploits the possible structure or sparsity of A. Choosing the U and L to be Toeplitz is important so that matrix-vector products involving U and L do not cost more than matrix-vector products involving .4. For a very sparse matrix A, with number of nonzero entries essentially linear instead of quadratic in a, this approach leads to an algorithm with running time 0~ (n1 log ft + nM(nlogft)) bit operations. For a dense input system, though, the running time is 0~(пл log/? + nM(n log ft)) bit operations. In order to prove that a small (i.e. logarithmic) number of perturbed rational solutions will converge1 to a Diophantine solution, entries in U and L are chosen from a small extension ring of the integers.
Our approach is to solve the perturbed rational systems using linear p-adic lifting as described in . In Section 5 we offer pseudo-code for thep-adic rational solver together with