模意义下的数学

填完了线性代数积性函数和莫比乌斯反演的坑之后,来补一下取模的笔记。

前置导入

什么是取模?

取模运算是求两个数相除的余数。
取模运算(“Modulo Operation”)和取余运算(“Remainder Operation ”)两个概念有重叠的部分但又不完全一致。主要的区别在于对负整数进行除法运算时操作不同。取模主要是用于计算机术语中。取余则更多是数学概念。
模运算在数论和程序设计中都有着广泛的应用,奇偶数的判别到素数的判别,从模幂运算到最大公约数的求法,从孙子问题到凯撒密码问题,无不充斥着模运算的身影。虽然很多数论教材上对模运算都有一定的介绍,但多数都是以纯理论为主,对于模运算在程序设计中的应用涉及不多。
(摘自百度

前置知识

mod\bmod 表示取模运算,也就是两数相除取其余数,在 C++ 中用 % 表示。amodba\bmod{b} 表示求 aa 除以 bb 的余数。下面这个式子 ab(modp)a\equiv b\pmod{p} 可以理解为 amodp=bmodpa\bmod{p}=b\bmod{p}。模运算与乘除运算同等级。但是本文有些地方为了方便就将模运算提升到比乘除高一个优先级的地方去了。

Anm\operatorname{A}_n^m 表示从 nn 个不同的元素中任意取出其中 mm 个元素的所有排列数的个数。其中一定有 mnm\leqslant n,若 m>nm>n,则 Anm=0\operatorname{A}_n^m=0

Cnm\operatorname{C}_n^m 表示从 nn 个不同的元素中任意取出其中 mm 个元素且并称一个组合(不区分顺序)的方案数。其中一定有 mnm\leqslant n,若 m>nm>n,则 Anm=0\operatorname{A}_n^m=0

充分条件:对于某一个命题,若条件成立的时候命题一定正确(命题正确的时候条件成立与否不管),则称这个条件是这个命题的充要条件。

必要条件:对于某一个命题,若命题成立的时候条件一定成立(条件成立的时候命题成立与否不管),则称这个条件是这个命题的必要条件。

充分必要条件:字面意思。

其他的出现的初中没有的符号都在前两期数学有关的博客中:线性代数积性函数和莫比乌斯反演

关于本文

  • 如果有误,洛谷私信。
  • 本文不带有任何网上找的例题,如果有需要自行查找。
  • 没有提到的变量范围请根据上下文情景判断,一般为有理数范围或正整数范围。
  • 本文仅为信息学竞赛基本需要掌握的内容,不代表其学科的全部。
  • 大部分内容源自校内讲课和网络资源。
  • 因为这里模块较多且涵盖许多有关模运算的算法,对模运算本身没有过多介绍。

模意义与逆元

模意义及其性质

小学的时候我们都做过整数除法,例如 27÷4=6327\div 4 = 6 \cdots 3。这里的 2727 是被除数,44 是除数,66 是商,33 是余数(都是老生常谈了吧)。模运算的表示为 273(mod4)27 \equiv 3 \pmod{4}。表示 272733 在模 44 意义下是相等的。因为 272733 显然本来不相等,但这两个数除以 44 的余数相同,就可以称 272733 在模 44 意义下是相等的。

确切地说,对于三个数 a,b,pa,b,p,如果 a÷pa\div pb÷pb\div p 的余数相同,就可以记作 ab(modp)a\equiv b\pmod{p}(有时为了方便阅读可能会有 b<pb<p 的规定),也就是 aabb 在模 pp 意义下是相同的。

算法竞赛下的很多题目,如果答案实在过大,但是善良的出题人并不想给这道题整一个高精度,就会运用到模运算。例如将答案对 998244353998244353 取模。意思就是你只要算出最终答案 AnsAns 除以 998244353998244353 的余数就可以了。当然除非出题人头脑发热,不然不可能会遇到求 AnsAns00 取模的。(不过对 11 取模也不是很可能的亚子)

不过如果你执意要用高精度计算出 AnsAns,然后写一个高精度取模的话……那你可能不是很理解一下三个性质:

(a+b)modp=(amodp+bmodp)modp(a+b)\bmod{p} = (a\bmod{p}+b\bmod{p})\bmod{p}

(ab)modp=(amodpbmodp)modp(a-b)\bmod{p} = (a\bmod{p}-b\bmod{p})\bmod{p}

(a×b)modp=(amodp×bmodp)modp(a\times b)\bmod{p} = (a\bmod{p}\times b\bmod{p})\bmod{p}

有了这样三个性质,一般的计算就可以一步一步对 pp 取模,以此避免溢出。但是千万不能想当然地整出一个 (a÷b)modp=(amodp÷bmodp)modp(a\div b)\bmod{p}=(a\bmod{p}\div b\bmod{p})\bmod{p}。举个例子:24mod13=9,12mod13=12,2mod13=224\bmod{13}=9,12\bmod{13}=12,2\bmod{13}=2,显然 9÷2129\div 2\neq 12

逆元

既然加减乘都满足这个东西,除法不满足。但是如果我们得到的式子需要用除法该怎么办?这个时候就相当于,我们需要找到一个数 kk,满足 a÷ba×k(modp)a\div b\equiv a\times k\pmod{p},此时这个 kk 就叫 bb 在模 pp 意义下的逆元。

如何求出这个逆元呢?如果你直到费马小定理,就可以很快算出来如何求一个逆元:当 pp 为质数时,有:

ap11(modp)a^{p-1}\equiv 1\pmod{p}

简单变一下形:

aap21(modp)a\cdot a^{p-2}\equiv 1\pmod{p}

这时我们发现 ap2a^{p-2}aa 相乘在模 pp 意义下恰好为 11。因此,aaap2a^{p-2} 在模 pp 意义下互为倒数,而除一个数等于乘一个数的倒数。因此对于一个数 kk,模 pp 意义下的 k÷ak\div a 可以转化为 k×ap2k\times a^{p-2}

从而我们得到了一个很重要的性质:若模数 pp 为质数,有 a÷bq×bp2(modp)a\div b\equiv q\times b^{p-2} \pmod{p}bb 的逆元就是 bp2b^{p-2}

如果 pp 不是质数的话,以上方法就行不通了。但是考虑到欧拉定理是对于所有数 pp 都适用的(不过也有一个前提就是互质,具体见下)。也就是对于互质的两个数 a,ma,m,有 aφ(m)1(modm)a^{\varphi(m)}\equiv 1\pmod{m}。费马小定理求逆元可以用快速幂,是 O(logn)\operatorname{O}(\log n) 的,后者需要求欧拉函数值,还要快速幂,求单个欧拉函数值需要 O(n)\operatorname{O}(\sqrt{n}),总的复杂度就是 O(nlogn)\operatorname{O}(\sqrt{n}\log n)。如果考虑预处理欧拉函数值,则是 O(n)\operatorname{O}(n) 预处理,单次操作 O(logn)\operatorname{O}(\log n)

扩展欧几里得

数论题中经常会要求你求出一个这样的不定方程的整数解(其中参数也是整数):

ax+by=cax+by=c

这时如果直接暴力枚举显然不现实。因此我们需要一个优秀的算法解决这类问题。

前置知识

在解这个方程之前,不妨先看一看怎么判断它有没有解。

裴蜀定理


对于不定方程 ax+by=cax+by=c,记 d=gcd(a,b)d=\gcd(a,b),则 x,yx,y 有整数解,当且仅当 dcd\mid c

裴蜀定理的一个重要的推论是:a,ba,b 两数互质的充分必要条件是存在一对整数 x,yx,y,满足 ax+by=1ax+by=1


辗转相除法


辗转相除法又叫做欧几里得算法。一般用于计算 gcd(a,b)\gcd(a,b),因为其求解 gcd\gcd 比暴力枚举公因数快得多,所以一般的算法竞赛都会使用这种算法。(如果不是,那可能就是用的扩展欧几里得)算法整体非常简单,只需要记住一个公式:

gcd(a,b)=gcd(b,amodb)\gcd(a,b)=\gcd(b,a\bmod{b})

直到 b=0b=0 时,aa 的值就是它们的最大公约数。公式看起来十分简单,而代码实现也就是递归,只有两行(当然如果你要杠说这是四行我也不管):

inline int gcd(int a,int b){
    if(!b) return a;
    return gcd(b,a%b);
}

这时我们已经可以对这个方程有一个无解与否的判断了:如果原方程没有整数解,则 gcd(a,b)c\gcd(a,b)\nmid c,而 gcd(a,b)\gcd(a,b) 的求法可以用辗转相处法简单实现。因此可以正式来看一看这个方程到底怎么解了。

扩展欧几里得实现

扩展欧几里得是一个非常强大的算法,它可以在计算出 ax+by=cax+by=c 的解的同时,求出 gcd(a,b)\gcd(a,b)(当然你也可以直接用辗转相除法求最大公约数)。

扩展欧几里得算法的基本思想,就是在通过辗转相除法求得 gcd(a,b)\gcd(a,b) 的同时,通过辗转相处过程中的一些式子,倒回去推导出 x,yx,y 的值。

先不妨假设我们已经找到了一组特殊解 (x0,y0)(x_0,y_0),并求出来了 k=gcd(a,b)k=\gcd(a,b),那么我们就可以用这一组特解表示所有的解(其中 tt 是参数):

x=x0+(b÷k)×t,y=y0(a÷k)×tx=x_0+(b\div k)\times t,y=y_0-(a\div k)\times t

为什么不是 x=x0+b×t,y=y0a×tx=x_0+b\times t,y=y_0-a\times t 呢?显然,tt 是一个整数,随着 tt 的变化,上面那两个式子可以得到的 (x,y)(x,y) 一定比这一个多一些。再者,按照上面这种方法,会发现 gcd(b÷k,a÷k)=1\gcd(b\div k,a\div k)=1。此时也还可以得到系数的最小值就是 b÷kb\div ka÷ka\div k

现在来考虑如何求出一组特解 x0,y0x_0,y_0。因为扩展欧几里得是从辗转相除的时候得出了一些式子,最后倒推回来做的。不妨考虑辗转相除法的结束状态,此时 b=0,a=gcdb=0,a=\gcd。此时只需要给 aa 一个为 11 的系数,也就是 1×a+y×b=c1\times a + y\times b=c,因为 bb00,所以 yy 随便是多少都可以,只要 x=1x=1,则 c=gcd(a,b)c=\gcd(a,b),根据裴蜀定理,此时一定有解。

此时考虑如何一步一步推回最初的状态。对于某一个状态 (a,b)(a,b),因为是倒着推的,所以我们已经求出了下面一个状态 (b,amodb)(b,a\bmod{b}) 的解,也就是一组 x0,y0x_0,y_0 使得:

b×x0+(amodb)×y0=gcdb\times x_0+(a\bmod{b})\times y_0=\gcd

显然 amodb=aab×ba\bmod{b}=a-\left\lfloor\dfrac{a}{b}\right\rfloor\times b,所以有:

gcd=b×x0+y0×(aab×b)=b×x0+a×y0ab×b×y0\begin{aligned} \gcd & = b\times x_0 +y_0\times (a-\left\lfloor\dfrac{a}{b}\right\rfloor\times b)\\ & =b\times x_0 + a\times y_0 - \left\lfloor\dfrac{a}{b}\right\rfloor\times b\times y_0 \\ \end{aligned}

对比我们要求的当前状态 (a,b)(a,b),使 ax+by=gcdax+by=\gcd,发现这里其实就是:

xy0x\leftarrow y_0

yx0ab×y0y\leftarrow x_0-\left\lfloor\dfrac{a}{b}\right\rfloor\times y_0

因此就得出一个状态是如何从另一个状态得出的了。代码实现如下:

int x1,x2,res; // x1 表示 x,x2 表示 y
inline int exgcd(int a,int b){
    if(!b){
        x1=1,x2=0;
        return a;
    }
    res=exgcd(b,a%b);
    int x3=x1;
    x1=x2,x2=x3-x2*(a/b);
    return res;
}

关于扩展欧几里得,它除了求解 ax+by=cax+by=c 的方程之外,还有很多别的用处。需要每一个人慢慢探索。例如:


引导例题——扩展欧几里得求逆元。给出 a,pa,p,求 aa 在模 pp 意义下的逆元,不保证 pp 是质数或 ppaa 互质。

我们发现记 xxaa 的逆元,也就有 a×x1(modp)a\times x\equiv 1\pmod{p},转换一下就可以将其转化为求不定方程 ax+py=1ax+py=1 的整数解,扩展欧几里得求即可。


扩展中国剩余定理

来趁热打铁地了解一下扩展中国剩余定理。不过在此之前还要看一下什么叫做中国剩余定理

孙子定理是中国古代求解一次同余式组(见同余)的方法。是数论中一个重要定理。又称中国余数定理。一元线性同余方程组问题最早可见于中国南北朝时期(公元5世纪)的数学著作《孙子算经》卷下第二十六题,叫做“物不知数”问题。(文章摘自百度)、

中国剩余定理需要解决这样一类问题:

求解下面这一个一元线性同余方程组的:

{xa1(modp1)xa2(modp2)xa3(modp3)xan(modpn)\begin{cases} x\equiv a_1\pmod{p_1} \\ x\equiv a_2\pmod{p_2} \\ x\equiv a_3\pmod{p_3} \\ \quad\vdots \\ x\equiv a_n\pmod{p_n} \end{cases}

其中 p1,p2,,pnp_1,p_2,\cdots,p_n 两两互质

前置知识

  1. 如果有 ab(modc)a\equiv b\pmod{c},则有 a÷kb÷k(modc)a\div k\equiv b\div k\pmod{c}
  2. 如果有 ab(modc)a\equiv b\pmod{c},则有 a+k×cb(modc)a+k\times c\equiv b\pmod{c}

中国剩余定理

举个例子,假设需要计算一个数,它除以 3322,除以 5533,除以 7722。这个数是多少?

假设我们对每一个方程写出一个解 x1,x2,x3x_1,x_2,x_3,满足 x12(mod3),x23(mod5),x32(mod7)x_1\equiv 2\pmod{3},x_2\equiv 3\pmod{5},x_3\equiv 2\pmod{7}

考虑从 x1x_1 角度出发,能不能找到某个 x2x_2,使得 x1+x22(mod3)x_1+x_2\equiv 2\pmod{3} 呢?根据前置知识的第二点,如果 3x23\mid x_2 就可以满足条件。同理,如果我要再找到某个数 x3x_3,使得 x1+x2+x32(mod3)x_1+x_2+x_3\equiv 2\pmod{3},则需要 3x33\mid x_3

同理,从 x2,x3x_2,x_3 的角度出发,总共可以得到下面这三个条件:

5x1,7x15\mid x_1,7\mid x_1

3x2,7x23\mid x_2,7\mid x_2

3x3,5x33\mid x_3,5\mid x_3

如果我们找到了符合条件的 x1,x2,x3x_1,x_2,x_3,那么把他们加起来就是原方程的一个解了(注意这个解不一定是最小的)。

先考虑第一个条件,如何在 5,7=35\operatorname{5,7}=35 的倍数中,找到一个数 x1x_1 使得,x12(mod3)x_1\equiv 2\pmod{3} 呢?

求出余数为 22 的如果不好做,那不妨转化为求出余数为 11 的。因为考虑到前置知识的第一点,我们只需要找到一个 x1x_1,使得 x11(mod3)x_1\equiv 1\pmod{3},那么我只需要把 x1x_1 乘上 22 就可以了。而求出一个 x1x_1,使得 x11(mod3)x_1\equiv 1\pmod{3},怎么求呢?很明显这个式子十分类似于扩展欧几里得求逆元的式子,只需要简单转变一下即可。

对于其他的情况,用类似这样的情况分析即可。最终答案为 xi\sum x_i

扩展中国剩余定理

还是前面那个方程组,考虑 p1,p2,,pnp_1,p_2,\cdots,p_n 不一定两两互质的情况。

这时候假设我们单独把前两个方程提出来:

{xa1(modp1)xa2(modp2)\begin{cases} x\equiv a_1\pmod{p_1} \\ x\equiv a_2\pmod{p_2} \end{cases}

可以简单转化成(其中 cc 是参数,因为方便所以把 pp 换成了 aa):

{x=a1×c1+b1x=a2×c2+b2\begin{cases} x=a_1\times c_1+b_1 \\ x=a_2\times c_2+b_2 \end{cases}

简单转化一下得到:

a1×c1a2×c2=b2b1a_1\times c_1-a_2\times c_2=b_2-b_1

如果把后面的 b2b1b_2-b_1 看作一个整体,这不就是 ax+by=cax+by=c 的形式吗?因此可以将整个不定方程直接使用扩展欧几里得求出来。我们可以只关心一个未知数的值,比如 xx,在这里也就是 c1c_1

得到 c1c_1 后,带入到 x=a1×c1+b1x=a_1\times c_1+b_1 里面,就可以求解出这个 xx,显然对于第二个方程这个 xx 也是成立的。

但是如何保证这个 xx 和后面的方程都满足呢?

假设我们现在求出的这个解叫做 xx',探究 xx' 和后面若干个方程的解的关系。这个时候再把第三个方程 x3b3(moda3)x_3\equiv b_3\pmod{a_3}。如果前面的 xx' 还可以带入到某一个方程,就可以与这个联立,进一步求解新的 xx 了。不难得出前两个方程合并之后可以得到:

xx(modlcm(a1,a2))x\equiv x'\pmod{\operatorname{lcm}(a_1,a_2)}

此时再与第三个方程联立,就可以又用扩展欧几里得求出这一个新的 xx,再将新的方程与第四个联立……不断循环下去,最终得到的 xx 就是一个符合条件的解了。

扩展中国剩余定理的代码如下,因为一个神奇的性质:a×b=gcd(a,b)×lcm(a,b)a\times b=\gcd(a,b)\times\operatorname{lcm}(a,b),因此求两个数最小公倍数也很简单。

inline int gcd(int a,int b){ // 辗转相除法
	if(!b) return a;
	return gcd(b,a%b);
}
inline int lcm(int a,int b){ // 求最小公倍数
	return a/gcd(a,b)*b;
}

int x1,x2,res;
inline int exgcd(int a,int b){ // 扩展欧几里得
    if(!b){
        x1=1,x2=0;
        return a;
    }
    res=exgcd(b,a%b);
    int x3=x1;
    x1=x2,x2=x3-x2*(a/b);
    return res;
}

inline int china(){ // 扩展中国剩余定理
	int a1=a[1],b1=b[1];
	for(register int i=2;i<=n;++i){
        int a2=a[i],b2=b[i],c=b2-b1;
        exgcd(a1,a2);

        int mod=a2/res;
        x1=((x1*c/res)%mod+mod)%mod;

        int LCM=lcm(a1,a2);
        b1=((a1*x1+b1)%LCM+LCM)%LCM;
        a1=LCM;
    }
    return b1%a1;
}

事实上,还可以证明这样找出的 xx 是符合同余方程组的最小非负整数解。而且显然,扩展中国剩余定理的思路可以运用到中国剩余定理上(即 p1,p2,,pnp_1,p_2,\cdots,p_n 两两互质的时候,扩展中国剩余定理仍然正确)。

BSGS 算法

考虑一个这样的问题:给定三个数 p,a,bp,a,b,其中满足 pp 是质数。求:

axb(modp)a^x\equiv b\pmod{p}

最小非负整数解 xx

求法也并不复杂,我们将 xx 拆分成 i×mji\times m-j 的形式(其中 i,ji,j 是参数,m=pm=\left\lceil\sqrt{p}\right\rceil )。也就是我们要求出:

ai×mjb(modp)a^{i\times m-j}\equiv b\pmod{p}

移项,得:

ai×mb×aj(modp)a^{i\times m}\equiv b\times a^j\pmod{p}

然后考虑一个比较奇怪的做法:开一个哈希表,从 0m0\sim m 枚举 jj,然后将 (b×aj)modp(b\times a^j)\bmod{p} 的值当做下标,jj 的值当做数值放进哈希表里面,随后从 1m1\sim m 枚举 ii,计算出 ai×mmodpa^{i\times m}\bmod{p} 的值。

如果这个值在哈希表里面存在,就取出哈希表里面的数值 jj,这时的 i×mji\times m-j 就是一个解。

因为 jj 不管怎么变,总有 jmj\leqslant m,因此整体上解的大小随 ii 的增大而增大,所以从小到达枚举的 ii 中,枚举出的第一个解,就是原方程的最小非负整数解。

不难看出这种算法并不复杂。如果每个数的快速幂可以 O(nlogn)\operatorname{O}(n\log n) 预处理出来再做,那么每一次询问的复杂度是 O(p)\operatorname{O}(\sqrt{p}) 的。如果快速幂不能预处理(值域过大等)只能当场算,设值域大小为 nn,就是 O(plogn)\operatorname{O}(\sqrt{p}\log n) 的。如果在此基础上并不打算写哈希表,用 STL 的 map (或 set 等)求的话,是 O(plog(n×p))\operatorname{O}(\sqrt{p}\log (n\times \sqrt{p})) 的。

不过对于中等数据来说这些差别算不上很大,复杂度大实现方法实现起来也更为简单,具体还是取决于个人习惯,以及题目的数据限制卡不卡。

map 实现哈希,而且现场用快速幂的代码如下:

// mp 是哈希数组
// mi 是快速幂,都会打吧,没必要放出来
inline int bsgs(int x,int z,int mod){ 
// 求 x 的 k 次方模 mod 等于 z 的最小的 k
    mp.clear();
    x%=mod,z%=mod; // 先模为敬
    int p=sqrt(1.0*mod)+1; // 手动 ceil

    if(!z && !x) return 1;
    if(!x && z) return -1;
    // 注意特殊情况的判断
    for(register int i=0;i<=p;++i)
        mp[z*mi(x,i,mod)%mod]=i; // 放进 hash
    x=mi(x,p,mod); // 先求一个 a^m
    for(register int i=0;i<=p;++i){
        int a=mi(x,i,mod); // (a^m)^i=a^(im)
        if(mp[a]){ // 在哈希里,表示找到了
            int b=mp[a];
            if(i*p>=b) return i*p-b;
        }
    }
    return -1; // 一个都没找到,无解
}

Lucas 定理

Lucas 定理用中文表述就是卢卡斯定理。一般用于快速求一个在模意义下的组合数。

虽然把它放在最后一个(写这篇文章时的最后一个,如果后面还有别的知识都是本篇文章的更新版本加的了),但是它可能比上面讲的很多知识都简单许多。确切地说,卢卡斯定理就是:

pp 为质数,有(这里为了简化式子就将斜杠当做整除了):

CnmCn/pm/p×Cnmodpmmodp(modp)\operatorname{C}_n^m\equiv \operatorname{C}_{n/p}^{m/p}\times \operatorname{C}_{n\bmod{p}}^{m\bmod{p}}\pmod{p}

其中有除法的那一个,可以递归实现,另一个直接套组合数公式求即可。(很简单吧)

代码实现如下:

/// fac[i] 表示 i 的阶乘
// inv[i] 表示 fac[i] 的逆元
inline int C(int n,int m,int p){ // 求组合数,套公式
    if(m>n) return 0;
    return (fac[n]*inv[m]%p*inv[n-m]%p)%p;
}
inline int lucas(int n,int m,int p){ // 卢卡斯定理
    if(!m) return 1; // 注意边界
    return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p; // 套公式
}