欧几里德算法

欧几里德算法。。。用来求最大公约数嘛。。。

咳，至于教练说欧几里德算法“很简单的东西”。。。嘛。。。理解起来确实有困难呢Orz。。。

不过好像最近（明明是一直）考的不少。。。写个博客整理一下。。。

欧几里德算法

首先是欧几里德算法。。。

定义

我们设（射）函数 $gcd(a,b)$ 为求解两个自然数的最大公因数的函数。例如：

$gcd(25,35) = 5$

那么。。。咋求（球）呢。。。

求法推演

首先考虑一下：

对于任意两个正整数 $a,b$ ，都有：

$a=kb+r$

$(k,r∈N)$

所以有：

$r=a\%b$

（在这里，%指的是取余运算）

然后我们假设 $c$ 是 $a$ 和 $b$ 的最大公约数，即

$c=gcd(a,b)$

然后，我们就能得到：

$c|a$

$c|b$

（在这里当然包括除了在这里的任何地方， $a|b$ 表示 $b$ 能够整除 $a$）

然后又因为上面那个式子，有：

$r=a-kb$

所以有：

$c|r$

整合一下上面的式子，我们可以得到：

$c=gcd(b,r)$

即

$gcd(a,b)=gcd(b,a\%b)$

编程实现

然后我们就可以编程咯。。。递归的方法是再好不过了。。。

但是还要考虑一个问题：什么时候停止递归。。。

我们由乘除法的性质，可以得到：0除以任何一个非0的数都得0，即：

$0÷x=0$

$(x∈(-∞,0)∪(0,+∞))$

所以0与一个数的最大公约数是这个数本身，即：

$gcd(a,0)=gcd(0,a)=a$

那就太好办了，反正 $a\%b$ 一定是在不断变小，并且总会有它等于0的时候，那么我们不妨将它设置为递归终点。

CPP代码（递归）

int gcd(int a,int b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}

这种写法就是辗转相除法。

当然为了防止某些情况爆栈（比如说高精度运算什么的），还可以不用递归来写。。。

CPP代码（迭代）

int gcd(int a,int b)
{
    while(b)
    {
        a=a%b;
        swap(a,b);
    }
    return a;
}

当然本质上这两种计算方式都是一模一样的23333

扩展欧几里德算法

然后就是邪恶的扩展欧几里得算法。。。

定义

扩展欧几里德算法旨在解决这样一个问题：对于以下二元一次方程：

$ax+by=gcd(a,b)$

$(a,b∈N)$

求解出其中的一组关于 $x$ 和 $y$ 的解。啊啊啊注意啊这里的一组并不是最小的那一组！！！

然后我们来看看。。。

算法推演

还是从最基本的谈起。。。

在这里提供两种推演方式，一个是自己弄得，另外一个好像是网上比较常见的。。。

推演方式1：

考虑情况1：

$a=b$的情况

我们可以得到：

$ax+ay=a$

所以：

$x+y=1$

所以：

$x=1$
$y=0$

或者

$x=0$
$y=1$

考虑情况2：

$b=0$的情况

当$b=0$时，必有：

$ax+0y=gcd(0,a)$

即：

$ax=a$

所以：

$x=1$
$y=-1,0,1,2,3,4,......,+∞$

这样，无论$y$取何值，都是一组有效的整数解。所以为了方便，我们取$y=0$，即：

$x=1$
$y=0$

考虑情况3：

$a,b ≠0$的情况

由那个二元一次方程，我们可以推导出：

$ax+by=bx'+(a\%b)y'$

（假设 $x',y'$ 是另外一组解）

把中间那个“$a\%b$”拆开（注意那个除号是整除）：

$ax+by=bx'+(a-a/b·b)y'$

然后稍微整理一下：

$ax+by=bx'+ay'+(a/b·b)y'$

然后再稍微整理一下，就得到：

$ax+by=ay'+b(x'-a/b·y')$

然后就可以得到：

$x=y'$
$y=x'-a/b·y'$

