模运算以及乘法逆元

Posted on Edited on In Views: Word count in article: 2k Reading time ≈ 7 mins.

基本的模运算

a,b,p均为整数

模p加法：\((a+b)\%p\)

模p减法：\((a-b)\%p\)

模p乘法：\((a*b)\%p\) ## Notes

同余式：正整数a,b对p取模，它们的余数相同，记作\(a \equiv b(mod\space p)\)

基本性质

1.若p=a-b,则\(a \equiv b\mod p\) 例如 \(11\equiv 4\mod 7\)

2.传递性：若\(a\equiv b\mod p\)且\(b\equiv c\mod p\),则\(a\equiv c\mod p\)

运算规则

\[ \begin{split} (a+b)\mod p&=(a\mod p+b\mod p)\mod p\\\\ (a-b)\mod p&=(a\mod p-b\mod p)\mod p\\\\ (a*b)\mod p&=(a\mod p*b\mod p)\mod p\\\\ (a^b)\mod p&=((a\mod p)^b)\mod p\\\\ \nonumber \end{split} \]

结合律

\[ \begin{split} ((a+b\mod p+c)\mod p&=(a+(b+c)\mod p)\mod p\\\\ ((a*b\mod p*c)\mod p&=(a*(b*c)\mod p)\mod p\\\\ \nonumber \end{split} \]

交换律

分配律

\[ (((a+b)\mod p)*c)\mod p=((a*c)\mod p+(b*c)\mod p)\mod p \]

重要定理

\(若a\equiv b\mod p\) 则对于任意的c,都有 \((a+c)\equiv(b+c)\mod p\)

\(若a\equiv b\mod p\) 则对于任意的正整数c,都有 \((a*c)\equiv(b*c)\mod p\)

\(若a\equiv b\mod p,c\equiv d\mod p\) 则对于任意的c,都有 \((a+c)\equiv(b+d)\mod p,(a-c)\equiv(b-d)\mod p,(a*c)\equiv(b*d)\mod p\)

乘法逆元

定义

如果ax≡1 (mod p),且gcd(a,p)=1（a与p互质），则称a关于模p的乘法逆元为x。

例子

例如：4关于1模7的乘法逆元为多少？

4X≡1 mod 7

这个方程等价于求一个X和K，满足

4X=7K+1

其中X和K都是整数。

若ax≡1(mod f), 则称a关于1模f的乘法逆元为x。也可表示为ax≡1(mod f)。

当a与f互素时，a关于模f的乘法逆元有解。如果不互素，则无解。如果f为素数，则从1到f-1的任意数都与f互素，即在1到f-1之间都恰好有一个关于模f的乘法逆元。

例如，求5关于模14的乘法逆元：

14=5*2+4

5=4*1+1

说明5与14互素，存在5关于14的乘法逆元。

1=5-4=5-(14-5*2)=5*3-14

因此，5关于模14的乘法逆元为3。

扩展欧几里得定理

欧几里得算法(辗转相除法)

即以除数和余数反复做除法运算，当余数为0时，取当前算式除数为最大公约数

证明

引理1：若d是a和b的公约数，那么d也是b和c的公约数（c=a%b）

引理2：若d是b和c的公约数(c=a%b)那么d也是b和c的公约数

定理：两个整数的最大公约数等于其中较小的那个数和两数相除余数的最大公约数。

证法：

a可以表示成\(a=kb+r\)(a,b,k,r均为正整数且r<b)

假设d是a，b的一个公约数，记作\(d|a,d|b\)，即a和b都可以被d整除

而\(r=a-kb\)，两边同时除以d，\(r/d=a/d-kb/d\)，由等式右边可知\(m=r/d\)为整数，因此\(d|r\)

因此d也是b，a mod b的公约数

因(a,b)和(b,a mod b)的公约数相等，则其最大公约数也相等

1
2
def gcd(a,b):
    return gcd(b,a%b) if b else a

时间复杂度为O(log n)

本体

设两个式子 \[ \begin{split} ax+by&=gcd(a,b)\\\\ bx'+(a\mod b)y'&=gcd(b,a\mod b) \nonumber \end{split} \] 由欧几里得算法知:\(gcd(a,b)=gcd(b,a\mod b)\)

所以 \(ax+by=bx'+(a\mod b)y'\)

而 \(a\mod b=a-kb\)

其中k=[a/b]

所以有\(ax+by=bx'+(a-kb)y'\)

展开移项得：\(ax+by=ay'+b(x'-ky')\)

根据对应系数关系，可以设

$ x=y'、y=x'-ky'$来进一步求一个可行解

递归地使用x'、y'表示“上一步”的x、y，就能递归地把问题转化为 \[ bx'+(a\mod b)y'=gcd(b,a\mod b) \] x和y会不断减小，类似gcd()的递归终点，当扩展欧几里得算法exgcd()的a\(a'x'+b'y'=gcd(a',b')\)中的\(b'\)为0时，可以得到递归终点的\(x'=1,y'=0\)，层层回溯套用前面等式 \[ \begin{split} x&=y'\\\\y&=x'-ky' \end{split} \] 就能得到\(ax+by=gcd(a,b)\)的一组可行解

算法

对于任意两个整数 a, b，必存在整数 x, y，使得 ax + by = gcd(a,b) 成立，求出整数解 x, y，可以得到 gcd(a,b).

具体步骤

1.定义 x1 = 1, y1 = 0, x2 = 0, y2 = 1

2.若 b ≠ 0，使用带余除法，用 b 除以 a 得到余数 r；否则转到第4步

3.用 b 代替 a, 用 r 代替 b， 令 x1, y1, x2, y2 = x2, y2, x1 - q·x2, y1 - q·y2，重复第2步

4.a 的值就是最大公约数 d，x1 和 y1 的值就是所求 x 和 y

代码实现

1
2
3
4
5
6
def exgcd(a,b):
    if(b==0):
        return 1,0,a
    x,y,g=exgcd(b,a%b)
    x,y=y,(x-(a//b)*y)
    return x,y,g

模拟矩阵运算的扩展欧几里得算法：https://www.cnblogs.com/kentle/p/14975039.html

求逆元

逆元：a关于模b得逆元整数d满足\(a*d\mod b\equiv 1\)

扩展欧几里得：求方程\(ax+by=gcd(a,b)\) 的一组可行解

关系：\((a*d) (mod\space b)\equiv1\)，等价于 ad-kb=1,其中k为未知整数

设d为x，-k为y，则\(ad-kb=1\)转化为\(ax+by=1\)

求出x即可得到a关于模b的逆元

当x为负值的时候

代码实现

1
2
3
4
5
6
7
8
def exgcdinv(a,b):
    if(b==0):
        return 1,0,a
    x,y,g=exgcdinv(b,a%b)
    t=x
    x=y
    y=t-a//b*y
    return x,y,g

时间复杂度为O(log n)（斐波那契复杂度）

补充

1.显然，当x求出负值时，由于乘法逆元应为正，经验证,若\(a*d\mod b\equiv 1\)，d=x+b（证明推导留待时候）

2.python自带的pow函数pow(a,b,c)，当b为-1的时候，可以求a关于模c的乘法逆元

存在即可求

费马小定理

前置知识

费马小定理

对于整数a与质数b，若a与b互质，则有 \[ a^{b-1}\mod b\equiv 1 \]

快速幂取模

两数之积取模等于对两数分别取模然后对积取模 \[ (a*b)\mod p=((a\mod p)*(b\mod p))\mod p \]

代码实现

1
2
3
4
5
6
7
8
def powmod(a,n,mod):
    ret=1
    while n>0:
        if n%2 ==1:
            ret=a*a % mod
        a=a*a % mod
        n//=2
    return ret

求逆元

\(a^{b-1}\mod b\equiv 1\) 等价于 \(a*a^{b-2}\mod b\equiv 1\)

显然\(a^{b-2}\)就是a mod b 的逆元

用b-2和b代替快速幂取模中的n和mod:

1
2
def Fermatinv(a,b):
    return powmod(a,b-2,b)

时间复杂度：大约O(log b)

适用范围：在模数b是质数时

Stein算法

只有整数的移位和加减法，便于大素数求公约数

前置知识

\(gcd(a,a)=a\)，即一个数和他自身的公约数是其自身

\(gcd(ka,kb)=kgcd(a,b)\)也就是最大公约数运算和倍乘运算可以交换，特殊的，当k=2时，说明两个偶数的最大公约数必然能被2整除

算法步骤

假设A1=A,B1=B,C1=1

for n=1,2,3...do:

1.若An=0或Bn=0，则max(An,Bn)*Cn为最大公约数，算法结束

2.若An和Bn均为偶数，则An+1=若An/2，Bn+1=Bn/2，Cn+1=2*Cn

3.若An为偶数，Bn为奇数，则An+1=An/2,Bn+1=Bn

4.若An为奇数，Bn为偶数，则An+1=An，Bn+1=Bn/2

5.若An，Bn均为奇数，则An+1=|An-Bn|/2

Bn+1=min（An,Bn）,Cn+`=Cn

代码实现

1
2
3
4
5
6
7
8
9
10
11
def Stein(a,b):
    if a<b:
        a,b=b,a
        if(b==0):return a
    if a%2==0 and b%2 == 0:
        return 2*Stein(a>>1,b>>1)
    if a%2==0:
        return Stein(a>>1,b)
    if b%2==0:
        return Stein(a,b>>1)
    return Stein(abs(a-b)>>1,min(a,b))