Solving linear equations

Question

I have to find out the integral solution of a equation ax+by=c such that x>=0 and y>=0 and value of (x+y) is minimum.

I know if c%gcd(a,b)}==0 then it's always possible. How to find the values of x and y?

My approach

for(i 0 to 2*c):
    x=i
    y= (c-a*i)/b
    if(y is integer)
    ans = min(ans,x+y)

Is there any better way to do this ? Having better time complexity.

Answer 1

Using the Extended Euclidean Algorithm and the theory of linear Diophantine equations there is no need to search. Here is a Python 3 implementation:

def egcd(a,b):
    s,t = 1,0 #coefficients to express current a in terms of original a,b
    x,y = 0,1 #coefficients to express current b in terms of original a,b
    q,r = divmod(a,b)
    while(r > 0):
        a,b = b,r
        old_x, old_y = x,y
        x,y = s - q*x, t - q*y
        s,t = old_x, old_y
        q,r = divmod(a,b)
    return b, x ,y

def smallestSolution(a,b,c):
    d,x,y = egcd(a,b)
    if c%d != 0:
        return "No integer solutions"
    else:
        u = a//d #integer division
        v = b//d
        w = c//d
        x = w*x
        y = w*y
        k1 = -x//v if -x % v == 0 else 1 + -x//v #k1 = ceiling(-x/v)
        x1 = x + k1*v # x + k1*v is solution with smallest x >= 0
        y1 = y - k1*u
        if y1 < 0:
            return "No nonnegative integer solutions"
        else:
            k2 = y//u #floor division 
            x2 = x + k2*v #y-k2*u is solution with smallest y >= 0
            y2 = y - k2*u
            if x2 < 0 or x1+y1 < x2+y2:
                return (x1,y1)
            else:
                return (x2,y2)

Typical run:

>>> smallestSolution(1001,2743,160485)
(111, 18)

The way it works: first use the extended Euclidean algorithm to find d = gcd(a,b) and one solution, (x,y). All other solutions are of the form (x+k*v,y-k*u) where u = a/d and v = b/d. Since x+y is linear, it has no critical points, hence is minimized in the first quadrant when either x is as small as possible or y is as small as possible. The k above is an arbitrary integer parameter. By appropriate use of floor and ceiling you can locate the integer points with either x as small as possible or y is as small as possible. Just take the one with the smallest sum.

On Edit: My original code used the Python function math.ceiling applied to -x/v. This is problematic for very large integers. I tweaked it so that the ceiling is computed with just int operations. It can now handle arbitrarily large numbers:

>>> a = 236317407839490590865554550063
>>> b = 127372335361192567404918884983
>>> c = 475864993503739844164597027155993229496457605245403456517677648564321
>>> smallestSolution(a,b,c)
(2013668810262278187384582192404963131387, 120334243940259443613787580180)
>>> x,y = _
>>> a*x+b*y
475864993503739844164597027155993229496457605245403456517677648564321

Most of the computation takes place in the running the extended Euclidean algorithm, which is known to be O(min(a,b)).