扩展欧几里德

问题:求解 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

  令 A = 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 的一组解 A,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

 

 

  

posted @ 2016-05-10 09:19  TensionRidden  阅读(279)  评论(0编辑  收藏  举报