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)).

Answer 2

First let assume a,b,c>0 so:

a.x+b.y = c
x+y = min(xi+yi)
x,y >= 0
a,b,c > 0
------------------------
x = ( c - b.y )/a
y = ( c - a.x )/b
c - a.x >= 0
c - b.y >= 0
c >= b.y
c >= a.x
x <= c/x
y <= c/b

So naive O(n) solution is in C++ like this:

void compute0(int &x,int &y,int a,int b,int c)  // naive
    {
    int xx,yy;
    xx=-1; yy=-1;
    for (y=0;;y++)
        {
        x = c - b*y;
        if (x<0) break;     // y out of range stop
        if (x%a) continue;  // non integer solution
        x/=a;               // remember minimal solution
        if ((xx<0)||(x+y<=xx+yy)) { xx=x; yy=y; }
        }
    x=xx; y=yy;
    }

if no solution found it returns -1,-1 If you think about the equation a bit then you should realize that min solution will be when x or y is minimal (which one depends on a<b condition) so adding such heuristics we can increase only the minimal coordinate until first solution found. This will speed up considerably the whole thing:

void compute1(int &x,int &y,int a,int b,int c)
    {
    if (a<=b){ for (x=0,y=c;y>=0;x++,y-=a) if (y%b==0) { y/=b; return; } }
    else     { for (y=0,x=c;x>=0;y++,x-=b) if (x%a==0) { x/=a; return; } }
    x=-1; y=-1;
    }

I measured this on my setup:

                      x        y      ax+by      x+y     a=50 b=105 c=500000000
[  55.910 ms]        10  4761900  500000000  4761910 naive
[   0.000 ms]        10  4761900  500000000  4761910 opt
                      x        y      ax+by      x+y     a=105 b=50 c=500000000
[  99.214 ms]   4761900       10  500000000  4761910 naive
[   0.000 ms]   4761900       10  500000000  4761910 opt

The ~2.0x difference for naive method times is due to a/b=~2.0and selecting worse coordinate to iterate in the second run.

Now just handle special cases when a,b,c are zero (to avoid division by zero)...