HihoCoder - 1297 數論四·擴展歐幾裏德
描述
小Hi和小Ho周末在公園溜達。公園有一堆圍成環形的石板,小Hi和小Ho分別站在不同的石板上。已知石板總共有m塊,編號為 0..m-1,小Hi一開始站在s1號石板上,小Ho一開始站在s2號石板上。
小Hi:小Ho,你說我們倆如果從現在開始按照固定的間隔數同時同向移動,我們會不會在某個時間點站在同一塊石板上呢?
小Ho:我覺得可能吧,你每次移動v1塊,我移動v2塊,我們看能不能遇上好了。
小Hi:好啊,那我們試試唄。
一個小時過去了,然而小Hi和小Ho還是沒有一次站在同一塊石板上。
小Ho:不行了,這樣走下去不知道什麽時候才匯合。小Hi,你有什麽辦法算算具體要多久才能匯合麽?
小Hi:讓我想想啊。。
提示:擴展歐幾裏德
提示:擴展歐幾裏德
小Hi:首先可以我倆現在的情況列出一個式子:
s1+v1*t=s2+v2*t-k*m (v1<v2)
也就是經過t時間過後,速度快的人剛好超過了速度慢的人k圈,且到達同一個位置。
將這個式子進行變換得到:
(v1-v2)*t+k*m=(s2-s1)
即原式子變成了形如"Ax+By=C"的情況,我們要求解的是一組(x,y)使得原公式成立。
小Ho:形如"Ax+By=C",也就是說我們令A=(v1-v2),B=m,C=(s2-s1),x=t,y=k。
小Hi:恩,沒錯。
求解該式子的算法我們稱為擴展歐幾裏德算法。
該算法分為兩個部分:
(1) 判定是否存在解
對於形如"Ax+By=C"的式子,其存在解的條件為C為A和B最大公約數的整數倍。
我們將A和B的最大公約數記為gcd(A,B)。因此其有解的條件是C=n*gcd(A,B)。
那麽我們應該如何來求解gcd(A,B)呢?
一個樸素的算法是枚舉1~min(A,B),最大的一個能同時被A,B整除的數即gcd(A,B)。顯然這個算法是非常沒有效率的。
為了求解gcd(A,B),歐幾裏德提出了一個輾轉相除法:
首先要證明一個定理:gcd(A,B) = gcd(B, A mod B)
證明: 假設A = k*B+r,有r = A mod B。不妨設d為A和B的一個任意一個公約數,則有A = pd, B = qd。 由r = A - k*B = pd - k*qd = (p - kq)*d,所以有d也為r的約數,因此d是B和A mod B的公約數。 由於對任意一個A和B的公約數都滿足這個性質,gcd(A,B)也滿足,因此有gcd(A,B)=gcd(B,A mod B)。
利用這個性質,我們可以得到算法:
A mod B = 0, 則B為gcd(A,B) A mod B ≠ 0, 則gcd(A,B) = gcd(B, A mod B)
通過不斷的模運算,數據的規模也越來越小,因此能夠快速得收斂到一個解。將其寫成偽代碼為:
gcd(A, B): If (A mod B == 0) Then Return B End If Return gcd(B, A mod B)
(2) 求解
在判定有解之後,我們需要在其基礎上再求出一組(x,y)。由於A,B,C均是gcd(A,B)的整數倍,因此可以將它們都縮小gcd(A,B)倍。即A‘=A/gcd(A,B),B‘=B/gcd(A,B),C‘=C/gcd(A,B)。
化簡為A‘x+B‘y=C‘,gcd(A‘,B‘)=1,即A‘,B‘互質。
此時,我們可以先求解出A‘x+B‘y=1的解(x‘,y‘),再將其擴大C‘倍,即為我們要求的最後解(x,y)=(C‘x‘, C‘y‘)。
那麽接下來我們來研究如何求解A‘x+B‘y=1:
假設A>B>0,同時我們設:
A * x[1] + B * y[1] = gcd(A, B) B * x[2] + (A mod B) * y[2] = gcd(B, A mod B)
已知gcd(A,B)=gcd(B, A mod B),因此有:
A * x[1] + B * y[1] = B * x[2] + (A mod B) * y[2] => A * x[1] + B * y[1] = B * x[2] + (A - kB) * y[2] // A = kB + r => A * x[1] + B * y[1] = A * y[2] + B * x[2] - kB * y[2] => A * x[1] + B * y[1] = A * y[2] + B * (x[2] - ky[2]) => x[1] = y[2], y[1] = (x[2] - ky[2])
利用這個性質,我們可以遞歸的去求解(x,y)。
其終止條件為gcd(A, B)=B,此時對應的(x,y)=(0,1)
將這個過程寫成偽代碼為:
extend_gcd(A, B): If (A mod B == 0) Then Return (0, 1) End If (tempX, tempY) = extend_gcd(B, A mod B) x = tempY y = tempX - (A / B) * tempY Return (x, y)
小Ho:那麽我只需要把A=(v1-v2),B=m,C=(s2-s1),x=t,y=k代入就可以得到t了麽?
小Hi:是的,在已知A,B,C的情況下,我們的確能夠順利求解出一組合法的(x,y)。
但是在求解過程中,我們並沒有保證x是最小的非負整數,它不能直接作為我們的解。
小Ho:那還需要做怎樣的處理麽?
小Hi:我們需要將(A‘,B‘,x‘,y‘)擴充為一個解系。
由於A‘B‘是互質的,所以可以將A‘x‘+B‘y‘=1擴展為:
A‘x‘+B‘y‘+(u+(-u))A‘B‘=1 => (x‘ + uB‘)*A‘ + (y‘ - uA‘)*B‘ = 1 => X = x‘ + uB‘, Y = y‘ - uA‘
可以求得最小的X為(x‘+uB‘) mod B‘,(x‘+uB‘>0)
同時我們還需要將X擴大C‘倍,因此最後解為:
x = (x‘*C‘) mod B‘
若x<0,則不斷累加B‘,直到x>0為止。
那麽最後,小Ho你來總結一下主體部分的偽代碼吧!
小Ho:好的,最後的代碼為:
solve(s1, s2, v1, v2, m): A = v1 - v2 B = m C = s2 - s1 If (A < 0) Then A = A + m // 相對距離變化 End If D = gcd(A, B) If (C mod D) Then Return -1 End If A = A / D B = B / D C = C / D (x, y) = extend_gcd(A, B) x = (x * C) mod B While (x < 0) x = x + B End While Return x
輸入
第1行:每行5個整數s1,s2,v1,v2,m,0≤v1,v2≤m≤1,000,000,000。0≤s1,s2<m
中間過程可能很大,最好使用64位整型
輸出
第1行:每行1個整數,表示解,若該組數據無解則輸出-1
Sample Input
0 1 1 2 6
Sample Output
5
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #include<cstring> 5 #include<algorithm> 6 7 using namespace std; 8 9 long long extend_gcd(long long a,long long b,long long &x,long long &y) 10 { 11 if(b==0) 12 { 13 x=1,y=0; 14 return a; 15 } 16 long long r=extend_gcd(b,a%b,y,x); 17 y-=a/b*x; 18 return r; 19 } 20 21 int main() 22 { 23 long long x,y,m,n,L,ans; 24 while(~scanf("%lld %lld %lld %lld %lld",&x,&y,&m,&n,&L)) 25 { 26 if(m<n) 27 { 28 swap(x,y); 29 swap(m,n); 30 } 31 32 long long a=m-n,b=y-x,X,Y; 33 if(b<0) 34 b+=L; 35 long long d=extend_gcd(a,L,X,Y); 36 if(b%d==0) 37 { 38 X%=L,X+=L,X%=L; 39 ans=X*(b/d)%(L/d); 40 } 41 else 42 ans=-1; 43 if(ans==-1) 44 printf("-1\n"); 45 else 46 printf("%lld\n",ans); 47 } 48 49 return 0; 50 }
HihoCoder - 1297 數論四·擴展歐幾裏德