算法学习之路|逆元取模（二）

原文链接：点击打开链接

摘要： 逆元取模的进一步了解

不胡扯了，一步步来看一下逆元是咋回事吧。

注：可能排版有点辣眼睛，因为要用到一些符号比如“*”,这个又和markdown语法冲突，懒得弄公式编辑器了，便用了代码格式来展示。

首先，补充，我觉得得先了解一些基础知识吧，下午也顺便看了一下：

扩展欧几里得算法

扩展欧几里得算法用来解决这个问题：给定两个非零整数a和b，求一组整数解(x,y),使得ax+by=gcd(a,b)成立（gcd(a,b)表示a,b的最大公约数，因为a,b一开始确定，所以gcd是个定值）。

这里上一段欧几里得算法代码：

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

欧几里得算法结束的时候a变量中存放的是gcd，b变量中存放的是0。

因此有a*1+b*0=gcd成立。此时，x=1,y=0。
证明：设 a>b。
1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；
2，ab!=0 时
设 ax1+by1=gcd(a,b);
bx2+(a mod b)y2=gcd(b,a mod b);
根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;
即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;（注意：这里的a/b是取整除）
根据恒等定理得：x1=y2; y1=x2-(a/b)*y2;
这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.
上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

代码：

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    int gcd=exgcd(b,a%b,x,y);//递归计算exgcd(b,a%b)
    int tem=x;//存放x的值
    x=y;
    y=tem-a/b*y;//这两步即使上面证明得到的恒等式，x1=y2; y1=x2-(a/b)*y2;
    return gcd;//因为x,y是引用传递的，所以最后也会得到x,y的值，也就是所求的一组解
}

扩展欧几里德算法的应用主要有以下三方面：
（1）求解不定方程；
（2）求解模线性方程（线性同余方程）；
（3）求解模的逆元；

（1）解决不定方程(ax+by=c)：

对于不定整数方程pa+qb=c，若 c mod Gcd(p, q)=0,则该方程存在整数解，否则不存在整数解。
上面已经列出找一个整数解的方法，在找到p * a+q * b = Gcd(p, q)的一组解p0,q0后，p * a+q * b = Gcd(p, q)的其他整数解满足：
p = p0 + b/Gcd(p, q) * t 
  q = q0 - a/Gcd(p, q) * t(其中t为任意整数)（提示：g别忘了gcd是最大公约数！想一想就明白了）
