《算法笔记》——第五章 exgcd 学习记录
- 扩展欧几里得算法(即ax + by = gcd(a,b)的求解)
- 方程ax+by=c的求解
- 同余式ax = c(mod m)的求解
- 逆元的求解以及(b/a)%m的计算。
本节数学证明较多，请读者认真学习。
扩展欧几里得算法
扩展欧几里得算法用来解决这样一个问题：给定两个非零整数a和b，求一组整数解(x,y)，使得ax+by = gcd(a,b)成立，其中ged(a,b)表示a和b的最大公约数。
通过相关定理可知解一定存在。为了讨论问题方便，记gcd = gcd(a,b)，其中a和b为初始给定的数值，因此可以认为在下面讨论的过程中gcd是一个固定的数。
回忆欧几里得算法，如下面的代码所示，它总是把gcd(a, b)转化为求解gcd(b,a%b)，而当b变为0时返回a，此时的a就等于gcd。
也就是说，欧几里得算法结束的时候变量a中存放的是gcd，变量b中存放的是0，因此此时显然有\(a*1 + b*0= gcd\)成立，此时有x=1、y=0 成立。
不妨利用上面欧几里得算法的过程来计算x和y。目前已知的是递归边界成立时为x=1、y=0，需要想办法反推出最初始的x和y。
当计算ged(a,b)时，有\(ax_1+by_1=gcd\)成立；而在下一步计算gcd(b,a%b)时，又有\(bx_2+(a\%b)y_2=gcd\)成立。
因此\(ax_1+by_1=bx_2+(a\%b)y_2\)成立。又考虑到有关系
\(a\%b=a-(a/b)*b\)成立(此处除法为整除)，因此\(ax_1 +by_1=bx_2 +(a-(a/b)*b)y_2\)成立。
将等号右边的式子整理后可得
对比等号左右两边可以马上得到下面的递推公式：
由此便可以通过\(x_2\)和\(y_2\)来反推出\(x_1\)和\(y_1\)了。于是只需要在达到递归边界、不断退出的过程中根据上面的公式计算x和y，就可以得到一组解。 代码如下:
int exgcd(int a,int b,int &x,int &y)
{
if(b == 0)
{
x=1,y=0;
return a;
}
int g=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return g;
}
由于使用了引用，因此当exgcd函数结束时x和y就是所求的解。
显然，在得到这样一组解之后，就可以通过下面的式子得到全部解：
\(k\)为任意整数。
为什么是这样呢?下面简单证明一下。
假设新的解为\(x+s_1\)、\(y-s_2\) ，即有\(a*(x+s_1)+b*(y-s_2)=gcd\)成立，通过代入\(ax+by=gcd\)可以得到\(as_1=bs_2\)，于是\(\frac{s_1}{s_2}=\frac{b}{a}\)成立。为了让\(s_1\)和\(s_2\)尽可能小，可以让分子和分母同时除以一个尽可能大的数，同时保证它们仍然是整数。显然，由于\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)互质，因此\(gcd\)是允许作为除数的最大数，于是\(\frac{s_1}{s_2}=\frac{b}{a}=\frac{b/gcd}{a/gcd}\)，得\(s_1\)和\(s_2\)的最小取值是\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)，证毕。
也就是说，\(x\)和\(y\)的所有解分别以\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)为周期。那么其中\(x\)的最小非负整数解是什么呢？从直观上来看就是\(x\%\frac{b}{gcd}\)。但是由于通过exgcd函数计算出来的x、y可正可负，因此实际上x%-
会得到一-个负数，例如\((-15)\%4=-3\)。考虑到即便\(x\)是负数，\(x\%\frac{b}{gcd}\)的范围也是在\((-\frac{b}{gcd},0)\)，因此对任意整数来说，\((x\%\frac{b}{gcd}+\frac{b}{gcd})\%\frac{b}{gcd}\)才是对应的最小非负整数解。
特殊地，如果\(gcd==1\)，即\(ax+by=1\)时，全部解的公式简化为下式，且x的最小非负整数解也可以简化为\((x\%b + b )\%b\)。
\(k\)为任意整数。
方程ax+by=c的求解
至此已经知道如何求解\(ax + by= gcd\)的解了，那么它有什么应用吗？一般来说，最常见的应用就是利用它来求解ax+by=c，其中c为任意整数。
首先，假设\(ax+by=gcd\)有一-组解\((x_0,y_0)\)，现在在其等号两边同时乘以\(\frac{c}{gcd}\)，即有\(a\frac{cx_0}{gcd}+b\frac{cy_0}{gcd}=c\)成立，因此\((x,y)=(\frac{cx_0}{gcd},\frac{cy_0}{gcd})\)是\(ax+by=c\)的一组解。但是显然这样做的充要条件是\(c\%gcd==0\)，否则第一步 在等号两边同时乘以\(\frac{c}{gcd}\)都无法做到。
于是\(ax + by=c\)存在解的充要条件是\(c\%gcd==0\)，且一组解\((x,y)\)等于\((\frac{cx_0}{gcd},\frac{cy_0}{gcd})\)。可能有些读者会觉得，如果要求全部解，只需要在\(ax+by=gcd\)全部解的基础上都乘以\(\frac{c}{gcd}\)即可，但事实上这只是一部分解而已。
为了获得全部解的公式，可以模仿之前的做法，假设新的解为\(x+s_1\)、\(y-s_2\)，然后将\(a*(x+s_1)+b*(y-s_2)=c\)与\(ax+by=c\)联立，发现同样可以得到\(\frac{s_1}{s_2}=\frac{b}{a}\)成立。于是因为同样的原因，\(s_1\)和\(s_2\)的最小取值仍然是\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)。因此\(ax+by=c\)的全部解的公式为
\(k\)为任意整数。
由此会发现这与\(ax+by=gcd\)全部解的公式是一样的，唯一不同的是初始解\((x,y)\)不同。因此对\(ax+by=c\)来说，其解\((x,y)\)同样分别以\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)为周期。那么为什么说在\(ax + by= gcd\)全部解的基础上都乘以\(\frac{c}{gcd}\)只能获得一 部分解呢?
这其实很简单，由于\(c \ge gcd\)(根据\(c\%gcd==0\)得到)，因此\(ax+by=gcd\)的全部解乘以\(\frac{c}{gcd}\)会导致周期放大为原先的周期放大为原先的\(\frac{c}{gcd}\)倍(\(x\)和\(y\)的周期会分别变成\(\frac{bc}{gcd^2}\)和\(\frac{ac}{gcd^2}\)。
而事实上周期应当是保持\(\frac{b}{gcd}\)与\(\frac{a}{gcd}\)不变的，于是导致漏解。
除此之外，可以得到和上面一样的结论：对任意整数来说，\((x\%\frac{b}{gcd}+\frac{b}{gcd})\%\frac{b}{gcd}\)是\(ax+by=c\)中\(x\)的最小非负整数解，一般来说可以让\(x\)取\(\frac{cx_0}{gcd}\)，其中\(x_0\)是\(ax+by=gcd\)的一个解。
并且，如果\(gcd == 1\)，那么全部解的公式可以化简为下式，且\(x\)的最小非负整数解可以简化为\((x\%b+b)\%b\)。
\(k\)为任意整数。
同余式\(ax \equiv c \pmod m\)
既然已经解决了\(ax +by=c\)，不得不提的就是同余式\(ax \equiv c \pmod m\)的求解。
先解释什么是同余式。对整数\(a\)、\(b\)、\(m\)来说，如果\(m\)整除\(a-b\)，即\((a-b)\%m=0\))，那么就说\(a\)与\(b\)模\(m\)同余，对应的同余式为\(a \equiv b \pmod m\)，\(m\)称为同余式的模。
例如\(10\)与\(13\)模\(3\)同余，\(10\)也与\(1\)模\(3\)同余，它们分别记为\(10 \equiv 13 \pmod 3\)、\(10 \equiv 1 \pmod 3\)。显然，每一个整数都各自与\([0,m)\)中唯一的整数同余。
此处要解决的就是同余式\(ax \equiv c \pmod m\)的求解。根据同余式的定义,有\((ax-c)\%m=0\)成立，因此存在整数\(y\)，使得\(ax-c=my\)成立。移项并令\(y=-y\)后即得\(ax +my=c\)。
由上面的结论，当\(c\%gcd(a,m)==0\)时方程才有解，且解的形式如下，其中\((x,y)\)是\(ax+my=c\)的其中一组解，可以先通过求解\(ax+my=gcd(a,m)\)得到\((x_0,y_0)\)，然后由公式\((x,y)=(\frac{cx_0}{gcd},\frac{cy_0}{gcd})\)直接得到。
\(k\)为任意整数。
虽然对方程\(ax + my=c\)来说，\(k\)可以取任意整数，但是对同余式来说会有很多解在模\(m\)意义下是相同的（由于只关心\(x\)，因此下面只考虑\(x\)）。
对同余式来说，只需要找出那些在模\(m\)意义下不同的解。因此考虑\(x'=x+\frac{b}{gcd(a,m)}*k\)，会发现当\(k\)分别取\(0、1、2、\cdots 、gcd(a,m)-1\)时，所得到的解在模\(m\)意义下是不同的，而其他解都可以对应到\(k\)取这\(gcd(a,m)\)个数值之一。由此可以得到结论：
设\(a,c,m\)是整数，其中\(m \ge 1\)，则
- 若\(c\%gcd(a,m)≠0\)，则同余式方程\(ax \equiv c \pmod m\)无解。
- 若\(c\%gcd(a,m)=0\)，则同余式方程\(ax \equiv c \pmod m\)恰好有\(gcd(a,m)\)个模\(m\)意义下不同的解，且解的形式为
逆元的求解以及\((b/a)\%m\)的计算
接着解决最后一个问题，假设\(a、m\)是整数，求\(a\)模\(m\)的逆元。
先解释什么是逆元（此处特指乘法逆元）。假设\(a、b、m\)是整数，\(m > 1\)，且有\(ab \equiv 1 \pmod m\)成立，那么就说\(a\)和\(b\)互为模\(m\)的逆元，一般也记作\(a \equiv \frac{1}{b} \pmod m\)或\(b \equiv \frac{1}{a} \pmod m\)。通俗地说，如果两个整数的乘积模\(m\)后等于1，就称它们互为\(m\)的逆元。
那么逆元有什么用处呢？对乘法来说有\((b*a)\%m = ((b\%m)*(a\%m))\%m\)成立，但是对除法来说，\((b/ a)\%m = (b\%m)/ (a\%m))\%m\)却不成立，\((b/a)\%m = ((b\%m)/a)\%m\)也不成立。
例如，如果要对\(12/4\)对\(2\)取模，采用\(((12 \% 2)/4) \% 2\)的做法会得到错误的结果\(0\)，而实际上应当是\(1\)。这时就需要逆元来计算\((b/a)%m\)。通过找到\(a\)模\(m\)的逆元\(x\)，就有\((b/a)\%m=(b*x)%m\)成立（只考虑整数取模，也即假设\(b\%a=0\)，即\(b\)是\(a\)的整数倍)，于是就把除法取模转化为乘法取模，这对于解决被除数\(b\)非常大（使得\(b\)已经取过模，不是原始值）的问题来说是非常实用的。
由定义知，求\(a\)模\(m\)的逆元，就是求解同余式\(ax \equiv 1 \pmod m\)，并且在实际使用中，一般把\(x\)的最小正整数解称为\(a\)模\(m\)的逆元，因此下文中提到的逆元都默认为\(x\)的最小正整数解。
显然，同余式\(ax \equiv 1 \pmod m\)是否有解取决于\(1\%gcd(a,m)\)是否为0，而这等价于\(gcd(a,m)\)是否为\(1\)：
- 如果\(gcd(a,m) \ne 1\)，那么同余式\(ax \equiv 1 \pmod m\)无解，\(a\)不存在模\(m\)的逆元。
- 如果\(gcd(a,m)=1\)，那么同余式\(ax \equiv 1 \pmod m\)在\((0,m)\)上有唯一解，可以通过求解\(ax+my=1\)得到。注意：由于\(gcd(a,m)=1\)，因此\(ax+my=1=gcd(a,m)\)，直接使用扩展欧几里得算法解出\(x\)之后就可以用\((x\%m+m)\%m\)得到\((0,m)\)范围内的解，也就是所需要的逆元。
下面的代码使用了扩展欧几里得算法来求解\(a\)模\(m\)的逆元，使用条件是\(gcd(a,m)=1\)，当然如果\(m\)是素数，就肯定成立了，可以放心使用。
int inverse(int a,int m)
{
int x,y;
int g=exgcd(a,m,x,y);
return (x%m+m)%m;
}
另外，如果m是素数，且a不是m的倍数，则还可以直接使用费马小定理来得到逆元，这种做法不需要使用扩展欧几里得算法。
费马小定理：设\(m\)是素数，\(a\)是任意整数且\(a \not\equiv 0\pmod m\)则\(a^{m-1} = 1 \pmod m\)。
使用费马小定理来推导逆元的方法非常简单：由\(a^{m-1} \equiv 1 \pmod m\)可知\(a*a^{m-2} \equiv 1 \pmod m\)，直接由逆元的定义便可以知道\(a^{m-2}\%m\)就是\(a\)模\(m\)的逆元，而这可以通过4.5.3节介绍的快速幂算法很容易求出来，因此不再给出代码。
顺便一提，当\(gcd(a,m) \ne 1\)时，扩展欧几里得算法和费马小定理均会失效，此时\(a\)模\(m\)逆元从概念上来说不存在，但是\((b/a)\%m\)仍然是有值的，此时应当如何求解呢？
再次强调，以下只考虑整数取模，即\(b\)是\(a\)的整数倍的情况。
假设\((b/a)\%m=x\)，因此存在整数\(k\)，使得\(b/a=km+x\)。等式两边同时乘以\(a\)，得\(b=kam+ax\)，于是有\(b\%(am)=ax\)。等式两边同时除以\(a\)，得\((b\%(am))/a=x\)，于是有\((b/a)\%m = (b\%(am))/a\)成立。
因此在\(a\)和\(m\)有可能不互素的情况下，可以使用公式\((b/a)\%m = (b\%(am))/a\) 来计算 \((b/a)%m\)的值，唯一要注意的是\(a\)和\(m\)的乘积可能会太大而导致溢出。因此一般来说尽量使用扩展欧几里得算法或者费马小定理求解，条件不成立时再采用这种办法。
