扩展欧几里德

问题：求解 s1 + v1*t = s2 + v2*t - k*m (v1<v2)

已知：s1, s2, v1, v2, m

 求解该式子的算法我们称为扩展欧几里德算法。

 

该算法分为两个部分：

(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)。

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

 

由于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为止。

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

 

解：（v1 - v2）*t + k*m = s2 - s1

　　令 Ａ = v1 - v2, B = m, C = s2 - s1;

　　当 A 小于 0 时， 相对距离变化， A += m

　　所以问题变成求解 A*x + B*y = C

　　令 D = gcd(A, B), 如果 C%D！= 0 则无解

　　A' = A/D, B' = B/D, C' = C/D

　　所以问题变成求解 A'*x + B'*y = C'

　　这时，gcd(A', B') = 1, 即 A', B' 互质

　　令 (x', y') = (x/C', y/C')

　　所以问题变成求解 A'*x' + B'*y' = 1

　　对于这个问题可以由递归解决。

　　因为 gcd(a, b) == gcd(b, a%b) == 1

　　所以 a*x[1] + b*y[1] = gcd(a, b)

　　　　 b*x[2] + (a%b)*y[2] = gcd(b, a%b)

　　所以，x[1] = y[2], y[1] = x[2] - k*y[2]　　（k = a / b）

// 扩展欧几里德    求 Ax + By = C 的一组解 
#include <iostream>
#include <string>
#include <algorithm>
#define ll long long
using namespace std;

typedef struct node{
    ll x, y;
}Node;

ll s1, s2, v1, v2, m;
// 辗转相除法求最大公约数 
ll gcd(ll a, ll b){
    if(a%b==0) return b;
    else return gcd(b, a%b);
}
// 求解 Ax + By = 1 的一组解 Ａ，B互质
// gcd(A, B) = gcd(B, A%B) = 1
// A*x[1] + B*y[1] = B*x[2] + (A%B)*y[2]
// x[1] = y[2], y[1] = x[2] - k*y[2]
Node extend_gcd(ll A, ll B){
    Node t;
    t.x = 0;
    t.y = 1;
    if(A%B == 0){
        return t;
    }
    Node tt = extend_gcd(B, A%B);
    t.x = tt.y;
    t.y = tt.x - (A/B)*tt.y;
    return t;
}
// 求 Ax + By = C 的一组解
ll solve(ll s1, ll s2, ll v1, ll v2, ll m){
    ll A = v1 - v2;
    ll B = m;
    ll C = s2 - s1;
    
    if(A < 0){
        A = A + m;// 相对距离变化
    }
    ll D = gcd(A, B);
    if (C % D){
        return -1;
    }    
    A = A / D;
    B = B / D;
    C = C / D;
    Node tt = extend_gcd(A, B);
    tt.x = (tt.x * C)%B;
    while(tt.x < 0)
        tt.x += B;
    return tt.x;
}
// 求解 （v1 - v2）*t + k*m = s2 - s1 
//    A*t + k*B = C
int main(){
    cin >> s1 >> s2 >> v1 >> v2 >> m;
    
    cout << solve(s1, s2, v1, v2, m) << endl;
} 

View Code