当我真正理解了扩展欧几里得定理

深碍√TFBOYSˉ_ 2022-08-07 09:52 111阅读 0赞

# [当我真正理解了扩展欧几里得定理][Link 1]转载：http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html #

首先、扩展欧几里得定理：对于两个不全为0的整数a、b，必存在一组解x,y，使得ax+by==gcd(a,b);

实现如下：

![ExpandedBlockStart.gif][]

![复制代码][copycode.gif]

int gcd(int a,int b)
    {
           
        int t,d;
        if(b==0)
        {
           
            x=1;
            y=0;   //不明处1
    return a;
        }
        d=gcd(b,a%b);
        t=x;
        x=y;
        y=t-(a/b)*y;  //不明处2
    return d;
    }

![复制代码][copycode.gif]

上面的程序中,x和y我是用全局变量保存的

我个人觉得第一次看到这个程序你会有以上两个不明白的地方（见注释），下面我分别解释

不明处1：由扩展欧几里得定理：ax+by==gcd(a,b)---式1，而此时b==0，也就是说gcd(a,0)==a。原式变为ax+by==a --> x==1,y==0。应该够清楚了吧

不明处2：这里先说明一下我的一些规则，x,y表示第一次递归时的值，x1,y1表示第二次递归时的值。那么

gcd(a,b)==gcd(b,a%b)，同时都代入式1，有ax+by==b\*x1+(a%b)\*y1。将右边变形一下

b\*x1+(a%b)\*y1==b\*x1+(a-(a/b)\*b)\*y1==a\*y1+b\*(x1-(a/b)\*y1)，最终得到ax+by==a\*y1+b\*(x1-(a/b)\*y1)

也就是说，上一深度的x等于下一深度的y1，上一深度的y等于下一深度的x1-(a/b)\*y1。 需要注意，上面推导时用的除法都是整型除法

到这里为止，我们便得到了不定式ax+by==gcd(a,b)的一组解，x、y。

那么对于一般的不定式ax+by==c，它的解应该是什么呢。很简单，x1=x\*(c/gcd(a,b)),y1=y\*(c/gcd(a,b))。很好理解吧~

再深入一点，就解出这么一组解其实一般来说是解决不了什么问题的。没有哪个ACM的题这么简单吧。。。比如我们现在要得到所有的解，那么这所有的解究竟是什么呢？

直接说吧,假设d=gcd(a,b). 那么x=x0+b/d\*t; y=y0-a/d\*t;其中t为任意常整数。

这个是怎么推导出来的，说实话我也不知道，就先这么记着吧！

好了，说了这么多，光说不练也无济于事

先做一个最简单的题目，pku 1061青蛙的约会

先说一下大概题意：有两只青蛙，一只在坐标x，另一直在坐标y，青蛙x一次跳跃可以前进m单位距离，青蛙y一次跳跃可以前进n单位的距离，两青蛙都在同一纬度，该纬度长度为L。两只青蛙同方向同时跳啊跳，问你最少跳多少次，它们才可以相遇，如果不能相遇，输出impossble

分析：假设跳了T次以后，青蛙1的坐标便是x+m\*T,青蛙2的坐标为y+n\*T。它们能够相遇的情况为（x+m\*T）-(y+n\*T)==P\*L，其中P为某一个整数，变形一下

得到(n-m)\*T+P\*L==x-y 我们设a=(n-m),b=L,c=x-y,T=x,P=y.于是便得到ax+by==c。激动啊，这不就是上面一样的式子吗~

直接套用扩展欧几里得函数，得到一组解x,y。由于问题是问最少跳多少次，于是只有x是我们需要的信息。那么再想，x是最小的吗？

答案是可能不是！那么如何得到最小解呢？ 我们考虑x的所有解的式子： x=x0+b/d\*t。x0是我们刚刚求到的，很显然右边是有个单调函数，当t为某一个与x正负性质相反的数时，可以得到最小的x。 令x的正负性质为正，那么x=x0-b/d\*t1 (t1==-t)。令x==0，那么t=x0\*d/b，最小的x等于x0减去t\*b/d。这里得到的x可能是负数，如果是负数，我们再为它加上一个b/d即是所求答案了！

按 Ctrl+C 复制代码

