求解線性方程

Question

我必須找出方程ax+by=c的積分解，使得x>=0和y>=0和值爲(x+y) is minimum。

我知道如果c%gcd(a,b)}==0那麼它總是可能的。如何找到x和y的值？

我的做法

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

有沒有什麼更好的方法來做到這一點？有更好的時間複雜性。

Answer 1

首先，讓我們假設a,b,c>0這樣：

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 

天真O(n)的解決辦法是在C++這樣的：

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; 
    } 

如果沒有解決找到返回-1,-1如果你想想方程那麼你應該認識到，最小的解決方案將是當x或y是最小的（哪一個取決於a<b條件），所以增加這樣的啓發式，我們可以只增加最小座標，直到找到第一個解。這將加速大幅整個事情：

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; 
    } 

我衡量這對我的設置：

     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 

幼稚方法次~2.0x差異是由於a/b=~2.0和選擇更糟糕的協調迭代中第二輪。

現在只是處理特殊情況時a,b,c爲零（由零避免師）...

Answer 2

使用Extended Euclidean Algorithm和linear Diophantine equations的理論就沒有必要進行搜索。這裏是一個Python 3的實現：

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) 

典型運行：

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

它的工作方式：先使用擴展歐幾里德算法找出d = gcd(a,b)和一個解決方案，(x,y)。所有其他解決方案的格式爲(x+k*v,y-k*u)，其中u = a/d和v = b/d。由於x+y是線性的，因此它沒有臨界點，因此在x儘可能小或者y儘可能小時在第一象限中被最小化。上面的k是一個任意的整數參數。通過適當使用floor和ceiling，您可以找到儘可能小的x或y儘可能小的整數點。只要拿最小的一個。

On編輯：我的原始代碼使用了應用於-x/v的Python函數math.ceiling。這對於非常大的整數是有問題的。我調整了它，以便只用int操作來計算天花板。它現在可以處理任意大數：

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

大部分的計算髮生在運行擴展歐幾里德算法，這是衆所周知的是O(min(a,b))。