solving linear diophantine equation with inequalities

Question

Is there an easy way to solve linear Diophantine equations with inequalities?

For example, say I have $a_1x_1 + a_2x_2 \equiv k \mod m$ where:

• $a_1, a_2, k, m$ are given

• I already know $b_1, b_2$ such that $a_1b_1 \equiv a_2b_2 \equiv 1 \mod m$, so if either of the $x_1$ or $x_2$ is given, I can calculate the other one:

  $x_1 = (k - a_2x_2)b_1 \bmod m \\ x_2 = (k - a_1x_1)b_2 \bmod m$

• I wish to solve for $x_1, x_2$ where $0 < x_1 < x_{max}$, $0 < x_2 < x_{max}$.

---

For example, $a_1 = 12345678987601, a_2 = 45678901, k = 135990606993339, m = 2^{48}-1, x_{max}=10^6.$

I can find $b_1 = 179768904497071, b_2 = 154376496085921$ easily, but other than a brute-force search, I'm not sure how to find the values $x_1$ and $x_2$ (in this case $x_1 = 50603, x_2 = 12481$ but that's because I created this example)

---

Note: I used tmyklebu's solution and implemented in Python with numpy and sympy:

import numpy as np
from sympy.matrices import Matrix

# thanks to user tmyklebu http://math.stackexchange.com/a/1151038/120

def reduce_lattice(L, verbose=False):
    """repeated pseudo-Gram-Schmidt orthogonalization
       until we get no more progress
    """
    N = L.shape[0]
    get_norm = lambda M: sum(np.amax(np.abs(M),1))
    Lnorm = None
    Lnew = L.copy()
    Lnew_norm = get_norm(Lnew)
    while Lnorm != Lnew_norm:
        for i in xrange(N):
            v = Lnew[i]
            num = np.dot(Lnew,v)
            den = np.sum(Lnew**2, axis=1)
            k = [int(np.round(x)) for x in num*1.0/den]
            k[i] = 0   # don't subtract multiples of self
            vnew = v-np.dot(k,Lnew)
            Lnew[i] = vnew
        Lnorm = Lnew_norm
        Lnew_norm = get_norm(Lnew)
    return Lnew

def solve_eqn(a1,a2,k,m,xmax):
    B = int(np.ceil(np.abs(xmax)))
    L = np.array([[m*B,0,0],
                  [a1*B,1,0],
                  [a2*B,0,1]], dtype=object)
    L1 = reduce_lattice(L)
    L2 = Matrix(L1.T)
    exact_coords=L2.LUsolve(Matrix(3,1,[k*B,0,0]))
    int_coords = Matrix(3,1,[int(c.round()) for c in exact_coords])
    soln = [int(c) for c in L2*int_coords]
    assert soln[0] == k*B
    return soln[1:]

and got:

>>> solve_eqn(a1=12345678987601, a2=45678901, 
              k=135990606993339, m=(1<<48)-1, xmax=1000000 )
[50603, 12481]

:-)

Answer 1

I'm going to pretend you want $|x_i| < x_{max}$ for $i \in \{1, 2\}$ instead of $0 \leq x_i < x_{max}$; it's straightforward to transform your problem into this form.

Consider the lattice $L$ in $\mathbb{R}^3$ spanned by $(m, 0, 0)$, $(a_1, 1/x_{max}, 0)$, and $(a_2, 0, 1/x_{max})$. These three vectors are called a "basis" for the lattice $L$. Every vector in this lattice is of the form $(a_1 x_1 + a_2 x_2 + z m, x_1/x_{max}, x_2/x_{max})$ for some integers $x_1$, $x_2$, and $z$.

If there is a vector in $L$ whose $\infty$-norm distance to $(k, 0, 0)$ is less than 1, then there are $x_1$, $x_2$, and $z$ such that $a_1 x_1 + a_2 x_2 + m z = k$ (since the first coordinate is always an integer), $|x_1| < x_{max}$, and $|x_2| < x_{max}$. Conversely, if a solution to your problem exists, it corresponds to a vector in $L$ whose $\infty$-norm distance to $(k, 0, 0)$ is less than 1.

---

So we've reduced your problem to that of finding the closest vector in $L$ to some point $p = (k, 0, 0)$. Algorithmically, the first step is to find a "good" basis for $L$ made of "short" and "nearly orthogonal" vectors. Then you express $p$ as a linear combination of those vectors. Since we're interested in landing on a box with side length 1, you can derive bounds on the coefficients of the vectors in the reduced basis and then check the small, finite number of sets of coefficients satisfying those bounds.

Finding a "good" basis for $L$ is a very interesting problem in higher dimensions. We're in three dimensions, though, so there's not a lot that can go wrong. Repeatedly look for one vector that can be made shorter by adding or subtracting integer multiples of the other two (this is a two-dimensional closest vector problem) and replace it with the shorter version.

Firstly, other, better, algorithms for these lattice problems exist, and I'm by no means an expert on this stuff. Secondly, there is quite probably another, simpler way to skin this cat. The lattice-based approach generalises straightforwardly to higher dimensions (i.e. more variables) and I think the formulation of this problem using lattices is particularly elegant.

---

Working your example, you begin with the lattice basis (I scaled everything by a million and all my vectors are row vectors)

281474976710655000000 0 0
12345678987601000000  1 0
45678901000000        0 1

You can reduce this (with a little work) to the basis

 5000000 -123921   1584701
-2000000 -2566579  5272172
 1000000 -6801061 -5903638

In this basis, the vector (135990606993339000000 0 0) has (approximate) coordinates (24644034888697.000457924580143 -5560521170833.99591073404456 1649390208186.00588890901). If you round those coordinates to their nearest integers and see what point the rounded coordinates correspond to, it's (135990606993339000000 50603 12481). (Note that just rounding all the coordinates doesn't always give you the closest lattice vector.)