\#include<iostream>  
\#include<string>  
\#include<cmath>  
\#include<algorithm>  
usingnamespace std;  
  
  
\_\_int64 x,y,a,b,c,d;  
\_\_int64 n,m,X,Y,L;  
  
  
\_\_int64 gcd(\_\_int64 a,\_\_int64 b)  
\{  
\_\_int64 t,d;  
if(b==0)  
\{  
x=1;  
y=0;  
return a;  
\}  
d=gcd(b,a%b);  
t=x;  
x=y;  
y=t-(a/b)\*y;  
return d;  
\}  
  
  
int main()  
\{  
while(scanf("%I64d%I64d%I64d%I64d%I64d",&X,&Y,&m,&n,&L)==5)  
\{  
a=n-m;  
b=L;  
c=X-Y;  
d=gcd(a,b);  
if(c%d!=0)  
\{  
printf("Impossible\\n");  
continue;  
\}  
x=x\*(c/d);  
y=y\*(c/d);  
  
  
/\*通解：  
x1=x+b/d\*t;  
y1=y-a/d\*t;  
t为任意整数  
\*/  
//找最小的x1，即求x+b/d\*t最小，那么只有t为某一个数时才最小  
//显然t必须与x正负相反才有最小，那么就看做x-b/d\*t,这个式子的最小值便是t=x/(b/d)时，注意这是整型除法  
\_\_int64 k=x\*d/b;  
k=x-k\*b/d;  
if(k<0)  
k+=b/d;  
printf("%I64d\\n",k);  
\}  
return0;  
\}