然后继续往下算，递归计算，直到遇到前面说的两种情况就好了。

推演方式2：

或者说你看不惯我的证明放肆方式，那么下面还有一个，是我同学当时写的一个手记，我用马克当写了一下：

以下转载自同学的手记

已知原方程 $ax+by=gcd(a,b)$
当$a=b$时 :
$x=1,y=0$ 或者 $x=0,y=1$
当$a>b$时 :
当$b=0$时 : $gcd(a,b)=a$ $x=1,y=0$
$a·b≠0$时：
设以下方程组：
$a_1x_1+b_1y_1=gcd(a_1,b_1)$
$a_2x_2+b_2y_2=gcd(a_2,b_2)$
$a_3x_3+b_3y_3=gcd(a_3,b_3)$
······
$a_ix_i+b_iy_i=gcd(a_i,b_i)$

那么$a_1=a$ $b_1=b$ $(x_1,y_1)$即一组解
因为$gcd(a_i,b_i)=gcd(b_i,a_i\%b_i)$
所以$a_i=b_{i-1}$ ，$b_i=a_{i-1}\%b_{i-1}$
所以$gcd(a_i,b_i)=gcd(b_i,a_i\%b_i)=gcd(a_{i+1},b_{i+1})$

又因为$gcd(a_i,b_i)=gcd(a_2,b_2)=······=gcd(a_i,b_i)$
所以$a_1x_1=b_1y_1=a_2x_2+b_2y_2$
所以$a_2=b_1$ $b_2=a_1\%b_2=a_1-(a_1/b_1)b_1$
所以$a_1x_1+b_1y_1$
$=b_1x_2+(a_1-(a_1/b_1)b_1)y_2$
$=a_1y_2+b_1(x_2-(a_1/b_1)y_2)$
所以$x_1=x_2$ ，$y_1=x_2-(a/b)y_2$
依次类推：
求得 $(a/gcd)x+(b/gcd)y=1$
所以$(a/gcd)(x+(b/gcd))+(b/gcd)(y-(a/gcd))=1$
$x$的周期是$b/gcd$

编程实现

按照上面所述的$a,b ≠0$的情况进行，直到遇到前面说的两种递归终点就好了。

CPP（定义）

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    int nx,ny;
    int q=exgcd(b,a%b,nx,ny);
    x=ny;
    y=nx-a/b*ny;
    return q;
}

在这里，函数的返回值即是$a$和$b$的$gcd$。

使用地址引用的方式传入$x$和$y$两个变量，执行完函数后，原先的$x$和$y$变量就已经被存入了一组解。

CPP（精简）

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    int q=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return q;
}

关于这个程序对于非递归终点的情况处理为什么如此精神焕发活力四射姬情满满奇怪，在这里做一下解释哈

为什么再下一层递归中要将本层递归的y作为下一层递归的x，将本层递归的x作为下一层递归的y呢。。。

别着急，先接着往下看。

你会注意到，将本层的$x$作为下一层的$y$，正好满足了上面的$x=y'$这一递推式。

然后，$y$被减去一个$a/b*x$，即是下面的递归结束后，返回的$x'$值被减去一个$a/b*y'$。

即是$y=x'-a/b·y'$，也满足了递推关系式。

所以总结一下：我们直接将本层递归的$x$和$y$变量交换传入下一层，达到了以下目的：

• 省略掉了中间变量占用的时间和空间复杂度

• 最大化利用了函数的地址引用

然而更推荐第一种写法。。。并不怎么容易写错。。。第二种看似精简，其实很容易出错的。。。（其实我本人很喜欢第二种的说。。。但是没办法，在考场上确实很容易写错）

关于最小正整数解的问题

先普及一下，最小正整数解指的是其中一个数是正整数就好啦，不用两个都是。。。

我们来想一下如何将$x$值变为最小正整数。

在原有方程左右两侧同除一个$gcd(a,b)$，得（以下将$gcd(a,b)$简化为$gcd$）：

