算法学习之路|逆元取模(二)

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

注:可能排版有点辣眼睛,因为要用到一些符号比如“*”,这个又和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的乘积可能会导致溢出。

参考:胡凡 曾磊《算法笔记》

上一篇:Git 最佳实践:分支管理


下一篇:不要在init和dealloc函数中使用accessor