至于pa+qb=c的整数解，证明略。直接看结果：
在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后，应该是得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b))，
p * a+q * b = c的其他整数解满足：
p = p1 + b/Gcd(a, b) * t
q = q1 - a/Gcd(a, b) * t(其中t为任意整数。注意：这里的a(b)/Gcd(a,b)并没有变，即周期不变）
p 、q就是p * a+q * b = c的所有整数解。

bool linear_equation(int a,int b,int c,int &x,int &y)
{
    int gcd=exgcd(a,b,x,y);//再次强调，ax+by=c存在解的充分必要条件是c%gcd==0
    if(c%gcd)
        return false;
    int k=c/d;
    x*=k; y*=k;    //求得的只是其中一组解,然后带回p = p1 + b/Gcd(a, b) * t，q = q1 - a/Gcd(a, b) * t求所有解
    return true;
}

特别说明几点：

- 比如,b/gcd是连在一起的，分子和分母！这里电脑格式问题，不能写成上下的格式，（b/gcd)为了看着舒服点，也省略括号了。
- 对于任意整数来说，（x%b/gcd+b/gcd)%b/gcd是ax+by=c中x的最小非负整数解，一般的让x取c*x0/gcd(即上文的p1，偷懒不改了~~）
- 如果gcd==1,那么x的最小非负整数解可以简化为(x%b+b)%b。同样的，通解公式可以化,too:
- x'=x+bK=c*x0+bk,
- y'=y-aK=c*y0-aK,(K为任意整数）

(2)求解模线性方程的方法（同余式 ax≡c(mod m)）：

根据同余式的定义，有(ax-c)%m=0成立，因此，存在正整数y,是得ax-c=my成立。移项并令y=-y，得ax+my=c。艾玛，是不是很熟悉?猜对了！不就回到一个问题了嘛！那行，默认你理解了不定方程的求解，我就直接讲了。

同样的，ax+by=c存在解的充分必要条件是c%gcd==0 。这里c%gcd(a,m)==0时方程才有解。解得形式和上一个问题得出的一模一样（不就是把b换成了m嘛！好像上个问题我是用p,q举例的~~其实是上问是复制粘贴~~侵删）。先通过求解ax+my=gcd(a,m)得到(x0,y0),（代码上文已经给出啦！拿去拿去~~）然后由上文推到出来的公式求出ax+my=c的一组解（就相当于线性代数里面的特解？具体我也记不清了，反正是那个意思），(x,y)=(c*x0/gcd(a,m),c*y0/gcd(a,m))。同样的，通解为：

x'=x+(m/gcd(a,m))*K

y'=y-(a/gcd(a,m))*K  (K为任意整数）

虽然对方程ax+my=c来说，K可以取任意整数，但对于同余式来说会有很多解在在模m的意义下是相同的（由于我们只要求x的值，y是你为了辅助假设出来的，因此下文只考虑x）。对于同余式来说只要考虑那些在模m的意义下不同的解。因此考虑x'=x+m/gcd(a,m)*K，会发现当K分别取0，1，2，.....，gcd(a,m)-1时，所得到的解在模m意义下是不同的，而其他的解都可以对应到K取这gcd(a,m)个数值之一。由此可以得到下面的结论：

设a,c,m是整数，其中m>=1,则：

若c%gcd(a,m)!=0,则同余式方程ax≡c(mod m)无解
若c%gcd(a,m)==0,则同余式方程ax≡c(mod m恰好有gcd(a,m)个模m意义下不同的解，且解的形式为：
x'=x+m/gcd(a,m)*K(x是特解，上文已给出求法）（K=0~gcd(a,m)-1）

代码：

bool modular_linear_equation(int a,int c,int m)
{
    int x,y,x0,i;
    int gcd=exgcd(a,m,x,y);
    if(c%gcd)
        return false;
    x0=x*(c/gcd)%m;   //特解
    for(i=0;i<gcd;i++)
        printf("%d\n",(x0+i*(m/gcd))%m);
    return true;
}

(3)逆元的求解以及(b/a)%m的计算

解释一下啥是逆元，怕你听不懂，我就用一句话说明白：

如果两个数的乘积模m后等于1，就称它们互为逆元。

ab≡1（mod m)，当然m>1,则a,b互为模m  的逆元。

那逆元有什么作用呢？对于乘法来说，有(b*a)%m=((b%m)*(a%m))%m成立，但除法却不成立。所以这时候需要用逆元来计算(b/a)%m。通过找到a模m的逆元x，就有(b/a)%m=(b*x)%m(乘ab(mod m),也就是乘1，所以等价，这里只考虑整数取模，也即假设b%a=0,即b是a的整数倍），于是就把除法取模转化为乘法取模，这样利用乘法取模的性质就可以把b的数值通过取模变小。

由定义可以知道，求解a模m的逆元，就是求解同余式ax≡1(mod m)并且在实际使用中，一般把x的最小整数解称为a模m的逆元，因此下文中提到的逆元都默认为x的最小正整数解。上文我已经强调好多次了，同样的，这里同余式ax≡1(mod m)是否有解取决于1%gcd(a,m)是否为0,而这个等价于：

若gcd(a,m)!=1,那么同余式ax≡1(mod m)无解，a不存在模m的逆元。
若gcd(a,m)==1，那么同余式ax≡1(mod m)在(0,m)上有唯一解（倒数啊，当然不能取0）,可以通过求解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是素数（别问我为什么当然，因为gcd(a,m)要等于1啊，也就是最大公约数要等于1），代码：

int inverse(int a,int  m){
    int x,int y;
    int gcd=exgcd(a,m,x,y);//求解ax+my=1
    return (gcd==1)?(x%m+m)%m;//a模m的逆元为(x%m+m)%m

补充一下：用费马小定理求逆元。

啥叫费马小定理？怕你听不懂，我还是用一句话讲一下：

假如m是质数，且gcd(a,m)=1，那么 a的(m-1)次方≡1（mod m）。即：假如a是整数，m是质数，且a,m互质(即两者只有一个公约数1)，那么a的(m-1)次方除以m的余数恒等于1。

使用费马小定理来推导逆元的方法很简单：由a的(m-1)次方≡1(mod m)可知a*a的(m-2)次方≡1(mod m)，直接由逆元的定义便知道a的(m-2)次方%m就是a模m的逆元，而这个可以用快速幂迅速求出来。

读到这里，该讲的都讲完了，代码也给了。有的读者也许会有疑问了，万一算(b/a)%m的时候a没有逆元呢！！！！那还瞎几把扯了这么多~~（吐血）。当然，还是有办法的。

当gcd(a,m)!=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有可能不互质的情况下可以用此公式来计算），也就转换过来啦。但要注意a和m的乘积可能会导致溢出。