$(a/gcd)x+(b/gcd)y=1$

然后在左右两侧同加一个$ab/gcd^2$，并移项，合并同类项，得到：

$(a/gcd)(x+(b/gcd))+(b/gcd)(y-(a/gcd))=1$

然后我们得到：每当$x$的值增加$b/gcd$时，都需要$y$的值对应减小$a/gcd$来维持等式。

所以，满足的解集：

$x,y$
$x+(b/gcd),y-(a/gcd)$
$x+2(b/gcd),y-2(a/gcd)$
$x+3(b/gcd),y-3(a/gcd)$

那么我们称，这个方程的解具有周期性，$x$的周期为$b/gcd$，$y$的周期为$a/gcd$。

也就是说，满足( $x+n(b/gcd),y-n(a/gcd)$ )就是有效的解。

然后我们对得出来的$x$值进行一下这样的操作：

1.对$b/gcd$取余
2.加上$b/gcd$
3.对$b/gcd$取余

这样得出来的$x$值一定是最小整数解。然后我们再推导一下怎么求$y$就行了。

为啥？因为如果我们对$b/gcd$取余，如果是负数，就会得出负的结果，不满足要求，所以要再加上$b/gcd$，使它变成正数。然后如果一开始是正数，那么一加就不是最小的整数解了，所以再.对$b/gcd$取余。这样无论是正数还是负数，都能得到最小正整数解。

如果想让$y$值是最小正整数解，只要对$y$进行相同的操作，但是注意是对$a/gcd$取余和相加就行了。即：对$y$进行如下操作：

1.对$a/gcd$取余
2.加上$a/gcd$
3.对$a/gcd$取余

然后再用相同的方法把$x$算出来就行了。

应用

求解不定方程（二元一次方程）的最小正整数解

这玩意绝逼是个坑。。。

内容

求解形如：

$ax+by=c$

形式的有整数解的方程。

这个东西。。。已经没有什么好说的了。。。一开始以为多么高端的一个东西。。。结果发现它只能解出有整数解的二元一次方程。。。

算法

算法。。。很简单。。。

一开始的时候，扩展欧几里德算法算出来的不是形如：

$ax+by=gcd(a,b)$

的方程嘛。。。

解方程的方法也超级简单，你只需要先算出$c$是$gcd(a,b)$的几倍，然后将扩展欧几里德求出来的$x$和$y$分别乘以这个求出来的倍数。。。就是答案了。。。（喵了个咪的就是这个东西害我以为欧几里得算法是螺旋上天法力无边结果发现原来只是这样一个东西）

证明正确性

然后我们要注意了。。。由于要解出的是整数解。。。所以就必须要满足：

$gcd(a,b)|c$

否则方程是木有整数解滴。。。

道洗呆？

先简化这个问题。首先来讲，从上面的解法中可以看出，右边的$c$变化，在有整数解的情况下，$x$和$y$只是对应的放大。

所以问题就成了为什么只有满足：

$n·ax+n·by=n·gcd(a,b)$

才有正整数解。。。

分情况讨论一下：

• $a$和$b$相等的情况：

这种情况下，方程的形式就变成了：

$n·a(x+y)=n·a$

的形式。所以如果右边一旦不是$n·a$的形式，就无法保证方程有整数解。

• $a$和$b$不相等的情况：

  • $a$和$b$的最大公因数等于1的情况：这种情况不管怎么样都有整数解，想想为啥。

  • $a$和$b$的最大公因数不等于1的情况：

这种情况下，我们设$a<b$。反正如果不满足的话互换他俩就是咯。

然后我们想：如果$a$和$b$的最大公约数不为1，那么一定满足：

$a=n·gcd(a,b)$
$b=m·gcd(a,b)$

所以在有整数解的方程$ax+by=c$中，一定满足：

$n·x·gcd(a,b)+m·y·gcd(a,b)=q·gcd(a,b)$

即：

$q=nx+my$

即：c的值一定是最大公约数的倍数，证毕。

编程实现

