Fastest Integer Search for Any Arbitrary X and Y

Question

With the application of Discrete Math, what's the fastest algorithm in python to solve this problem:

With the equation ax + by = d, where a, b and d are user inputs, search for the integer value pairs for x and y in the range |L, R| (including L and R) that satisfies the equation.

L and R are user inputs as well.

Output all possible values for x and y in increasing order. Print none if there are no possible pairs.

Case 1:

a = 1
b = 5
d = 40
L = 0
R = 10

Result:
0 8
5 7
10 6

Case 2:

a = 14
b = 91
d = 53
L = 1
R = 100

Result:
none

Here's my code, but I believe there is a much faster way in searching. This is too inefficient.

a = int(input())
b = int(input())
d = int(input())
L = int(input())
R = int(input())

isNone = True
for x in range(L, R+1):
    for y in range(L, R+1):
        if (a*x) + (b*y) == d:
            print(x, y)
            isNone = False
if isNone: print("none")

Can there be an algorithm with O(1)? What's the fastest way?

Answer 1

The identity I imagine you want to apply is as follows (straight from wikipedia, though a similar phrasing should be found in most discrete math texts or could be proven on your own):

The simplest linear Diophantine equation takes the form ax + by = c, where a, b and c are given integers. The solutions are described by the following theorem:

This Diophantine equation has a solution (where x and y are integers) if and only if c is a multiple of the greatest common divisor of a and b. Moreover, if (x, y) is a solution, then the other solutions have the form (x + kv, y − ku), where k is an arbitrary integer, and u and v are the quotients of a and b (respectively) by the greatest common divisor of a and b.

This nearly immediately answers the question, especially since the standard proof of it uses something called the euclidean algorithm. In the interest of simplicity we're going to do the following:

1. Use the euclidean algorithm in the forward direction to find g = gcd(a, b).
2. Backsolve through the euclidean algorithm to find _x, _y such that _x*a + _y*b == g.
3. If d isn't a multiple of g then there can't possibly be any solutions, so exit early.
4. Otherwise, x, y = _x*(d//g), _y*(d//g) is a possible solution. Use that to find all solutions in the desired range.

def gcd(a, b):
    # forward euclidean algorithm
    q, r, x, qs = None, b, a, []
    while x%r:
        (q, r), x = divmod(x, r), r
        qs.append(q)

    # save the gcd for later
    g = r

    # backsolve euclidean algorithm
    if not qs:
        return 1, 1-a//b, g
    theta, omega = 1, -qs[-1]
    for q in reversed(qs[:-1]):
        theta, omega = omega, theta - q * omega
    
    # theta * a + omega * b == g
    # g might be negative, but we don't care about a canonical solution
    return theta, omega, g

def idivide_zero(a, b):
    # integer division a/b, round toward 0 instead of round down
    q = a // b
    if q < 0 and b*q != a:
        q += 1
    return q

def bounded_solutions(a, b, d, L, R):
    _x, _y, g = gcd(a, b)
    if d%g:
        return

    # a*x + b*y == d
    x, y = _x*(d//g), _y*(d//g)

    # solutions are of the form (x+k*v, y-k*u)
    u, v = a//g, b//g

    # The next trick is to find all solutions in [L, R].
    # Basically, we need L <= x+k*v <= R and L <= y-k*u <= R.
    # Note that valid choices of k exist in a contiguous interval, so
    # we only have to find the lower and upper bounds to be able to
    # quickly enumerate all options.
    xb = sorted(idivide_zero(b-x, v) for b in (L, R))
    yb = sorted(idivide_zero(y-b, u) for b in (L, R))
    m, M = min(xb[0], yb[0]), max(xb[1], yb[1])
    for k in range(m, M+1):
        yield x+k*v, y-k*u

a = int(input())
b = int(input())
d = int(input())
L = int(input())
R = int(input())

empty = True
for x, y in bounded_solutions(a, b, d, L, R):
    print(x, y)
    empty = False
if empty:
    print('none')

Code is untested. It's right in spirit, but there may be some debugging left.