zoukankan      html  css  js  c++  java
  • 扩展欧几里德

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

      

    转载请注明出处:http://www.cnblogs.com/ygdblogs
  • 相关阅读:
    前端工具Gulp的学习
    研究javascript中的this
    如何让引入ES6的html文件运行起来
    windows用命令方式查看文件内容
    windows中用'ls'命令查看项目目录
    一步步理解ajax
    【拥抱ES6】搭建一个ES6环境
    npm还是cnpm
    【聊一聊】css中的经典布局——圣杯布局
    【聊一聊】css中的经典布局——双飞翼布局
  • 原文地址:https://www.cnblogs.com/ygdblogs/p/5476395.html
Copyright © 2011-2022 走看看