再来看一题难一点的 pku 2142 the balance [http://poj.org/problem?id=2142][http_poj.org_problem_id_2142]

还是先说一下题意吧：一个家伙有一种天平，这种天平只有两种重量的砝码a和b，现在要称出重量为c的物品，问你至少需要多少a和b，答案需要满足a的数量加上b的数量和最小，并且他们的重量和也要最小。（两个盘都可以放砝码）

分析：

这题我刚刚开始还以为是动规或者搜索，也算是碰了一鼻子的灰吧。

假设a砝码我们用了x个，b砝码我们用了y个。那么天平平衡时，就应该满足ax+by==c。x,y为正时表示放在和c物品的另一边，为负时表示放在c物品的同一边。

于是题意就变成了求|x|+|y|的最小值了。x和y是不定式ax+by==c的解。

刚刚上面已经提到了关于x,y的所以解的同式，即

x=x0+b/d\*t

y=y0-a/d\*t

你是不是下意识的想要穷举所有解，取|x|+|y|最小的？显然是行不通的，仔细分析：|x|+|y|==|x0+b/d\*t|+|y0-a/d\*t|，我们规定a>b(如果a<b，我们便交换a b)，从这个式子中，我们可以轻易的发现：|x0+b/d\*t|是单调递增的，|y0-a/d\*t|是单调递减的，而由于我们规定了a>b，那么减的斜率边要大于增的斜率，于是整个函数减少的要比增加的快，但是由于绝对值的符号的作用，最终函数还是递增的。也就是说，函数是凹的，先减小，再增大。那么什么时候最小呢？很显然是y0-a/d\*t==0的时候，于是我们的最小值|x|+|y|也一定是在t=y0\*d/a附近了，我是在t点左右5个点的范围内取最小的(据说左右一个点都可以，不过我试了一下wa了)。

![ExpandedBlockStart.gif][]

![复制代码][copycode.gif]

#include<iostream>
    #include<string>
    #include<queue>
    #include<algorithm>
    usingnamespace std;
    #define min(a,b) a<b?a:b
    #define max(a,b) a>b?a:b
    
    int x,y;
    int a,b,c;
    
    int gcd(int a,int b)
    {
           
        int d;
        if(b==0)
        {
           
            x=1;y=0;
            return a;
        }
        else
        {
           
            d=gcd(b,a%b);
            int t=x;
            x=y;
            y=t-(a/b)*y;
        }
        return d;
    }
    
    int main()
    {
           
        int d;
        //freopen("e:\\in.txt","r",stdin);
    while(scanf("%d%d%d",&a,&b,&c)!=EOF)
        {
           
            if(!a &&!b &&!c)
                break;
            int flag=0;
            if(a<b)
            {
           
                flag=1;
                swap(a,b);
            }
            d=gcd(a,b);
            x=x*c/d; y=y*c/d;
            int t=y*d/a,ans=0xfffffff;
            int x1,y1,x2,y2;
            for(int i=t-5;i<=t+5;i++)
            {
           
                x2=x+(b/d)*i;
                y2=y-(a/d)*i;
                if(abs(x2)+abs(y2)<ans)
                {
           
                    ans=abs(x2)+abs(y2);
                    x1=x2;
                    y1=y2;
                }
            }
            if(!flag)
                printf("%d %d\n",abs(x1),abs(y1));
            else
                printf("%d %d\n",abs(y1),abs(x1));
        }
        return0;
    }

![复制代码][copycode.gif]

嘿嘿，每次搞懂新问题都很开心~

再来个更难一点的 PKU 2891 Strange Way to Express Integers [http://poj.org/problem?id=2891][http_poj.org_problem_id_2891]

大概题意是，给你很多同余关系，要你求一个数，使得这个数c mod a1==b1 ,c mod a2==b2,,,,,, c mod ak==bk（a1,a2,a3....ak之间可能不互质）

这个其实就是中国剩余定理的欧几里得解法，如果(a1,a2,,,,,ak)两两互质的话，那么就可以直接套用中国剩余定理的公式了，可是这题可能不互质。那么就只能用扩展欧几里得来解了，对于头两个关系式，c % a1==b1, c % a2==b2，可以联立方程组 a1\*x+a2\*y==b2-b1 ，很熟悉，直接套用扩展欧几里得算法得到一组解x,y，再取x的最小值，带入式子c%a1==b1中，便得到了答案 c 。

这是只有两个式子的，很简单，那么对于剩下的第三个第四个怎么办呢？刚刚得到的x已经满足了头两个同余式，对于新来的第三个式子，他们得到的数必须要前面两个式子的要求，同时还要满足第三个新式子的要求，那么我们便将前两个式子合并成一个式子，用这个新式子与第三个式子联立求解，然后依次类推，知道解出了所有式子的解。

那么如何合并呢？

新式子公式是这样的：c1 mod lcm(a1,a2) ==c 这里c1是要求的新答案，c是前两个式子求的的一个解，lcm(a1,a2)是a1，a2的最小公倍数

你可以理解为新的答案c1即要mod a1,又要mod a2,所以就mod lcm(a1,ba)。而余数c 是刚前面求得的答案，也就是说c1不管怎么折腾，都还会保留一个c，即不会破坏前面的成果！

![ExpandedBlockStart.gif][]

![复制代码][copycode.gif]

#include<iostream>
    #include<string>
    usingnamespace std;
    
    __int64 b,b1,c,c1,d,x,y,A,B,C;
    
    __int64 gcd(__int64 a,__int64 b)
    {
           
        __int64 d,t;
        if(!b)
        {
           
            x=1;
            y=0;
            return a;
        }
        d=gcd(b,a%b);
        t=x;
        x=y;
        y=t-(a/b)*y;
        return d;
    }
    
    int main()
    {
           
        int k,i,flag;
        freopen("in.txt","r",stdin);
        while(scanf("%d",&k)!=EOF)
        {
           
            flag=0;
            scanf("%I64d%I64d",&b,&c);
            if(k==1)
            {
           
                printf("%I64d\n",b+c);
                continue;
            }
            for(i=1;i<k;i++)
            {
           
                scanf("%I64d%I64d",&b1,&c1);
                if(flag)
                    continue;
                A=b; B=b1; C=c1-c;
                d=gcd(A,B);
                if(C%d!=0)
                {
           
                    flag=1;
                    continue;
                }
                x=x*C/d; //取得x的最小值
                x=x-(x*d/B)*(B/d);
                if(x<0)
                    x+=B/d;
                c=x*b+c; b=(b*b1)/d; //按前面的分析头两合并式子
            }
            if(!flag)
            {
           
                printf("%I64d\n",c);
            }
            else
            {
           
                printf("-1\n");
            }
        }
        return0;
    }

![复制代码][copycode.gif]

[Link 1]: http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html
[ExpandedBlockStart.gif]: /images/20220805/1e0081dee60c41d9a6062c8676c44847.png
[copycode.gif]: /images/20220805/f86e0bf06ee0453ca3bbbb7dcaf3c6cb.png
[http_poj.org_problem_id_2142]: http://poj.org/problem?id=2142
[http_poj.org_problem_id_2891]: http://poj.org/problem?id=2891