原理很简单，遵循算法按照定义来就好了。。。

这里就不再重复写exgcd的代码了，直接引用上面给出的就好啦。

CPP代码

bool eqa_solve(int &x,int &y,int a,int b,int c)
{
    int q=exgcd(a,b,x,y);
    if(c%q!=0) return false;
    int tms=c/q;
    x*=tms;
    y*=tms;
    return true;
}

解释用法：将要存储结果的变量填写在$x$和$y$的位置，然后按照方程$ax+by=c$的形式填写$a、b、c$三个位置，返回值为方程是否有整数解。

关于最小正整数解的问题

跟扩展欧几里德算法还是这玩意：

$(a/gcd)(x+(b/gcd))+(b/gcd)(y-(a/gcd))=1$

这个方程的解同样具有周期性，$x$的周期为$b/gcd$，$y$的周期为$a/gcd$。

也就是说，满足( $x+n(b/gcd),y-n(a/gcd)$ )就是有效的解。

然后我们对得出来的$x$值进行一下这样的操作：

1.对$b/gcd$取余
2.加上$b/gcd$
3.对$b/gcd$取余

这样得出来的$x$值一定是最小整数解。然后我们再推导一下怎么求$y$就行了。

如果想让$y$值是最小正整数解，只要对$y$进行相同的操作，但是注意是对$a/gcd$取余和相加就行了。即：对$y$进行如下操作：

1.对$a/gcd$取余
2.加上$a/gcd$
3.对$a/gcd$取余

然后再用相同的方法把$x$算出来就行了。

求乘法逆元

定义

先说什么是乘法逆元。

一般来讲，如果要运算加法、减法、乘法、乘方，都应该满足以下式子：

$(a+b)\%c=(a\%c+b\%c)\%c$

$(a-b)\%c=(a\%c-b\%c)\%c$

$(a·b)\%c=(a\%c·b\%c)\%c$

$a^b\%p =(a\%p)^b\%p$

然而这里出现了一个问题：

如果是除法，并不满足$(a/b)\%c=(a\%c/b\%c)\%c$，不信你代个数试试：

$(6/3)\%3=2≠(6\%3/3\%3)\%3=RuntimeError$

那怎么实现对除法的取余呢？

这里就引入乘法逆元这个东西。

他可以达到这样的效果：

$(a/b)\%k=(a·c)\%k$

你可以将它简单地理解为类似于倒数的东西，只不过是再对倒数取余而已，即：

$(b·c)\%k=1$

所以注意，逆元是针对一个数而言的，并不是针对一个表达式。

求法

使用扩展欧几里德算法，但是好像只能在所求的数和取余的数互质才行。

因为：

$(b·c)\%k=1$

所以我们可以得到

$b·c+n·k=1$

所以只有在$b$和$k$互质的情况下，才能得到：

$b·c+n·k=gcd(b,k)$

然后用这个方程求解出$c$的值，就是所求的乘法逆元。

编程实现

CPP代码

int mod_inverse(int a)
{
    int x,y;
    exgcd(a,mod,x,y);
    return (mod+x%mod)%mod;
}

当然，这个还是求得最小正整数解。

解线性同余方程

定义

额。。。同余方程是什么。。。

同余方程 $ax≡b$ $(mod$ $m)$ (也就是 $ax \%m = b\%m$)

也就是说，求解

$ax+my=b\% m$
（y是某个可爱的未知数）

解法

太好办了。。。

使用欧几里德算法解不定方程，将方程：

$ax+by=c$

中的$a$设置为上面的$a$，$b$设置为$m$，$c$设置为$b\% m$

然后解方程就行啦。。。当然同样的，如果要解的不定方程没有整数解的话。。。那自然同余方程就没有整数解啦。。。

编程实现

直接带入求方程就行了。。。

CPP代码

bool congruence_solve(int a,int &x,int b,int m)
{
    int y;
    return eqa_solve(x,y,a,m,b%m);
}

返回值代表该方程是否被解出。