同余代数

Luogu & cnblogs

数论中的同余代数和数论函数是数论在 OI 中的主要应用,本篇笔记总结了同余代数相关知识点。

尽管近年来的重大考试当中少有涉及同余代数的题目,但熟练掌握相关的各种算法和它们的推导过程对我们的思维水平有很大的锻炼。

数论和抽象代数的联系是非常密切的,以至于同余代数几乎完全属于抽象代数的范畴。如果能够从抽象代数(尤其是循环群)的角度自顶向下地学习同余代数(在 OI 中),那么它将非常简单且容易理解。例如:

  • 循环群在同余代数中的最直接应用就是模 p 下的加法(见 1.10 小节)。
  • 中国剩余定理的另一种形式 Z/nZZ/n1Z×Z/n2Z××Z/nkZ(见 3.3 小节)。
  • 考察 (Z/nZ)× 的结构性质给出了原根的存在性定理和在模意义下开 N 次根的算法(见 6.3 小节)。

打星号的小节(不含例题)需要一些抽象代数的前置知识,见 OI 的群论部分或抽象代数课程笔记。链接附在了参考资料部分。

定义与记号

Abstractness is the price of generality.

为了更好地理解同余理论相关算法,读者需熟知基本概念,如素数的定义,算数基本定理,同余符号及其含义,最大公约数等,并掌握基本算法如快速幂,辗转相除法求最大公约数(Euclid 算法)等。

说明

在模 1 下讨论整数同余毫无意义,下文所有同余式的模数默认 p>1,有 a1a(modp)

定义

  • 整除:若非零整数 a 是整数 b 的因数,则称 a 整除 bba 整除,记作 ab。如 248221064
  • 同余:若整数 a,b 模正整数 p 的余数相等,则称 a,b 在模 p 下同余,记作 ab(modp)。如 1557(mod7)17(mod8)
  • 最大公约数:整数 ab最大公约数 为最大的整数 d 满足 d 整除 ab,记作 gcd(a,b)。如 gcd(4,6)=2gcd(0,n)=n
  • 互质:若整数 a,b 满足 gcd(a,b)=1,则称 a,b 互质,记作 ab。一般 a,b 均为非负整数。如 381n4⊥̸6
  • 剩余类:模 p 同余的所有数构成的等价类被称为模 p剩余类。模 p 下它们等价。
  • 完全剩余系:{0,1,,p1} 构成模 p完全剩余系,记为 Zp
  • 简化剩余系:所有小于 n 且与 n 互质的正整数构成模 n简化剩余系,又称 既约剩余系缩系,记为 (Z/nZ)。这个集合的大小为 φ(n)

注意区分因数和因子:因数表示一个数的约数,因子表示乘积的一项。

记号

  • 素数集:记 P 表示素数集。P={2,3,5,7,11,13,17,19,23,29,}
  • 质因数次数:n 所含质因数 p 的次数记作 νp(n)pνp(n)npνp(n)+1n
  • 各位数字和:np 进制下的各位数字之和记作 sp(n)

基本知识(非常重要)

  • 给定 a,b,使得 abx 的最小的 xagcd(a,b),且所有这样的 x 恰为 agcd(a,b) 的倍数。对每个质因数单独分析即得。特别地,若 bx0(moda)ab,则 x0(moda)
  • dm,在 1m 中使得 gcd(m,x)=dx 的个数等于在 1md 中使得 mdxx 的个数,即 φ(md)。因为 gcd(m,x)=d 等价于 gcd(md,xd)=1

1. 基础知识

OI 中的同余理论主要是 在模意义下解方程

1.1 同余的性质

了解同余的性质是学习同余数论的基础。同余符号的根本含义:ab(modp) 等价于 a=b+kp (kZ),前者是后者的简要写法。

在同余式中,加法和乘法的交换律与结合律以及乘法分配律依然成立。以下两条性质是进行模意义下计算的依据。

  • 两整数相加(减)再对 p 取模等于两个数先对 p 取模再相加(减),再对 p 取模。
  • 两整数相乘再对 p 取模等于两个数先对 p 取模再相乘,再对 p 取模。

证明

a=q1p+r1b=q2p+r2,其中所有数都是整数且 0r1,r2<n,即 r1=amodpr2=modp。因为 kpmodp=0,所以

(a+b)modp=(q1p+r1+q2p+r2)modp=((amodp)+(bmodp))modp

乘法同理。

和等式一样,同余式的基本性质:若 ab(modp),则 a+cb+c(modp)acbc(modp)

需要注意,即使 acbc(modp),也不一定有 ab(modp)cp 是命题成立(c 有乘法逆元)的充要条件。关于同余式的除法,见 1.3 小节。

1.2 Fermat 小定理

要求模数 p 为质数。

引理

p 是质数时,其因数只有 1p。因此,若两个数相乘是 p 的倍数,则其中至少有一个是 p 的倍数。

这个引理很重要

a 不是 p 的倍数时,不存在 xy1x,y<p 满足 xaya(modp),否则 xyp 的倍数,与 1x,y<p 的限制矛盾。

考虑 1p1 所有数,它们乘以 a 之后在模 p 下互不相同且不为 0,所以仍得 1p1 所有数。于是

i=1p1ii=1p1aiap1i=1p1i(modp)(ap11)i=1p1i0(modp)

因为 i=1p1i 不是 p 的倍数,所以

ap11(modp)

该结论称为 Fermat 小定理。根据推导过程,它适用于 p 是质数且 a 不是 p 的倍数的情形。

1.3 乘法逆元

定义

为了支持同余式的除法运算,需要引入乘法逆元的概念:对整数 a,若 ax1(modp),则称 xa 在模 p 下的 乘法逆元,记作 (a1)p。模 p 同余式中出现乘法逆元,除非在符号的右下角标记模数,否则认为逆元在模 p 下,可直接写 a1

乘法逆元不一定存在,也不一定唯一。在实际运算中,只需要任意一个。

考虑同余式 axb(modp),若希望将同余式两侧同时除以 a,则只需乘以 a 的乘法逆元:

aa1xba1(modp)xba1(modp)

求法

根据 Fermat 小定理,当 p 为质数时,aap21(modp),即

a1ap2(modp)

据此快速幂求出一个数在模质数下的乘法逆元。

p 不是质数时,a 存在乘法逆元当且仅当 ap。我们将在第二章展开讨论。

以下给出几个模 质数 p 下求乘法逆元的常见技巧。

前缀逆元

增量法,设 1i1 的逆元已知。

p=ki+r (0r<i),即 kr 分别是 pi 做带余除法的商和余数。ki+r0(modp),两边同时除以 ir,得

i1kr1pi(pmodi)1

时间复杂度线性。模板题

离线 O(1) 逆元

求任意 n 个数 ai (1ai<p) 的逆元。

sa 的前缀积。通过 Fermat 小定理算出 sn1 后,计算 si1si+11×ai+1 得到前缀积的逆元,则 ai1si1×si1

时间复杂度 O(n+logp)模板题

在线 O(1) 逆元

首先 O(k) 预处理一段前缀 [1,k] 的逆元。对询问的 x,找到 u 使得 xumodpk,则 x1=u(xumodp)1

不能对每个 x 都预处理对应的 u,所以我们只能对一些 x 预处理。考虑将 1p1 分成若干组,每组的 u 可共用或根据组代表元的 u 快速计算。

最简单的分组方式是整除。设阈值 Bx=qB+r,则 xuqBu+ru(modp),考虑让 ru 较小且 qBumodp 较小。设 u 的阈值为 Bu,则 ru<BBu。令 qBumodpru 差不多,则时间复杂度至少为 Ω(BBu+pB)。取 B=Bu=p13 较优秀。

目标:对每个 q 求出 0<up13 使得 qBumodp=v,满足 |v|p23。抽屉原理保证解存在:对 qBu1v1qBu2v2(modp)(u1u2,v1v2) 即一组解。

枚举每个 u,考虑 q0 增加到 pB,每次增加 1v 增加 Bup23。若当前 q 对应的 v 不合法,因为合法 v 的区间长度大于 p23,所以可立即找到下个合法 q

对每个 u,共跳过 O(p23Bup) 次(qBu 的最大值除以 p),每次跳过带来 O(p23Bu) 个合法 q(合法区间长度除以步长),共有 O(p13) 个合法的 q

复杂度 O(n+p23)模板题

1.4 Lagrange 定理

我们知道代数学基本定理:n 次多项式在复数域 C 上有且仅有 n 个可重根。那么,模意义下的多项式是否仍然满足该性质呢?

Lagrange 定理:对于 Zp (pP) 上的整系数多项式 f(x)=i=0naixi,其中 pan,则 f(x)0(modp) 在模 p至多n 个不同解。

证明

考虑数学归纳法。当 n=0 时定理成立,因为 a00(modp) 无解(注意 pa0)。

n>0 时,假设方程有 n+1 个不同解 x0,x1,,xn。因为 f(x0)0(modp),所以 f(x)(xx0)g(x)(modp),其中 g(x) 是不超过 n1 次的多项式。

因为 xix0 (1in)(xix0)g(xi)0(modp),故 x1,x2,,xn 均为 g(x) 的根,矛盾。

解的数量并非恰为 n,而是不超过 n,例如 x220(mod3) 无解,这是 Lagrange 定理和代数学基本定理的区别之一。所以 f(x) 在模 p 下并不一定可以被分解为 n 个一次式的乘积。

  • 注意:p 不是质数时结论不成立,如 x21(mod8) 有解 x=1,3,5,7

1.5 Wilson 定理

一般形式

乘法逆元是相互的。若 ax 的乘法逆元,则 x 也是 a 的乘法逆元。这启发我们思考:当 p 是质数时,1p1 所有数能否两两配对互为逆元?

pPp>2,考虑 x21(modp) 有解 x=±1(Lagrange 定理保证不会有其它解存在),于是 2p2 两两配对互为逆元,即 i=2p2i1(modp)。因此 p!1(p1)1(modp)。对 p=2 成立。

考虑 p 不是质数的情况:

  • p=q2 时,考虑 q2q。它们相乘后模 p 等于 0。因此,若 2q<p,即 p 为大于 4 的完全平方数时,(p1)!0(modp)
  • p=4 时,3!2(mod4)
  • p 为大于 4 的非完全平方数时,令 qp 的最小质因数,则 qpq,所以 (p1)!0(modp).

综上,我们得到 Wilson 定理(p1)!1(modp) 当且仅当 p 是质数。

扩展形式

尝试将 Wilson 定理扩展至素数幂的情形。

考虑 pk 以内所有与 p 互质的数的乘积模 pk,记为 (pk!)p

仍然考虑求解 x21(modpk),求出所有逆元为本身的数,并将不为 x 的其它所有数与其逆元两两配对。因为我们只考虑与 p 互质的数,所以逆元存在。因此只需考虑所有 x 的乘积。

因式分解得 (x+1)(x1)0(modpk)

p>2 时,x+1x1 不可能均是 p 的倍数,必然有 x+10(modpk)x10(modpk)。故仍有 (pk!)p1(modpk)

p=2 时,有平凡解 ±1。除此以外,x+1x1 可能同时为 2 的倍数。但它们不可能同时为 4 的倍数。若 x+1x1 同时为 2 的倍数,那么将 (x+1)(x1)0(mod2k) 两侧除掉不能被 4 整除的数(模数除以 2,因为被除掉的数恰好贡献一个质因数 2),得到 x±10(mod2k1)。当 p=2 时,方程有四个解 ±12k1±1

注意 k=1 时四个解均重合,此时 1!1(mod2)k=2 时两对解重合,此时 1×31(mod4)。否则 1×(1)×(2k1+1)×(2k11)1(mod2k)

综上,得到 Wilson 定理的扩展形式

(pk!)p{1p=2k31otherwise(modpk)

在 exLucas 中用到了该结论。

1.6 Kummer 定理

阶乘的素数幂次

给定正整数 n 和质数 p,求 νp(n!)

提出 1n 当中所有 p 的倍数,将它们全部除以 p。这一步消掉的质因数个数为 np。上一步操作后,得到 1np。将所有 p 的倍数提出来,然后全部除以 p。这一步消掉的质因数个数为 np2。以此类推,直到 npk<p。此时已经没有 p 的倍数了。有如下结论:

νp(n!)=i=1logpnnpi

该结论由 Legendre 在 1808 年提出。

考虑换一种方法求和。尝试对 np 进制下的每一位,求出它对每次提出因数个数的贡献之和。令 c=logpn,设 n=i=0caipi (0ai<p),则

νp(n!)=i=0cj=1caipipj=i=0caij=0i1pj=i=0caipi1p1=(i=0caipi)(i=0cai)p1=nsp(n)p1

其中用到了等比数列求和公式 i=0kpi=pk1p1

通过上述推导,得到 νp(n!) 的另一种表示方法:n 减去 np 进制下的各位数字和,再除以 p1。特别地,当 p=2 时,ν2(n!)=npopcount(n)

组合数的素数幂次

根据组合数公式 (nm)=n!m!(nm)!,有

νp((nm))=nsp(n)(msp(m))(nmsp(nm))p1=sp(m)+sp(nm)sp(n)p1

考虑 nmmp 进制下相加得到 n 的过程。每产生一次进位,各位数字之和就会减少 p1,因此 p(nm) 中的幂次等于 p 进制下 nmm 需要进位的次数。

该结论称为 Kummer 定理

1.7 步长与子环

结论

在长为 n 的环上每一步走 k 条边,形成 d=gcd(n,k) 个子环,每个环的环长为 nd

按顺序给环上的每个点标号 0,1,2,,n1。从 r 开始每一步走 k 条边,经过的点的编号与 r 在模 d 下同余,因为 n,k 均为 d 的倍数,r 加上 kn 取模不改变 rmodd

引理

nk 时,从环上任意一点出发,能够走遍所有点。

证明

即证 (r+ck)modn0c<n 互不相同。假设存在 0i<j<nr+ikr+jk(modn),移项得 (ji)k0(modn)。因为 kn,所以 nji,这与 0i,j<n 矛盾。

说明

这个引理本质上和最开始的基本知识等价。如果跳 x 步回到原来的位置,那么 nkx,推出 xngcd(n,k) 的倍数。但 nk,所以 x 无论如何都是 n 的倍数。

即当 kn 时,ckmodn 取遍 0n1,这是 Fermat 小定理的引理的推广,可以进一步推出 Euler 定理。

不妨设 0r<d。考虑 0n1 当中所有模 dr 的数 r,r+d,,r+(nd1)d,将它们重标号为 0,1,,nd1,形成长为 nd 的子环。

因为 dn,k 的最大公约数,所以 ndkd。根据引理,在这个长为 nd 的子环上每一步走 kd 条边(相当于在原环上走 k 条边),能够走遍子环上的每一个点。

于是子环长为 ndngcd(n,k)。每个模 d 的剩余系形成一个长为 nd 的子环,子环个数为 d=gcd(n,k)

例题:P5582P5887P6187

1.8 Euler 定理

前置知识:Euler 函数(见数论 II)。

一般形式

回顾 Fermat 小定理的证明过程,我们发现 p 是质数的条件只是为了保证 axmodp 互不相同。由 1.7 小节的引理,这个条件可以仅由 ap 保证。同时我们知道,两个模 p 互质的数相乘后模 p 仍与 p 互质。因此,我们尝试将费马小定理扩展至更一般的情况。

Kpp 的简化剩余系。对于任意 x,yKp,若 xy,则 axayp 的余数不相等且仍与 p 互质。因此,在模 p 下,Kp 中的每个数乘以 a 后仍与 Kp 相等。故 xKpxxKpax(modp)

等式两边同时除以 xKpx 得到 Euler 定理

aφ(p)1(modp)

在计算与模数互质的某个数的幂时,指数可以对模数的 Euler 函数取模。可以用来化简公式或减小常数,前提是模数是定值或其 Euler 函数容易计算。

p 是质数时,φ(p)=p1,与 Fermat 小定理一致。

扩展形式

gcd(a,p)1 时,不断乘以 a 会让 aip 的 gcd 不断增加,直到某个特定的 ar 之后 gcd 不再变化,因为此时 gcd(ar,p) 的每个质因数受到的限制都是 p 提供的。将 gcd 除掉之后得到 Euler 定理的条件,于是有 exEuler 定理

ab{ab mod φ(p)gcd(a,p)=1abgcd(a,p)1, b<φ(p)ab mod φ(p)+φ(p)gcd(a,p)1, bφ(p)(modp)

证明

d=gcd(ar,p)。在模 p 下考虑 ar,ar+1,,它们均为 d 的倍数。因为 apd,所以它们除以 d 之后会取遍所有与 pd 互质的数。根据 Euler 定理,akmodpd 有循环节 φ(pd),因此 (ar+kmodpd)×dar+kmodp 也有长度为 φ(pd) 的循环节。而 φ(pd)φ(p),所以从 ar 开始,对 k0ar+kmodp 有长度为 φ(p) 的循环节。

还需证明 rφ(p)。感性理解就是如果 gcd 变化,一定会乘上一个因子,所以 rO(logp) 级别的,会比 φ(p) 小。具体怎么证明呢?我们发现 r 就是 a 的每个质因数幂次 qiki 对应的 ri 的最大值,所以只需对 a 是质数幂次 qk 的情况证明。设 pc 个质因数 q,则 r=ck。因此,显然地,我们只需对 a 是质数 q 的情况证明,因为此时 r 最大。

p=qcp,则 r=c。由 Euler 函数的性质,φ(p)φ(qc)=qc1(q1)c

1.9 同余式的除法

对于 ab(modp),一种除法是将 a,b 同时除以 d=gcd(a,p) 以化简同余式(如果 b 无法除尽,则同余式不可能成立),此时 p 也应当除以 gcd(a,p)。另一种除法是在 ap 时,将同余式两侧同时除以 a,通过乘以 a 的乘法逆元实现,此时 p 保持不变。

读者需要区分这两种对同余式的基本变形方法。同余式的本质为存在 k 使得 a+kp=b,第一种除法相当于将等式两侧同时除以 d 得到 ad+kpd=bd,因为等式左侧是整数,所以右侧也必须是整数,即 db。第二种除法相当于 ax+kpx=bx,其中 ax1(modp),于是存在 k 使得 1+kp=bxmodp,即 1bx(modp),这是对同余式的简单变形。

*1.10 循环群和直积

循环群 是只由一个元素生成的群 G=x={xnnZ}。大小相同的有限循环群是同构的。大小为 p 的循环群(p 阶循环群)记为 Zp,则 xp=1(这里的 1 是单位元)。

Zp 的直观理解是一个 p 个点的环 01p10,编号为 i 的点对应群内的元素 xi,于是元素之间的乘法运算(定义一个群时的二元运算通常看成元素之间的乘法运算,但注意这并不是整数乘法,尽管它们的运算规则几乎一致)可以看成在环上移动。乘以 xi 相当于在环上跳 i 步。为什么要想象成一个环呢?因为这样就可以用 1.7 小节的结论得到一些不那么显然的结论。

Zp 在同余代数中的最直接应用就是模 p 下的代数加法。因为 xp=1,所以 Zp 上的运算关于 x 的指数是在模 p 下的。这也是为什么模 p 的完全剩余系记作 Zp={0,1,,p1}:用 i 代替 xi,元素的乘法定义为模 p 下的加法。

特别地,当同时需要模 p 下的加法和乘法时,会把这个代数结构写成两个环(环同时定义了乘法和加法,其中加法构成交换群)的商 Z/pZ(Z/pZ)× 表示所有和 p 互质的数模 p 构成的乘法群。

定义两个群 G,H直积G×H={(x,y)xG, yH},其上的二元运算定义为 (x1,y1)(x2,y2)=(x1x2,y1y2)。简单地说,就是用一个二元组把 GH 绑起来,但是 GHG×H 上的运算仍保持一定的独立性。

2. 二元线性不定方程

扩展 Euclid 算法(exgcd)用于求解形如 ax+by=c二元线性不定方程,其中 a,b,c,x,y 都是整数,x,y 是未知数。其得名于和求最大公约数的 Euclid 算法相似的流程。

2.1 Bezout 定理

在求解 ax+by=c 之前,我们需要先判定它是否有解。

首先考虑一些较弱的结论(必要条件)。无论 x,y 的取值如何,左式一定是 d=gcd(a,b) 的倍数。因此,当 dc 时,方程无解。于是,可以将方程两侧同时除以 d,此时 ab

根据直观感受,如果 ab,那么 ax+by 能组合成任何整数。只需证明 axmodb 可以取到 [0,b1] 的所有整数。相当于从 0 开始在长为 b 的环上每次跳 a 步。因为 ab,根据 1.7 小节的引理,可以经过环上的所有点。

因此,对任意整数 c,存在 x 使得 axc(modb)。取 y=caxb 即可。

综上,我们得到了 Bezout 定理:二元整数线性不定方程 ax+by=c 有解当且仅当 gcd(a,b)c

2.2 扩展 Euclid 算法(exgcd)

欲求解 ax+by=c,设 d=gcd(a,b),只需求解 ax+by=d,并将解乘以 cd。根据 Bezout 定理,方程有解。

注意到等式右端等于 x,y 的系数的最大公约数。回顾 Euclid 算法,其使用了关键结论 gcd(a,b)=gcd(b,amodb) 以减小问题规模。利用该思想,尝试将方程也缩至更小的规模上。

因为 d=gcd(a,b)=gcd(b,amodb),所以 bx+(amodb)y=d 有解。解这个方程得到 (x,y),则

ax+by=d=bx+(amodb)y=bx+(ab×ab)y=ay+b(xaby)

对比等式两侧,令 x=yy=xaby 即得一组 (x,y)。不断递归缩小问题规模直到平凡情况。当 b=0 时,根据辗转相除法,a=d,此时 x=1y=0 即为一组解。

扩展 Euclid 算法的时间复杂度是辗转相除法的 O(logV)

void exgcd(int a, int b, int &x, int &y) {
  if(!b) return x = 1, y = 0, void();
  exgcd(b, a % b, x, y);
  int _y = x - (a / b) * y;
  x = y, y = _y;
}

也可以写成

void exgcd(int a, int b, int &x, int &y) {
  if(!b) return x = 1, y = 0, void();
  exgcd(b, a % b, y, x), y -= a / b * x;
}
  • 注意:exgcd 求得的解为 ax+by=d 的一组特解。为得到原方程 ax+by=c 的特解,还需将 x,y 乘以 cd

2.3 扩展

通解形式

ax+by=c 有解,则显然有无穷多组解。我们需要通解的形式。

使用扩欧求得特解 (x0,y0)。对任意一组解 (x,y)ax+by=ax0+by0,即 (axax0)+(byby0)=0。根据 Δ(ax)+Δ(by)=0,设 Δ=|Δ(ax)|=|Δ(by)|,则 a,bΔ,即 lcm(a,b)Δ。因此 Δx 是必须是 lcm(a,b)a=bd 的倍数,Δy 同理。这是必要条件。

x 增加 bdy 减少 ad,方程仍成立,因此条件充分,即 Δxbd 的倍数是通解的充要条件。故通解形如

{x=x0+bdky=y0adk(k\Z)

特解的数值范围

关于 exgcd 求得特解的数值范围,详见 博客。结论如下:对于非平凡情况(a,b0ab),有 |x||b2d||y||a2d|

应用

  • 解一元线性同余方程 axb(modp)。看成二元线性不定方程 ax+py=b 使用 exgcd 求解。根据 Bezout 定理,方程有解当且仅当 gcd(a,p)b

  • 当模数不是质数时,Fermat 小定理不再适用,但模意义下的乘法逆元仍可能存在。a 在模 p 下存在逆元当且仅当 ap,因为求逆元相当于解 ax1(modp)

  • axb(modp)gcd(a,p)=db 时无解,否则在 0p1 中有 d 个整数解。这是因为其通解可写为 x0+kpd 的形式,其中 0x0<pd:将方程两侧同时除以 d,则 adx0bd(modpd)

2.4 例题

P5656【模板】二元一次不定方程 (exgcd)

题目要求 x 是正整数,先求出 x 的最小正整数值 x0=xmodbd,此时 ymax=cax0b。若 ymax>0 则有正整数解,且 ymax 对应正整数解中最大的 y,否则没有正整数解。对于 y0xmax 同理。

#include <bits/stdc++.h>
using namespace std;
#define int long long
int T, a, b, c;
void exgcd(int a, int b, int &x, int &y) {
	if(!b) return x = 1, y = 0, void();
	exgcd(b, a % b, y, x), y -= a / b * x;
}
signed main() {
	cin >> T;
	while(T--) {
		scanf("%lld %lld %lld", &a, &b, &c);
		int d = __gcd(a, b), x, y, xx, xy, yx, yy;
		if(c % d) {
      puts("-1"); continue;
    }
		exgcd(a, b, x, y);
		x *= c / d, y *= c / d;
		xx = x % (b / d);
		if(xx <= 0) xx += b / d;
		xy = (c - a * xx) / b;
		yy = y % (a / d);
		if(yy <= 0) yy += a / d;
		yx = (c - b * yy) / a;
		if(xy <= 0) cout << xx << " " << yy << "\n";
		else cout << (yx - xx) / (b / d) + 1 << " " << xx << " " << yy << " " << yx << " " << xy << "\n";
	}
	return 0;
}

UVA12775 Gift Dilemma

cgcd(a,b,c)200 说明 c200。枚举 c 前的系数 z,再求 ax+by=pcz 的非负整数解的个数。

z=pc,则时间复杂度为 O(z),因为只需调用一次 exgcd。

CF1728E Red-Black Pepper

先通过一遍贪心对每个 i 求出选 i 个红辣椒的最大值。

对每组询问,exgcd 求出使得 x,y 均非负的 x 的最小整数解和最大整数解,若不存在则无解,否则三分或在极值点两侧找到最近的两个解。

时间复杂度 O((n+q)logn)代码

[模拟赛] 你还没有 AK 吗

给定 n,m,求整数 x[0,n]y[0,m],执行任意正整数次 xx+yyx+y 后使得 x=X 的初始值 x,y 的个数。

多组数据。T105n,m,X1018

f1=f2=1fi=fi1+fi2,则相当于求多少对 (x,y) 满足 f2k+1x+f2k+2y=X

Fibonacci 数列增长很快,考虑直接枚举 k。问题转化为求关于 x,y 的二元一次方程 ax+by=X 有多少组解满足 x[0,n]y[0,m](Fibonacci 数列相邻两项互质)。使用 exgcd 求解,但 O(Tlog2X) 无法通过。

不难发现只有本质不同的 ka,b,直接预处理它们 exgcd 的结果,复杂度 O(TlogX)

进一步地,有 f2k+121(modf2k+2),所以 xf2k+1 即可求得一组特解,不需要使用 exgcd,但需要 __int128

*P3518 [POI2011] SEJ-Strongbox

v 是密码,则所有 n 且是 gcd(v,n) 的倍数的数也是密码,因为 kvmodn 取到了所有这样的数。

证明

d=gcd(v,n)dtt 不是密码,则 kvt(modn) 无解。根据 Bezout 定理,它等价于 gcd(v,n)tdt,矛盾。

进一步地,若 u,v 是密码,则 u=gcd(u,n)v=gcd(v,n) 是密码。由裴蜀定理,gcd(u,v) 在模 n 下能被 ux+vy 表出,所以 gcd(u,v) 是密码。

因此,设密码集合为 S,则 d=gcd(n,gcduSu)S。显然,S 恰由 d 的所有倍数组成。

考虑枚举这个 d=gcd(v,n),若合法则答案即 maxnd,即我们需要找到最小的合法的 d

di=gcd(mi,n)d 必须是密码与 ngcddk 的因数,其次任何 di (1i<k) 不能是 d 的倍数。对于后者的限制,相当于在 n 的所有因数形成的图上,一个点向它的因数连边,能被某个 di 到达的因数是不合法的。

首先给所有 di (1ik) 打上标记。从大到小枚举 n 的每个因数 c,若 c 被打上标记,则 cpj 也应被打上标记,其中 pj 表示能整除 cn 的质因数,表示若 c 是某个 di 的因数,则 cpj 也是。

剩下没有被打标记的 n 的因数 d,若 d 能被 dk 整除则合法。找到最小的这样的 d,则答案为 nd

打标记的过程可以使用哈希表实现,时间复杂度 O(n+klogn+d(n)ω(n)),其中 ω(i) 表示 i 的质因数个数。代码

*P3421 [POI2005] SKO-Knights

题目要求我们找到两个向量 d,e,使得它们的 整系数 线性组合得到的线性空间(张成)等于给出的所有向量的整系数张成。

span(S) 表示集合 S 的整系数张成,即

span(S)={aiSciai:ciZ}

Sol 1

如果我们能将三个向量 a,b,c 合并成等价的两个向量 d,e,就能解决本题。根据基础的线性代数知识,不改变张成的初等行变换有三种:

  • 交换:对应本题即任意交换 a,b,c
  • 倍乘:对应本题即乘以任意非 0 整数。
  • 倍加:对应本题即向量加上任意整数倍其它向量。

当然,上述性质仅在 实系数 线性组合的前提下成立。当要求系数为 整数 时:

  • 交换显然不改变张成。
  • 倍乘变换 可能改变 张成,因为整数乘法不可逆。如 span{(0,1),(1,0)} 显然不等于 span{(0,2),(1,0)}
  • 倍加变换 不改变 张成。对任意 p=y(a+xb)+zbspan(a+xb,b),它仍是 a,b 的整系数线性组合 ya+(xy+z)b。对任意 p=ya+zbspan(a,b),我们也总可以将其表示为 a+xbb 的整系数线性组合 y(a+xb)+(zxy)b

很好!如果整系数倍加变换不改变整系数张成,就可以使用 辗转相减法

具体地,考虑两个向量 a,b,我们能够在不改变其整系数张成的前提下,将 b 的某一维变为 0。只需对 a,b 在这一维进行辗转相减法即可。

考虑三个向量 a,b,c,我们对 a,bx 这一维进行辗转相减,再对 a,cx 这一维进行辗转相减。此时 xb,xc 均为 0,意味着 b,c 共线。因此,对 b,cy 这一维进行辗转相减,即 ybgcd(yb,yc),此时 c 变成零向量,对 span(a,b,c) 没有贡献,可以直接舍去,即此时 span(a,b,c)=span(a,b)

综上,我们在 O(nlogV) 的时间内解决了问题,其中 V 是值域。代码

注意题目限制坐标绝对值不超过 104。辗转相减法得到最终的 a 时,xayb 是原坐标的 gcd,显然满足坐标限制。但 ya 不一定,因此要对 yb 取模,相当于做了一步辗转相减(此时 xb=0 所以不需要改变 xa)。最终得到的坐标绝对值不超过 102,比题目限制优一个数量级。

从上述过程中,我们发现一个有趣的性质:整系数下,两个向量可以在不改变其张成的前提下,使其中一个向量的某一维变成 0,而另一个向量的对应维变成原来两个向量在这一维上的 gcd。这也是解决本题的关键性质。

Sol 2

让我们深入地钻研一下题解区的其它做法。

实际上两者的核心思想是等价的:将两个向量的其中一个的某一维变成 0,另一个的对应维变成 gcd。不同点在于,线性代数做法是从倍加变换和辗转相除法推得该性质,即我们使用正确性有保证的方法,最终得到的结果是这样的性质。而题解区的做法是首先令这一条件成立,再根据得到的线性方程组,使用 exgcd 和一些数学推导求解。这体现了 exgcd 与辗转相除的本质联系。

不妨设对于向量 a,b,我们要将 xb 变为 0。则 xa 只能等于 dx=gcd(xa,xb)(不管 y 坐标,由 Bezout 定理,能达到的 x 坐标恰为 dx 的所有倍数)。因此,对于方程 xax+xby=dx,使用 exgcd 求得一组特解 (x0,y0),则 ya=yax0+yby0,因为 a=x0a+y0b

对于 b,观察到 a 的通解为 (x0+kxbdx,y0kxadx),所以保持 xa=dx 不变,ya 的最小变化量为 yaxbdxybxadx,即让 k1ya 产生的变化量。因为 xb=0,所以 yb 只能等于这个最小变化量。可以验证这样得到的 a,b 符合要求。

综上,我们将 a,b 写成了 a,b,其中 xa=gcd(xa,xb)ya 用 exgcd 求解,xb=0yb=yaxbxaybgcd(xa,xb)。对 a,cy 坐标做类似的操作,再合并 b,c 即可。

Sol 3

另一种求解 yb 的思路(来源 TheLostWeak):a,b 可以整系数线性组合出 a,b

考虑 a。因为 xa=gcd(xa,xb),所以为了使横坐标等于 xa,需要将 a 乘以 xagcd(xa,xb)xaxa 的系数,得到 (xa,xayaxa)。因为可以通过向量 b=(0,yb) 得到 (xa,ya),因此 yb yaxayaxa 。同理 yb ybxbyaxa。因此

yb=gcd(yaxaxaya,ybxaxbya)xa

yb 小于上述 gcd,将导致存在 pspan(a,b)pspan(a,b)。若 yb 大于上述 gcda,b 无法整系数线性组合出 a,b 中的至少一个。

抛砖引玉:联立后两种方法描述 yb 的式子,稍作化简后得到等式

yaxbxayb=gcd(yaxaxaya,ybxaxbya)

CF516E Drazil and His Happy Friends

不妨设 nm

如果一个男生被一个女生变得快乐,那么要么这个女生一开始就是快乐的,要么这个女生先被别的男生变得快乐。对于前者,总共只有 O(g) 种情况。对于后者,对应天数可以表示为某个男生变得快乐的时间加上若干倍的 n。因此,设 fi 表示第 i 个男生变得快乐的时间,那么 f(i+kn)modm 可以被 fi+kn 更新。

i(i+n)modm 连边成环,用 exgcd 找到初始 fi 有值(fi<n)的男生在环上的位置。考虑每相邻两个这样的男生之间的贡献,即相邻两个这样的男生 ij 中,j 的前驱的 f 对答案产生贡献。

男生的 f 的最大值和女生的 f 的最大值的较大值即为答案。

n,m 不互质时,对每个模 d=gcd(n,m) 的同余类分别求解即可。

时间复杂度 O((b+g)log(b+g))代码

*P3543 [POI2012] WYR-Leveling Ground

区间操作转化为差分数组 di=aiai1 (1in+1) 的端点操作,其中 an+1=0。一次 +a 和一次 a 抵消,一次 +b 和一次 b 抵消。设 A=adB=bd

首先考虑对每个差分值单独求解,解不定方程 ax+by=di。设 d=gcd(a,b),根据 Bezout 定理,若 ddi 则无解。否则用 exgcd 求得 ax+by=d 的特解,再乘以 did 得到 ax+by=di 的特解。

先不考虑可行性。因为 axi+byi=di 的贡献为 |xi|+|yi|,而最终答案为所有贡献之和除以 2,所以先找到 |xi|+|yi| 最小的特解。若 xiyi 均不为最小非负整数或最大负整数可行解,则调整法可证 |xi|+|yi| 必然不是最优。因此我们只需检查 xiyi 分别取到最小非负整数和最大负整数的情况。

当前 12(|xi|+|yi|) 是答案的下界,但不能保证可行性,还需满足 X=xi=0。当 X=0Y=0,所以只需考虑 X

根据特解的形式

{x=x+kBy=ykA

X<0 时,我们需要进行 XB 次将某个 xi 加上 B,并将对应的 yi 减去 AXB 一定是整数,因为差分数组之和为 0。类似的,X>0xi 减去 Byi 加上 A。称这样的操作为对 i 进行一次调整,调整的代价等于新的 |xi|+|yi| 减去原来 |xi|+|yi|

容易证明调整一个数的代价随着调整次数增加仅在 xiyi 改变符号时增加,其它时候不变,因此我们可以用堆维护这个过程:每次取出代价最小的 i 并弹出,在需求次数与不改变符号的限制下尽可能多地调整。若经过调整后需求次数为 0,则算法结束,否则将新的调整代价压入堆。时间复杂度 O(nlogn),因为 xi,yi 中任意一个改变符号才会改变代价,最多发生 2n 次。

一个神奇的性质是 |XB|O(n) 级别,这使得我们可以从堆中取数时仅进行一次调整就塞回去。证明(来自 评论区):当 |xi|+|yi| 取到最小值时若 xi 为最小非负整数或最大负整数,则 |xi|B,则 |X|nB。若 yi 为最小非负整数或最大负整数则 |yi|A,即 |XA+diB|nA,根据 di=0|X|nB

忽略问题限制得到基本解法再调整,巧妙结合反悔贪心和 exgcd。代码

*P7325 [WC2021] 斐波那契

解决本题需要的知识:对任意整数 m,普通 Fibonacci 数列在模 m 下纯循环,且循环节长度为 O(m)。此外,fp1fp

F0=aF1=b 的第 pFp=afp1+bfp,独立 a,b 前的系数可得该结果。

首先特判 a,b=0 的情况。令 b=b,则问题等价于

afp1bfp(modm)

自然,我们希望将 b 移到方程左侧,fp1 移到方程右侧,这样可以预处理所有 fpfp1 的值。但是 m 并不是质数,因此需要化简。

  • 化简的核心思路是,同余方程 AB(modm) 等价于 AdBd(modmd),其中 A,B 均为代数式,dA,B,m。直观理解很显然,感性说明就是 d 在该同余方程中完全没有起到任何作用。

d=gcd(a,b,m),等式两侧同时除以 d,得到

afp1bfp(modm)

其中 a=adb=bdm=md

欲将等式两侧同时除以 b,即保证 b 在模 m 下有逆元,则需 bm。但并不一定满足。

d=gcd(b,m)。根据既约性(gcd(a,b,d)=1),da,因此必然有 dfp1,因为方程右侧是 d 的倍数。因此,

afp1dbfp(modm)

其中 b=bdm=md。此时可以将 b 除过去。

a(b)m1fp1dfp(modm)

欲将等式两侧同时除以 fp1d,则需 fp1dm,因为 Fibonacci 数列的性质 fpfp1,而等式右侧是 gcd(fp1d,m) 的倍数,所以 gcd 只能等于 1

以上推导过程说明我们需要枚举 d,dp,并将 (d,d,fp(fp1d)m1modm) 作为三元组塞入 map。根据一开始的结论,复杂度是 m 的所有因数的因数和,再乘以 maplog,无法接受。此时需要利用使得方程有解时这三个变量必须满足的性质降低复杂度。

注意到 dfp1fp1dmd,这说明枚举 d,p 之后,d 只能等于 gcd(fp1,md)。时间复杂度变成 O(σ(m))m 的约数和,大约是 mloglogm 量级(证明见 ycx's blog),总复杂度 O((mloglogm+n)logm)代码

注意:当 fpfp1 等于 0 时,直接忽略,因为 fpfp1 不可能同时为 0a,b=0 的情况特判)。

总结一下,我们在化简过程中,假设了 dd 两个仅与 a,b,m 有关的变量,因此预处理时需要枚举 dd。但根据方程有解时 dd 所具有的性质,我们得以降低复杂度。这种思想方法(在回答多组询问时,为了解决问题,需要设和询问参数相关的变量,那么在预处理时枚举这些变量的所有情况,即可考虑到给定询问参数的所有情况)是很巧妙的。

3. 线性同余方程组

前置知识:扩展 Euclid 算法。

形如

{xa1(modm1)xa2(modm2)xak(modmk)

的方程组称为 线性同余方程组,其应用非常广泛(从三模 NTT 到高次剩余)。最常见的用处是将模任意合数简化为模指数幂。特别地,当模数不含平方因数时(square-free number,μ(p)0),可以简化为模质数。

exCRT 比 CRT 更简单且适用范围更广泛,因此笔者在合并线性同余方程组时一般使用 exCRT。

3.1 中国剩余定理(CRT)

中国剩余定理(Chinese Remainder Theorem,CRT)用于求解 模数两两互质 的线性同余方程组。对任意 ij,均有 mimj

CRT 的思想是求出若干个数,使得这些数依次满足每个同余方程,且在其它模数下均为 0,则这些数的和即为所求。可以理解为 ZM 上的 Lagrange 插值,其中 M=mi

我们尝试为每个同余方程构造 xiai(modmi),记为条件一,且对于 ji 均有 xi0(modmj),记为条件二。

为满足条件二,xi 必须是 jimjci=Mmi 的倍数,令 xi=aici

为满足条件一,需要给 xi 乘以 ci 在模 mi 下的逆元。因为 cimi 互质,所以 di=(ci1)mi 存在。注意 mi 不一定是质数,所以需要 exgcd 求逆元。于是 xi=aicidi

对每个模数进行上述流程,则 (i=1kaicidi)modM 即为所求。时间复杂度线性对数。模板题 P1495 代码,注意使用 __int128

#include <bits/stdc++.h>
using namespace std;
constexpr int N = 10 + 5;
int n, a[N], b[N];
long long M = 1, ans;
void exgcd(int a, int b, int &x, int &y) {
	if(!b) return x = 1, y = 0, void();
	exgcd(b, a % b, y, x), y -= a / b * x;
}
int inv(int x, int p) {
  return exgcd(x, p, x, *new int), (x % p + p) % p;
}
int main() {
	cin >> n;
	for(int i = 1; i <= n; i++) cin >> a[i] >> b[i], M *= a[i];
	for(int i = 1; i <= n; i++) ans = (ans + (__int128)b[i] * (M / a[i]) * inv(M / a[i] % a[i], a[i])) % M;
	cout << ans << endl;
	return 0;
}

3.2 扩展中国剩余定理(exCRT)

扩展中国剩余定理(exCRT)用于求解一般形式的线性同余方程组。除了求解的问题相似,它和 CRT 并没有太大联系。

当模数不互质时,考虑合并两个同余方程 xa1(modm1)xa2(modm2)。因为 x 可以表示成 a1+pm1,所以 a1+pm1a2(modm2),即 pm1+qm2a2a1(modm2),exgcd 求解 p 即可。也可以理解为 pa2a1m1(modm2) 并使用 exgcd 求逆元。合并后的方程为

xa1+pm1(modlcm(m1,m2))

根据 Bezout 定理,若 gcd(m1,m2)a2a1,则方程组无解。时间复杂度 O(klogV)。模板题 P4777 代码。

3#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll n, a, b, c, d;
void exgcd(ll a, ll b, ll &x, ll &y) {
  if(!b) return x = 1, y = 0, void();
  exgcd(b, a % b, y, x), y -= a / b * x;
}
int main() {
  cin >> n >> a >> b;
  for(int i = 1; i < n; i++) {
    cin >> c >> d;
    ll e = __gcd(a, c), f = (d - b % c + c) % c, inv;
    c /= e, f /= e;
    exgcd(a / e, c, inv, *new ll), inv = (inv % c + c) % c;
    b += (__int128)f * inv % c * a, a *= c, b %= a;
  }
  cout << b << endl;
  return 0;
}

注意考虑 exgcd 解出来是负数的情况。在得到解之后要进行取模。

*3.3 抽象代数中的 CRT

当且仅当 nm 时,Z/nmZ(Z/nZ)×(Z/mZ)。这个表达式的含义是,存在二元组 (x1,x2)x 的一一对应,其中 x1[0,n1]x2[0,m1]x[0,nm1],使得对应关系保持结构。设 f(x1,x2) 表示对应的 x,那么

f(x1,x2)+f(y1,y2)=f(x1+y1,x2+y2),f(x1,x2)f(y1,y2)=f(x1y1,x2y2)

实际上,一个双射由 x(xmodn,xmodm) 给出。

需要指出的是,抽象代数中的 CRT 并没有给出根据 (x1,x2)x 的方法,且其一般形式要复杂一些。

3.4 例题

P4774 [NOI2018] 屠龙勇士

设第 i 次的攻击力为 ci,multiset 模拟确定。找到最小的 x 使得对所有 i 都有 cixai(modpi)cixpi

先不管第二个条件,设 di=gcd(ci,pi),若 diai 则无解,否则将 ai,ci,pi 同时除以 di,此时 (ci,pi)=1,把 ci 除过去就是标准的线性同余方程组。exCRT 求解得到 x

为满足第二个条件,先求出 ei=aici 表示打败 i 至少需要多少次攻击。若最终求得 x<maxei,那么将 x 加上 maxeixp×p,其中 p 是最终的同余方程的模数。

时间复杂度 O(nlogp)代码

P4621 [COCI2012-2013#6] BAKTERIJE

细菌的运动情况会在不超过 4n2=104 次后循环。先 BFS 找环,特判前 104 秒,然后 4k 枚举每个细菌进入 (x,y) 时的方向,从而唯一确定一个线性方程组,exCRT 即可。

时间复杂度 O(4kklogV)代码

P2480 [SDOI2010] 古代猪文

通过组合数学可知题目要求 gdn(nd)mod999911659。模数是质数,由 Fermat 小定理,指数对 999911658 取模。分解质因数后发现是 square-free number 且每个质因数较小。根据 CRT 拆成组合数对小素数取模,使用 Lucas 定理计算。

时间复杂度 O(d(n)logn)代码

4. 组合数取模

前置知识:组合数,二项式定理(组合数学)。

4.1 Lucas 定理

Lucas 定理是连接组合数学和数论的桥梁。

根据 1.6 小节的结论,当 p 是质数时,(nm)modp=0 当且仅当 nmp 进制下需要借位,这等价于 np 进制下的每一位均不小于 mp 进制下的对应位。

i<j 时,(ij) 定义为 0。这启发我们思考:如果将 n,m 表示为 p 进制数 n=i=0caipim=i=0cbipi(不足的位补 0),那么是否有 (nm)i=0c(aibi)(modp)

引理

p 是质数时

(x+y)pxp+yp(modp)

证明

使用二项式定理将上式拆开,当 1i<p 时,(pi)p 的倍数。

推论

(1+x)p1+xp(modp)

Lucas 定理

p 是质数时,

(nm)(npmp)(nmodpmmodp)(modp)

证明

n=dp+r (0r<p),则

(1+x)n= (1+x)dp(1+x)r (1+xp)d(1+x)r(modp)

由二项式定理,上式 xm 前的系数即 (nm)

因为 (1+xp)d 展开后 x 的次数均为 p 的倍数,(1+x)r 展开后 x 的次数小于 p,因此 xm 项只能由 (1+xp)d 展开后的 xm/pp 项与 (1+x)r 展开后的 xmmodp 项相乘得到。

另一种形式

n,mp 进制下为 i=0caipii=0cbipi,则

(nm)i=0c(aibi)(modp)

显然,两种形式是等价的。

根据 Lucas 定理,计算大组合数对 p 取模可拆成计算小组合数。O(p) 预处理阶乘及其逆元做到 O(1) 计算组合数,时间复杂度 O(logpn)。模板题 P3807 代码。

#include <bits/stdc++.h>
using namespace std;
constexpr int N = 1e5 + 5;
int T, n, m, p, fc[N], ifc[N];
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % p;
    a = 1ll * a * a % p, b >>= 1;
  }
  return s;
}
int binom(int n, int m) {
  return 1ll * fc[n] * ifc[m] % p * ifc[n - m] % p;
}
int lucas(int n, int m) {
  if(n % p < m % p) return 0;
  if(n < p) return binom(n, m);
  return 1ll * lucas(n / p, m / p) * binom(n % p, m % p) % p;
}
void solve() {
  cin >> n >> m >> p;
  for(int i = fc[0] = 1; i < p; i++) fc[i] = 1ll * fc[i - 1] * i % p;
  ifc[p - 1] = ksm(fc[p - 1], p - 2);
  for(int i = p - 2; ~i; i--) ifc[i] = 1ll * ifc[i + 1] * (i + 1) % p;
  cout << lucas(n + m, n) << endl;
}
int main() {
  cin >> T;
  while(T--) solve();
  return 0;
}

p 进制下考虑组合数取模是很有帮助的。计算大组合数对 p 取模的值,可以将 n,mp 进制下的每一位分开考虑,再将所得结果相乘。特别地,若 np 进制下至少有一位小于 m,则 (nm)0(modp),因为对 ai<bi(aibi)=0。这符合 1.6 小节的结论。

Lucas 定理还可以快速计算组合数奇偶性。(nm) 是奇数当且仅当 m 在二进制下每一位均不大于 n,当且仅当 nm 的按位与等于 m,当且仅当 nm 的按位或等于 n

4.2 扩展 Lucas 定理(exLucas)

当模数 p 不是质数时,也有快速计算 (nm)modp 的方法。

模合数的问题的第一步一定是 CRT,将问题简化为模指数幂的情况。具体地,设 p=qiki,分别计算答案模 qiki 的结果,再使用 CRT 合并。这就是 CRT 的神力!

考虑 (nm)modqk,自然地把组合数写成 n!m!(nm)!。因为 m!(nm)! 含有质因数 q,所以无法求逆元。但因为模数只含质因数 q,所以我们考虑把 q 和其它质因数分开来算:求出组合数含有多少 q,以及分子和分母除掉其所包含的所有 q 之后的结果,即变形为

n!qvq(n!)m!qvq(m!)(nm)!qvq((nm)!)×qvq(n!)vq(m!)vq((nm)!)

使用 Kummer 定理计算组合数的质因数幂次(vq(n!)vq(m!)vq((nm)!))。此时分子和分母均不含 q,分别计算,再对分母求逆元即可。问题转化为 n!qvq(n!)modqk,也就是阶乘去掉所有 q 之后对 qk 取模。

d=nqdk=nqkr=nmodqk。既然除掉了所有 q,自然想到将所有数分成 q 的倍数和不是 q 的倍数讨论。

对于 q 的倍数,有 q,2q,,dq,我们给每个数除掉 q,问题变成求 d!qvq(d!)。这是子问题。

对于不是 q 的倍数,把它们全部模掉 qk,得到 dk1qk 内和 q 互质的数,以及一组 1r 内和 q 互质的数。回忆扩展 Wilson 定理的符号,(n!)q 表示 1n 所有和 q 互质的数的乘积,所以只需计算 ((qk!)qdk×(r!)q)modqk。预处理 qk 以内所有与 q 互质的数的前缀积,这部分可以 O(logn) 计算(快速幂)。进一步地,根据扩展 Wilson 定理,(qk!)qmodqk 等于 11,所以不需要快速幂。

综上,

n!qvq(n!)d!qvq(d!)×(qk!)qdk×(r!)q(modpk)

n<q 时,直接返回预处理结果,否则使用递归式计算。时间复杂度 O(logqn)

  • 注意:递归的子问题是除以 q,而循环节数量是除以 qk
  • 注意:区分 n!qvq(n!)(n!)q。前者只是将所有 p 去掉,而后者将含有质因数 p 的所有数去掉了。

综上,我们在 O(ω(p)logn+qiki) 的时间内求出一个大组合数模任意数 p 的值。若模数固定且预处理,则可以做到单次 O(ω(p)logn)。模板题 P4720 代码。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
void exgcd(int a, int b, int &x, int &y) {
  if(!b) return x = 1, y = 0, void();
  exgcd(b, a % b, y, x), y -= a / b * x;
}
int inv(int x, int p) {
  return exgcd(x, p, x, *new int), (x % p + p) % p;
}
ll n, m;
int p, tp, ans, curp = 1, q, qk, fc[1000005];
ll num(ll n) {
  ll s = 0;
  while(n) s += n /= q;
  return s;
}
int calc(ll n) {
  if(!n) return 1;
  return 1ll * calc(n / q) * (tp && (n / qk & 1) ? qk - 1 : 1) % qk * fc[n % qk] % qk;
} // 注意 n / q 和 n / qk!
int solve() {
  tp = (qk & 1) || qk <= 4; // 扩展 Wilson 定理.
  for(int i = fc[0] = 1; i < qk; i++) fc[i] = 1ll * fc[i - 1] * (i % q ? i : 1) % qk;
  int pw = num(n) - num(m) - num(n - m);
  int res = 1ll * calc(n) * inv(1ll * calc(m) * calc(n - m) % qk, qk) % qk;
  while(pw--) res = 1ll * res * q % qk; // 根据 Kummer 定理, 质因子在组合数中的个数为对数级别.
  return res;
}
int main() {
  cin >> n >> m >> p;
  for(int i = 2; i <= p; i++) {
    if(i * i > p) i = p;
    if(p % i == 0) {
      q = i, qk = 1;
      while(p % i == 0) qk *= i, p /= i;
      ans += 1ll * (solve() - ans % qk + qk) * inv(curp, qk) % qk * curp;
      ans %= (curp *= qk); // exCRT 合并更方便.
    }
  }
  cout << ans << endl; 
  return 0;
}

4.3 例题

[BZOJ2445] 最大团

定义一张 n 个点的图是 “好图” 当且仅当它能被分成若干大小相等的完全图。结点有标号。记 n 个点的好图个数为 G(n)。求 mG(n)mod999999599n,m2×109

模数是质数,Fermat 小定理转化为求 G(n)mod(109402)

被分成的完全图大小 dn 的因数,考虑枚举 dn 并计算完全图大小为 d 时的方案数。

根据组合意义,从 n 个节点中选出 d 个点连成完全图,再从剩下 nd 个节点中选出 d 个点连成完全图,以此类推。最后剩下的 d 个节点连成完全图。由于完全图之间无序,而我们枚举完全图的过程有序,因此要除以 nd 的阶乘。记 c=nd,则

ans=1c!(nd,d,,d)=n!(d!)c×c!

观察 99999598,将其分解质因数后得到 2×13×5281×7283,模数是 square-free number。考虑分别求上式对这四个质数取模后的结果,再 CRT 合并。问题转化为计算答案对质数 p 取模后的结果,其中 p 较小。

根据 1.6 小节的结论,先求出答案含有多少个质因子 p,若个数不为 0 则直接返回 0。否则问题转化为计算一个数的阶乘去掉所有质因子 p 后模 p 的结果,使用 exLucas 即可。代码

*P5598 【XR-4】混乱度

S(l,r) 表示 a 的区间和,根据组合意义得到

f(l,r)=f(l,r1)(S(l,r)ar)=f(l,r)(S(l,r)S(l,r1))

固定 l,考虑 r 在增大的过程中,S(l,r) 不断增加。根据 Lucas 定理,一旦 S(l,r)p 进制下进位,那么之后的 f(l,r)modp 均为 0。因此,在经过至多 plogpV 个非零的 ar 后,f(l,r) 对答案不再有贡献。

预处理每个非零的 arp 进制下有哪些位非零以及这一位上的值,对每个 l,维护当前 S(l,r)p 进制下每一位的值,可快速计算 (S(l,r)S(l,r1))。考虑到 ar=0f(l,r1)=f(l,r),所以对每个位置预处理下一个非零的位置即可快速跳过等于 0 的位置。

时间复杂度 O(nplogpn)代码

5. 离散对数问题

离散对数问题即求解 离散对数方程

axb(modp)

其中 x 即模 p 下的 logab。解可能不存在,也可能存在多个,我们只要求出任意一个。

5.1 大步小步算法(BSGS)

大步小步算法(Baby Step Giant Step,BSGS)要求 ap

BSGS 使用了 meet-in-the-middle:设块长为 Mx=AMB,则 aAMbaB(modp)。因为 ap,所以若 aAMbaB(modp),则 aAMBb(modp)。枚举 aAMbaB 所有可能的取值,并通过哈希表判断是否存在 A,B 使得这两个值相等。时间复杂度 O(max(A,B))

根据 Euler 定理,若有解,则在 [0,p2] 上有解,所以 1BM1Ap1M。将阈值 M 设为 p 得到最优复杂度 O(p)。模板题 P3846 代码。

#include <bits/stdc++.h>
using namespace std;
int BSGS(int a, int b, int p) {
  int A = 1, B = sqrt(p) + 1;
  a %= p, b %= p;
  unordered_map<int, int> mp;
  for(int i = 1; i <= B; i++) {
    A = 1ll * A * a % p;
    mp[1ll * A * b % p] = i;
  }
  for(int i = 1, AA = A; i <= B; i++) {
    if(mp.find(AA) != mp.end()) {
      return i * B - mp[AA];
    }
    AA = 1ll * AA * A % p;
  }
  return -1;
}
int main() {
  int p, b, n;
  cin >> p >> b >> n;
  int ans = BSGS(b,n,p);
  if(~ans) cout << ans << "\n";
  else puts("no solution");
  return 0;
}

根据枚举顺序可知以上代码求得 logab 的最小的非负整数解。解的循环节长度为最小的正整数 y 使得 ay1(modp),即 δp(a),见 6.1 小节。

5.2 扩展大步小步算法(exBSGS)

对于更一般的离散对数方程,没有 ap 的条件,该怎么办呢?

从已知到未知。既然可以解决 ap,那么我们尝试把 a,p 凑成互质。同余式两侧同时除以 d=gcd(a,p),方程化为 adax1bd(modpd)。若 db 显然无解。

注意到 adpd,所以前者在模后者下存在逆元。p 可能不是质数,需要使用 exgcd 求逆元。移项,方程变为 ax1bd×(ad)1(modpd),等式右边为常数。

因为 apgcd(a,p) 可能不互质,所以重复上述操作直到 a,p 互质,此时直接使用 BSGS 求解即可。设提出的 a 的数量为 k,则求出的解需要加上 k。注意在每次操作开始前特判 b=1,此时答案恰等于已经提出的 a 的数量。

显然 k=O(logp),时间复杂度 O(p)。模板题 SP3105 代码。

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
using namespace std;
using namespace __gnu_pbds;
int a, p, b;
void exgcd(int a, int b, int &x, int &y) {
	if(!b) return x = 1, y = 0, void();
	exgcd(b, a % b, y, x), y -= a / b * x;
}
int inv(int x, int p) {return exgcd(x, p, x, *new int), (x % p + p) % p;}
int BSGS(int a, int b, int p) {
	a %= p, b %= p;
	int A = 1, B = sqrt(p) + 1;
	gp_hash_table <int, int> mp;
	for(int i = 1; i <= B; i++) mp[1ll * (A = 1ll * A * a % p) * b % p] = i;
	for(int i = 1, AA = A; i <= B; i++, AA = 1ll * AA * A % p) if(mp.find(AA) != mp.end()) return i * B - mp[AA];
	return -1;
}
int exBSGS(int a, int b, int p) {
	int d = __gcd(a, p), cnt = 0;
	a %= p, b %= p;
	while(d > 1) {
		if(b == 1) return cnt;
		if(b % d) return -1;
		p /= d, b = 1ll * b / d * inv(a / d, p) % p;
		d = __gcd(p, a %= p), cnt++;
	}
	int ans = BSGS(a, b, p);
	return ~ans ? ans + cnt : -1;
}
int main() {
	cin >> a >> p >> b;
	while(a) {
		int ans = exBSGS(a, b, p);
		if(~ans) cout << ans << "\n";
		else puts("No Solution");
		cin >> a >> p >> b;
	}
	return 0;
}

5.3 例题

P2485 [SDOI2011] 计算器

操作 1 快速幂,操作 2 exgcd 求逆元,操作 3 BSGS。

P3306 [SDOI2013] 随机数生成器

x 的坐标左移一位,显然有 xi=aix0+bj=0i1aj。由等比数列求和公式,

xi=aix0+ai1a1×b=ai(x0+ba1)ba1

对同余式进行简单变形后得到离散对数问题,使用 BSGS 求解。时间复杂度 O(Tp)代码

6. 阶与原根

对于含乘法的同余式 abc(modm),如果能找到 g,x,y,z 使得 gxa(modm)gyb(modm)gzc(modm),那么根据 Euler 定理,同余式等价于 x+yz(modφ(m))。也就是说,我们希望找到一个共同的底,将幂次的乘法化为指数上的加法。

在模数是合数时,一般会将问题化为只考虑和模数互质的数。m (m2) 的简化剩余系在模 m 的乘法下构成群(封闭、结合律、逆元、单位元)。借用群相关的定义,得到阶和原根(循环群的生成元)的概念。

am,如果 a 的幂次模 m 可以取到整个简化剩余系,那么称 am 的原根。可以发现 a 就是上述问题的一个较为可行的底,因为所有和 m 互质的数都可以用它的幂次表示。

6.1 阶及其性质

定义

ax1(modm) 的最小正整数 x 称为 am,记作 δm(a)amδm(a) 存在的充要条件。

必要性

am 不互质,则 axm 不互质,模 m 不可能得 1

充分性

am,根据 Euler 定理,x=φ(m) 即一解,因此阶存在。

以下设 am

性质

阶本质上在刻画 a 的幂次的循环节,a 的幂次模 m 是纯循环。如果 axay(modm),则 ayx1(modm)。由此我们知道 a0,a1,,aδm(a)1 在模 m 下一定是互不相同的,否则会和 δm(a) 的最小性矛盾。于是 a 的幂次模 m 一定是两两不同的 a0,a1,,aδm(a)1 的不断循环。

更具体地,axay(modm) 当且仅当 xy(modδm(a))

性质 1

ax1(modm) 当且仅当 δm(a)x

证明

ax1(modm),则 akx1(modm)。因此若 δm(a)x,则 ax1(modm)

δm(a)x,设 r=xmodδm(a),则 axar(modm)。因为 0<r<δm(a),所以 ar1(modm)

考虑一个数自身相乘时,它的阶如何变化。如果我们将模 ma 的幂次看成一个长为 δm(a) 的环,那么 a 的幂次相当于在环上每次跳一步,ak 的幂次相当于在环上每次跳 k 步。使用 1.7 小节的结论即可。

性质 2

δm(ak)=δm(a)gcd(δm(a),k)=lcm(δm(a),k)k

证明

注意到 (ak)x=akx

对于第一部分,根据性质 1,δm(ak) 为最小的 x 使得 δm(a)kx。于是 xmin=δm(a)gcd(δm(a),k)(见最开头的基本知识)。

对于第二部分,同时是 kδm(a) 的倍数的最小正整数为 lcm(δm(a),k),所以使得 akx1(modm) 的最小的 kxlcm(δm(a),k)。于是 xmin=lcm(δm(a),k)k

以上两部分也可以根据 ab=gcd(a,b)×lcm(a,b) 相互推导。

考虑两个数相乘时,它们的阶如何变化。和一个数自身的幂次 δm(ak)δm(a) 不一样,乘积的阶可能小于其中每个因子的阶,也可能大于。我们需要知道这种变化是由什么引起的。从循环节的角度考虑会带给我们一些启发:乘积的阶变小是由于其两个因子的阶的公共因数,乘积的阶变大是由于其两个因子的阶的非公共因数。

如果 δm(a)=δm(b),那么 δm(ab) 可以是 δm(a) 的任何因数。δm(g)=δm(g2)=δm(g4)=δm(g14)=15,而 δm(gg14)=1δm(gg4)=3δm(gg2)=5δm(gg)=15

性质 3

δm(ab)=δm(a)δm(b) 当且仅当 δm(a)δm(b)

证明

必要性是显然的。

充分性的证明思路是注意到对所有 δm(a) 的倍数或 δm(b) 的倍数但不为 lcm(δm(a),δm(b))=δm(a)δm(b) 的倍数 x,对于 axmodmbxmodm,一定恰有一个等于 1,另一个不等于 1,所以 x 一定不是 δm(ab) 的倍数。

两个最小循环节分别为 δm(a)δm(b) 的序列相乘,一定有 L=δm(a)δm(b) 的循环节,最小循环节显然是 L 的因数。如果最小循环节不等于 L,那么存在 L 的质因数 p 使得 (ab)L/p1(modm)。不妨设 pδm(b),则 Lp=kδm(a),其中 k=δm(b)p。因为 δm(a)δm(b),所以 kδm(a) 不是 δm(b) 的倍数。但 1(ab)kδm(a)bkδm(a)(modm),矛盾。这说明最小循环节等于 L

对于更一般的情况呢?读者可以先结合以上两种特殊情况自行思考。

性质 4

对于质因数 p,如果它在 δm(a)δm(b) 的次数相等,那么它在 δm(ab) 的次数是任意的,但不会超过它在 δm(a) 的次数。如果不等,那么它在 δm(ab) 的次数是它在 δm(a)δm(b) 的次数的较大值。

可以这样类比:设 pa 的次数为 x,在 b 的次数为 y。如果 x=y,那么 pa+b 的次数 zx 可以任意大。如果 xy,那么 z=min(x,y)。至于为什么是加法,考虑到阶本身是描述指数循环节的,而相乘转化到指数上就是加法。

如果 m 有原根 g,那么一切都是好理解的:δm(ab)=φ(m)gcd(φ(m),logga+loggb),指数在分母上解释了一个是任意小一个是任意大,一个是较大值一个是较小值。

性质 4 描述了乘积的阶与每个因子的阶的关系。

性质 2 和性质 4 给予我们操作阶的空间,并引出了以下性质。

性质 5

给定 a,bm,存在 c 使得 δm(c)=lcm(δm(a),δm(b))

证明

对每个质因数 p,如果 pδm(a),δm(b) 的次数相等且不为 0,那么让 aap。由性质 2,因为 pδm(a),所以 δm(ap)=δm(a)p。处理后 δm(a)δm(b) 没有次数相同的质因数。根据性质 4 可知 δm(ab)=lcm(δm(a),δm(b))

离散对数的通解形式:当 am 时,使得 axb(modm)x 若存在,则其循环节恰为 δm(a)。可以求出最小和次小的解 x1,x2,则通解形式为 x=x1+k(x2x1) (k\N)。也可以先把阶求出来。

求法

法一:使用 BSGS 求 ax1(modm) 的最小正整数解。

法二:根据性质 1,对 x=kδm(a)ax1(modm)。考虑令 x=φ(m),对于 φ(m) 的每个质因子 p,用 x 不断试除 p 直到 x 无法整除 pax/p1(modm)。时间复杂度为分解质因数加上 O(log2m)

6.2 原根

以下设 m2am

定义

δm(a)=φ(m),则称 a 为模 m原根。并不是所有数都有原根。存在原根的模 m 的缩系同构于循环群:原根等价于循环群的生成元。

我们需要原根是因为它可以 生成m 的缩系,即它的若干次幂在模 m 下取遍了所有与 m 互质的数。这使得我们在考虑解与模数互质且模数存在原根的同余方程时可以用原根的若干次幂代替未知数,在求解高次剩余时非常有用。

简单地说,它可以将乘法转化为加法。

性质

判定原根的一般方法是直接计算阶。如果 δm(a)φ(m),那么一定存在 φ(m) 的质因数 p 使得 φ(m)pδm(a) 的倍数,即 aφ(m)/p1(modm)

原根判定定理

am 的原根当且仅当对任意 φ(m) 的质因子 p,均有 aφ(m)/p1(modm)

除了判断一个数是否是原根,我们还需要判断一个数是否有原根。

原根存在定理

m 有原根当且仅当 m=2,4,pα,2pα,其中 p奇质数

证明较复杂,不感兴趣的读者可以跳过。

2,4 有原根是显然的。

奇质数

对奇质数 p,由引理 5,存在 g 使得对任意 1i<pδp(i)δp(g),所以 xδp(g)1(modp)p1 个不同的解。由 Lagrange 定理,δp(g)=p1

奇质数幂

考虑数学归纳。

g 是奇质数 p 的原根。如果 gp11(modp2),那么

(g+p)p1gp1+(p1)gp2p1pg11(modp2)

g 是奇质数幂 pα 的原根,1k<p,且 gpα1(p1)1+kpα(modpα+1)(当 α=1 时,上述说明保证这样的 g 存在),那么

gpα1(p1)1+(k+cp)pα(modpα+2)

于是

gpα(p1)(1+(k+cp)pα)p1+kpα+1(modpα+2)

g 是奇质数幂 pα+1 的原根。由数学归纳法可知 g 是任意 pα (α1) 的原根。

  • 请思考以上推导过程在哪里对 p=2 不适用。

gpα 的原根,考虑 gmodpα 可知 δ2pα(g)δpα(g)=φ(pα)=φ(2pα)。所以 g2pα 的原根。

必要性

首先,如果 m 有原根,那么 m 的所有因数都有原根。所以只需证明 84pp1p2 没有原根。8 没有原根是显然的。

对任意 x4p,因为 φ(p)φ(4)=2 的倍数,所以 xφ(p)1(modp)xφ(p)1(mod4)。由 CRT,xφ(p)1(mod4p)

对任意 xp1p2,因为 φ(p1),φ(p2) 都是 2 的倍数,所以 xφ(p1p2)/21(modp1)xφ(p1p2)/21(modp2)。由 CRT,xφ(p1p2)/21(modp1p2)

如果一个数有原根,那么其缩系上的乘法等价于模 φ(m) 下的加法:将缩系的每个数变成原根的幂次,两个数相乘是同底数幂相乘,指数相加。这提供了更多阶的性质。

性质 6

m 有原根,则使得 δm(a)=la 的个数为

{0lφ(m)φ(l)lφ(m)

证明

根据性质 1,第一部分显然。考虑第二部分。

gm 的原根。因为任意 am 均可表示为 gk (0k<φ(m)),所以答案即使得 δm(gk)=lk 的个数。根据性质 2,δm(gk)=δm(g)gcd(δm(g),k),因此 φ(m)gcd(φ(m),k)=l,即 gcd(φ(m),k)=φ(m)l。由基本知识,k 的个数为 φ(φ(m)φ(m)/l)=φ(l)

形象地,将 gkmodmgk+1modm 看成长度为 φ(m) 的环。在环上每次走 k 步形成子环,则 gk 的阶等于子环环长。由 1.7 小节,环长为 ngcd(n,k)。为使 φ(m)gcd(φ(m),k)=l,则 gcd(φ(m),k)=φ(m)l

性质 6 反映了欧拉函数的性质 dnφ(d)=nlφ(m)φ(l)=φ(m),即模 m 的阶等于 l 的数的个数之和(模 m 下存在阶的数的个数)等于和 m 互质的数的个数,而和 m 互质是在模 m 下存在阶的充要条件。

原根个数定理

若正整数 m 有原根,则原根个数为 φ(φ(m))

求法

求一个数的所有原根:求出任意一个原根 g,根据性质 2,gkmodm (kφ(m)) 也是 m 的原根。相当于在长为 φ(m) 的环上每次走 k 步,因为 kφ(m),所以一定走遍整个环。

中国数学家王元证明了质数的最小原根的数量级为 O(m0.25+ϵ)(1959,需要指出的是大 O 符号仅表示上界)。结合原根存在定理的推导过程,可以证明任何有原根的数最小原根不超过这个数量级。再根据原根判定定理,可以 O(m) 预处理每个数的最小质因子后 O(m0.25ω(m)logm) 求出 m 的最小原根,从而求出 m 的所有原根。模板题 P6091 代码。

#include <bits/stdc++.h>
using namespace std;
constexpr int N = 1e6 + 5;
int ksm(int a, int b, int p) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % p;
    a = 1ll * a * a % p, b >>= 1;
  }
  return s;
}
void solve() {
  int n, d;
  cin >> n >> d;
  int x = n, phi = n;
  for(int i = 2; i <= x; i++) {
    if(x % i == 0) {
      phi = phi / i * (i - 1);
      while(x % i == 0) x /= i;
    }
  }
  vector<int> pr;
  x = phi;
  for(int i = 2; i <= x; i++) {
    if(x % i == 0) {
      pr.push_back(i);
      while(x % i == 0) x /= i;
    }
  }
  for(int i = 1; i < n; i++) {
    if(__gcd(i, n) != 1) continue;
    bool ok = 1;
    for(int it : pr) {
      ok &= ksm(i, phi / it, n) != 1;
    }
    if(ok) {
      static int is[N], g[N];
      memset(is, 0, N << 2);
      is[i] = 1;
      for(int j = 2; j <= phi; j++) {
        if(__gcd(j, phi) == 1) {
          is[ksm(i, j, n)] = 1;
        }
      }
      int cnt = 0;
      for(int j = 1; j <= n; j++) {
        if(is[j]) g[++cnt] = j;
      }
      cout << cnt << "\n";
      for(int p = 1; p <= cnt / d; p++) {
        cout << g[p * d] << " ";
      }
      cout << "\n";
      return;
    }
  }
  cout << "0\n\n";
}
int main() {
  int T;
  cin >> T;
  while(T--) solve();
  return 0;
}

补充:因为 nφ(n)O(loglogn) 级别,所以 nφ(φ(n))=O((loglogn)2)。如果每次在 1m1 中均匀随机选择一个数,那么期望随机 O((loglogm)2) 次得到一个原根,比 O(m0.25) 更好。

原根和单位根

m 有原根 g,设 n=φ(m),则 m 的简化剩余系与 n 次单位根 同构

对于任何 am,其可以写成 gk。任何一个 n 次单位根也可以写成 ωnk。同时 gn1(modp)omegann=1。因此,我们可以将单位圆应用在 m 的简化剩余系上,为原根提供了一个形象的表示方法,如下图。

qlTzex.png

图中,单位根的 09 次幂等分单位圆。ωni×ωnj=ωn(i+j)mod10 对应 gi×gjg(i+j)mod10(mod11),环上运算(两个数相乘)的循环往复(逆时针旋转)体现了 “模” 的操作。

因此,对于特定模数,我们可以使用原根代替 n 次单位根进行快速傅里叶变换 FFT。当需要使用 2k 次单位根时,若 2kφ(p)p 有原根,可使用原根代替 2k 次单位根,即快速数论变换 NTT。常见模数 998244353 的 Euler 函数 φ(998244353)=222×238

类比原根和单位根,得到以下性质。

性质 7

p 是奇质数或奇质数的幂,gp 的原根,则 gφ(p)/21(modp)

证明

gx1(modp),则 xφ(p)。如果 xφ(p)2,那么 (gx)2=g2x1(modp),和 gx1(modp) 矛盾。因此 logp1=φ(p)2

性质 7 在二次剩余部分有用。

*6.3 (Z/nZ)× 的结构

关于记号,见 1.10 小节。

根据 CRT,我们只需要考虑 (Z/pkZ)× 的结构。

根据原根的定义,(Z/pkZ)× 是循环群当且仅当 pk 有原根。对奇质数 p,已知 pk 有原根,所以 (Z/pkZ)× 一定是循环群。实际上,证明 (Z/pkZ)× 和证明 pk 有原根等价。笔者在查阅资料时找到了一篇分析 (Z/nZ)× 的结构的论文,其中直接证明 (Z/pkZ)× 是循环群的推导过程和以上证明 pk 有原根是几乎一致的。

考虑 (Z/2kZ)× 的结构。当 k=1k=2 时平凡。读者可以先猜一猜结果会是什么样。由有限生成阿贝尔群基本定理,(Z/2kZ)× 同构于若干循环群的直积,如果有一个元素的阶是 2k2,那么 (Z/2kZ)× 只能是 Z2×Z2k2,因为 |(Z/2kZ)×|=2k1

事实上,5 就是这样的元素。

引理 1

k\N52k1 恰含有 k+2 个质因数 2

证明

k=0 时命题成立。52k+11=(52k1)(52k+1),前一项恰因子含有 k+2 个质因数 2,是 4 的倍数,所以后一项因子不是 4 的倍数。因此 52k+11 恰含有 k+3 个质因数 2

引理 2

k3δ2k(5)=2k2

证明

因为阶是 2k1 的因数但不能等于 2k1(没有原根),所以 δ2k(5)2k2。如果 52k31(mod2k),那么 52k31 含有至少 k 个质因数 2,和引理 1 矛盾。因此 δ2k(5)=2k2

现在我们知道 5Z2k2 这一项的生成元,那么 Z2 这一项的生成元是什么?找一个阶为 2 的元素,也就是找一个不是 1 但平方是 1 的数。我们知道这样的数有 3 个,1 是其中之一,那么 1 是否合法呢?

引理 3

k2m\N5m1(mod2k)

证明

如果成立,那么 5m1(mod4),但 5m1m1(mod4),矛盾。

上述推导说明 (Z/2kZ)×Z2×Z2k2,其中生成元分别是 15。于是得到以下结论。

性质 8

对任意 k3a2,存在 c[0,1]d[0,2k21] 使得 a(1)c5d(mod2k),且这种表示唯一,即每个 a[0,2k1] 的奇数唯一对应一组 (c,d)

性质 8 是求模 2k 下的 N 次剩余的关键结论。

现在让我们站在更高的视角看看为什么 2k (k3)4p (p2)pq (p,q2) 没有原根:由 CRT,对 nm(Z/nmZ)×(Z/nZ)××(Z/mZ)×,于是

(\Z/2k\Z)×Z2×Z2k2(\Z/4p\Z)×(\Z/4\Z)××(\Z/p\Z)×Z2×Zp1(\Z/pq\Z)×(\Z/p\Z)××(\Z/q\Z)×Zp1×Zq1

由 3.3 小节,如果 nm 不互质,那么 Zn×Zm 一定不是循环群,所以这些数对应的模意义下的乘法群都不是循环群,自然没有原根。

6.4 例题

P5605 小 A 与两位神仙

前置知识:Pollard-rho(数论 III)。

模数 p 是奇质数的幂,存在原根 g

尝试将 xy 写成 gXgY。根据欧拉定理,xay(modp) 等价于方程 XaY(modφ(p))。根据 Bezout 定理,后者有解的充要条件是 gcd(X,φ(p))Y,这等价于 gcd(X,φ(p))gcd(Y,φ(p))。根据阶的性质 δp(gk)=δp(g)gcd(δp(g),k),我们发现 δp(x)=φ(p)gcd(φ(p),X),因此条件又等价于 δp(y)δp(x)

直接求阶即可,需要使用 Pollard-rho 分解 φ(p),时间复杂度 O(nlog2p+p0.25)

本题告诉我们的结论:对于奇质数幂 pαx,ypxay(modpα) 有解当且仅当 δpα(y)δpα(x)

*P6730 [WC2020] 猜数游戏

和上一题差不多的思路。

分成 ap 互质和 ap 不互质两种情况考虑。并且将每个数的贡献拆开来考虑,即求出有多少种情况使得 a 对答案有贡献。

当且仅当 S 中不存在一个数能生成 a 时,a 对答案有贡献。因此对每个 a,答案加上 2nc,其中 c 是能生成 a 的数的个数。

ap 不互质时,因为 p 是奇质数的幂,所以在不超过 log 次后 a 的幂变为 0,因此直接考虑每个数能生成哪些数即可得到每个数能被多少个数生成。使用哈希表存储,该部分时间复杂度 nlogp

ap 互质时,直接使用上一题的结论,对每个数求阶后做高维后缀和即可。注意若两个数能互相生成,即阶相等,则需要钦定一个顺序,使得前面的数钦定生成后面的数。

7. N 次剩余问题

N 次剩余问题即求解形如 xna(modp) 的方程。注意未知数 x 在底数上,而离散对数问题的未知数在指数上。

7.1 定义与符号

pa,若存在整数 x 使得 xna(modp),则称 apn 次剩余

特别地,对 pa,若存在整数 x 使得 x2a(modp),则称 ap二次剩余,反之称 ap二次非剩余。注意,这两个定义基于 a 不是 p 的倍数,也就是说,如果 ap 的倍数,那么 a 既不是 p 的二次剩余,也不是 p 的二次非剩余。

Legendre 符号:记作 (ap)。当 ap 的二次剩余时,该值为 1;当 ap 的二次非剩余时,该值为 1;当 ap 的倍数时,该值为 0

7.2 二次剩余的分布

奇质数

p 是奇质数时,二次剩余非常容易理解。如果 ap 的二次剩余,那么 ax2(modp)。用原根表示 ax,设 agd(modp)xgk(modp),那么 gdg2k(modp),即 d2k(modp1)。因为 p1 是偶数,所以 d 必须是偶数,否则无解。当 d 是偶数时,k 有两个解(1.3 小节应用部分),每个 k 对应一个 xgk(modp)

也就是说,原根的偶数次幂和奇数次将 1p1 分成两类,偶数次幂对应二次剩余,奇数次幂对应二次非剩余。因此,二次剩余和二次非剩余的数量均为 p12。这和原根的选择无关,也就是说,任意原根的偶数次幂集合都是相同的。

进一步地,如果 x0x1 都是 a 的平方根,那么由平方差公式,(x0+x1)(x0x1)0(modp)。因为 p 是质数,所以 x0+x10(modp),即 x0x1 互为相反数。从原根的角度分析,x0x1 写成以 g 为底的幂时在指数上相差 p12,而 g(p1)/21(modp)。另一种角度是如果 x2a(modp),那么 (x)2a(modp)

类比:正数(二次剩余)的平方根有两个且互为相反数,负数(二次非剩余)没有平方根。

奇质数幂

将奇质数的结论稍微推广即可,读者可以先自行尝试推导结论。

pk 有原根 g,所以对于 ap,依然可以用奇质数的分析方法得到 logga 是偶数时是二次剩余,是奇数时是二次非剩余。同样地,其平方根有两个且互为相反数。

设非零整数 a=pca。如果 a 是二次剩余,那么存在 x=pdx 使得 ax2(modpk)。显然 c=2d 必须是偶数。此外,还要看 a 是否为 pkc 的二次剩余。如果是,设 xa(x)2(modpkc) 的解,那么 pc/2x 是原方程的解,即 a 是二次剩余。否则 a 是二次非剩余。

以上推导可以推广到 N 次剩余的情况。实际上,这是求解 N 次剩余的一部分。

2 的幂

p=2k

首先考虑 a 是奇数的情况。注意到任何奇数的平方模 81,所以 ax2(modp)x 是奇数的必要条件为 a1(mod8)。那么这个条件是否充分呢?

假设 y 是方程的解,那么 (xy)(x+y)0(modp)。因为 y 是奇数,所以 2y 不是 4 的倍数。因此,x+yxy 至少有一个不是 4 的倍数。如果 x+y 不是 4 的倍数,那么 xy 必须是 2k1 的倍数,即此时 x 恰有两个解。类似地,如果 xy 不是 4 的倍数,那么此时 x 同样恰有两个解。综合两种情况,如果 ax2(modp) 有解,那么至多有 4 个解。这个分析和扩展 Wilson 定理非常类似。

一共有 2k1 个奇数,且至多有 2k3 个奇数的平方,但一个平方数最多有 4 个根,那么唯一的情况是所有模 81 的数都是二次剩余,且它们的平方根都有 4 个。综上,如果 a 是奇数,那么 a 是二次剩余当且仅当 a1(mod8)

对于 a 是偶数的情况,和奇质数幂的推导类似。a 是二次剩余当且仅当 a=4ka,其中 a 是奇二次剩余。

任意合数

根据 CRT 拆成模质数幂。显然地,对 paap 的二次剩余当且仅当对任意拆出的质数幂 qka 不是 qk 的二次非剩余。注意并非 aqk 的二次剩余,区别在于 a 可以是 qk 的倍数。

7.3 二次剩余的判定:Euler 准则

7.2 小节给出了 2k 的二次剩余的判定准则,以及其它情况是如何归约到判定奇质数的二次剩余。因此,本节只关注奇质数的情况。

给定奇质数 p,则 p 有原根 g

对于整数 a,我们希望知道 a 在以 g 为底时的指数是不是偶数。那么有没有什么办法不用求出指数就能确定它的奇偶性呢?提示:利用阶的性质 7。

注意到偶数的 p12 倍是 p1 的倍数,但奇数的 p12 倍不是 p1 的倍数。因此,如果 k 是偶数,那么 gk(p1)/21(modp);如果 k 是奇数,那么根据阶的性质 7,gk(p1)/21(modp)。而 gk(p1)/2a(p1)/2

Euler 准则

p 是奇质数。

ap12(ap)(modp)

简单地说,a 是二次剩余还是二次非剩余等价于 ap12modp,等价于 a 能否表示成 p 的一个原根的 偶数 次幂。读者应仔细体会其中的等价关系。

Euler 准则同时告诉我们 Legendre 符号关于 a 是完全积性函数。即对于任意 a,b

(ap)(bp)=(abp)

笔者的理解:由 Fermat 小定理,对奇素数 pap11(modp)。指数除以 2 相当于开平方根,此时 1p1 会因为自身性质的不同而分化成两类,一类是二次剩余 a(p1)/21(modp),另一类是二次非剩余 a(p1)/2=1(modp)。这种分化,用单位根的知识来解释,就是对于 p1 次单位根 ωp1k,当 k 为偶数时 (ωp1k)(p1)/2=1,反之为 1。原根和单位根的联系恰好解释了 agk(modp)k 的奇偶性对 a 是否是二次剩余的影响。

7.4 模意义下开平方根

考虑求解 ax2(modp),以下只考虑 ap 的二次剩余的情况。

7.4.1 模质数下开平方根

p 是奇质数(模 2 开平方根是平凡的),解方程 ax2(modp) (a0),则 x 即模 p 下的 x

由 7.2 小节的分析,方程恰有两个解,只需求出任意一个解,通过取相反数得到另一个。

使用 BSGS,显然存在根号时间的算法:先求出 logga。因为 a 是二次剩余,所以 logga 是偶数,则 glogg(a)/2 即为所求。

特殊情况的解法:由 Euler 准则,a(p1)/21(modp)。特别地,如果 p12 是奇数,那么 a(p+1)/2a(modp)p+12 是偶数。令 x=a(p+1)/4 即可。

模意义下扩域

ω 是二次非剩余,i2ω(modp)。这样的 i=ω 在模 p 下的整数域中不存在。就像 1 在实数域上不存在使得数学家们定义复数域 {a+bia,b\R} 一样,我们定义数域 Zp[i]={a+bia,b[0,p1]},可以理解为模 p 下以 w 为虚数单位的复数。

该数域满足很多常见的代数性质,如乘法结合律,乘法分配律等,在此不一一证明。其上的运算可以看成模 p 下的多项式的运算,只不过 x2=ω

这个技巧常见于模意义下的二次非剩余强制开方。例如 5 在模 109+7 下没有二次剩余,所以,在计算斐波那契数列模 109+7 时,需要将一个数表示为 a+b5

由 Euler 准则,ip1ω(p1)/21(modp),所以 ipi(modp)

Cipolla

我们希望找到一个数的偶数次幂等于 a。考虑 Lucas 定理的引理,

(b+x)p+1(bp+xp)(b+x)(modp)

b 是模 p 下的整数,x 是二次非剩余 ω 开根,记为 i,则

(b+x)p+1(bi)(b+i)b2ω(modp)

也就是说,如果能找到正整数 b 使得 ω=b2a 是二次非剩余,那么 (b+i)(p+1)/2 就是 aZp[i] 上的二次剩余。

  • (x+y)pxp+yp(modp)Zp[i] 上的正确性:使用二项式定理展开后,考虑到对任意 a+biZp[i]p(a+bi)0+0i0(modp)
  • 假如存在 (c+di)2a(modp),对比两侧虚部得到 2cd0(modp)。如果 d0,那么 c=0,即 d2ωa(modp)。因为 a 是二次剩余,所以 ω 是二次剩余(Legendre 符号的完全积性),矛盾。因此 d=0,即 aZp[i] 上的二次剩余也一定是在 Zp 上的二次剩余(p 的二次剩余),也即 (b+i)(p+1)/2 的虚部为 0
  • 考虑 b2a 是二次剩余的 b 的数量,即存在 x[1,p1] 使得 b2ax2(modp)(bx)(b+x)a(modp)。固定 b+x[1,p1] 得到 bx,可以唯一求出一组解 (b,x),而固定 b 之后对应的 x 恰有两个(开根有两个解)。去掉 x=0 的两种情况(b2a(modp)),可知 b 的数量为 (p1)22。因此,使得 b2a 是二次非剩余的 b 的数量为 pp322=p12。这说明期望随机两次 b 即可得到 ω

综上,Cipolla 算法的时间复杂度为 O(logp)。模板题 P5491 代码。

#include <bits/stdc++.h>
using namespace std;
mt19937 rnd(1064);
int T, n, p, b, w;
struct comp {
  int a, b;
  comp operator + (comp c) {
    return {(a + c.a) % p, (b + c.b) % p};
  }
  comp operator * (comp c) {
    return {(1ll * a * c.a + 1ll * b * c.b % p * w) % p, (1ll * a * c.b + 1ll * b * c.a) % p};
  }
};
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % p;
    a = 1ll * a * a % p, b >>= 1;
  }
  return s;
}
comp ksm(comp a, int b) {
  comp s = {1, 0};
  while(b) {
    if(b & 1) s = s * a;
    a = a * a, b >>= 1;
  }
  return s;
}
void solve() {
  cin >> n >> p;
  if(!n) return puts("0"), void();
  if(ksm(n, p - 1 >> 1) == p - 1) return puts("Hola!"), void();
  do b = rnd() % p; while(ksm(w = (1ll * b * b - n + p) % p, p - 1 >> 1) != p - 1);
  int x = ksm({b, 1}, p + 1 >> 1).a;
  cout << min(x, p - x) << " " << max(x, p - x) << "\n";
}
int main() {
  cin >> T;
  while(T--) solve();
  return 0;
}

除了 Cipolla 以外,还有能够简单扩展到模质数幂的 Tonelli-Shanks 算法,此处不做介绍,详见参考资料。

7.4.2 任意模数下开平方根

奇质数幂

设模数为 pka=pca。特别地,由 7.2 小节,a 是二次剩余当且仅当 c 是偶数且 a 是二次剩余,可以通过 ax2(modpkc) 的解乘以 pc/2 得到 ax2(modpk) 的解,所以接下来只考虑 apkap 的情况。

处理模质数幂的第一步一般是先把答案在模质数下求出来。设 xa 在模 p 下的平方根,则 x2a0(modp)。设 aa 在模 pk 下的平方根。将同余式两侧 k 次方,得

(x2a)k(x+a)k(xa)k0(modpk)

因为 ap,所以 a 不是 p 的倍数。于是 x+axa 不可能都是 p 的倍数,那么一定恰有一个是 p 的倍数。

a 看成未知数待定,展开 (x+a)k 得到 u+va,那么展开 (xa)k 的结果一定是 uva。于是

(u+va)(uva)u2v2a0(modpk)

又因为 u+vauva 恰有一个是 p 的倍数,所以 va 一定不是 p 的倍数,即 vp,于是 vpk。根据上式,u2v2a(modpk),令 yuv1(modpk),则 y2u2v2a(modpk)

求出一个平方根后,取相反数得到另一个,即得所有平方根。

时间复杂度 O(logp)

2 的幂

类似地,我们只考虑 a 是奇数的情况。

如果 ax2(mod2k+1),那么 ax2(mod2k)。因此,如果 xa 在模 2k+1 下的平方根,那么它一定是 a 在模 2k 下的平方根。

因此,考虑 a 在模 2k 下的所有平方根 x1,,xt,只需检查 xixi+2k 是否为 a 在模 2k+1 下的平方根即可。由 7.2 小节,t4

显然地,该算法可以求出所有平方根。如果认为平方根的数量是常数,则时间复杂度 O(logp)

也可以使用类似开 N 次方根的算法,时间复杂度较劣。见 7.5 小节。

任意模数

使用 CRT 合并即可。如果要求出所有平方根,当 a 是二次剩余时 x2a(modpk) 的所有解的求法已经在上文介绍了,还需考虑当 x20(modpk) 的所有解。这个问题的答案是显然的:x=pkx,其中 2kk

*7.5 任意模数开 N 次根

在掌握了足够多二次剩余相关性质和开平方根的算法之后,我们考虑解决更一般 N 次剩余问题:xna(modp)。其中,a=0 的情况可类似开平方根推导,不再赘述。

奇质数

gp 的原根,ap。使用 BSGS 求出 d=logga。如果 xgkn 次根,那么 nkd(modp1),解这个同余方程即可。

奇质数幂

设模数为 pka=pca。如果 x=pdxn 次根,那么 c=ndn 的倍数。除掉这部分之后解 (x)na(modpkc),类似模奇质数用原根 + BSGS 即可。

*2 的幂

类似地,容易简化至 a2k3 的情况。

开平方根时的递推算法在这里不适用了,因为递推过程中解的数量可能非常大,但最终它们都无法成为 n 次根。这导致就算保证了最终答案的数量,递推算法的复杂度依然无法保证,见 这个帖子

根据 6.3 小节的结论,所有 a2 可以写成唯一的 (1)c5d。我们使用 BSGS 求出这样的表示:如果 log5a 存在,那么 c=0d=log5a。否则 c=1d=log5(a)。如果 x=(1)c5dn 次根,那么有同余方程组

{ncc(mod2)ndd(mod2k2)

两个方程独立,分别求解后再合并即可。

任意模数

使用 CRT 合并即可。如果认为平方根的数量是常数,则时间复杂度 O(p)

7.6 总结

N 次剩余相关的知识比较难,为了更好地理解这些结论和算法,当模数有原根时,抓住原根等价于单位根的性质,想象一个环是很有帮助的。每个数编号为其以 g 为底的离散对数,两个数相乘相当于环上的编号相加。一个数的编号给出了它的很多性质,编号的奇偶性决定了其能否表示成一个数的平方,也引出了 Euler 准则。这对理解高次剩余同样有帮助。是时候再次贴出这张图了。

qlTzex.png

接下来我们介绍了模意义下开根,这个问题有三个维度:

  • a=0 / a 是剩余。
  • 模质数 / 奇质数幂 / 2 的幂 / 任意合数。
  • n=2 / n2

对于不同的情况有复杂度和复杂程度不同的解法。n=2 时可以做到 O(logp),而 n2 时本文仅介绍了 O(p)

在模数没有原根的时候,我们以 6.3 小节 (Z/nZ)× 的结构为基础,同样能够把开方操作转化为指数上简单的同余方程。一个例子是模 2 的幂下开 N 次根。

限于篇幅和 OI 方面的应用,还有很多内容没有介绍,如二次互反律。感兴趣的读者可以查阅资料自行学习。

7.7 例题

*P6610 [Code+#7] 同余方程

根据 μ(p)0 不难想到使用 CRT 将问题拆成模奇质数。每个奇质数的解的个数的积即为所求。设 p 是奇质数。

问题转化为求 x 能被拆成多少组 a,b 的和,使得 a,b 均为二次剩余或 0。其中,每个二次剩余的贡献系数是 2(一个二次剩余可表示为两个不同的数的平方),0 的贡献是 1,每组 (a,b) 的贡献即 ab 分别的贡献系数之积。

我们需要数学公式而非判断式来描述答案,因为判断式不可化简。

二次剩余对应 20 对应 1,二次非剩余对应 0,这使得我们想到 Legendre 符号:注意到 (ap)+1 等价于使得 x2a(modp)x 的个数,因此答案为

a+bx((ap)+1)((bp)+1)

1p1的二次剩余和二次非剩余个数相等,所以 i=0p1(ip)=0。再根据 bxa(modp) 和 Legendre 符号的完全积性,上式化简为

p+a=0p1(axa2p)

接下来是一步巧妙转化。考虑到给 axa2 除以 a2 并不影响二次剩余性,因此答案即

p+a=1p1(xa1p)

显然,对于 x0a[1,p1]xa 取遍 [1,p1],即 xa1 取遍 [0,p2]。故 x0 时,答案等于

p(1p)

根据 Euler 准则,上式即 p(1)p12

对于 x=0,答案 p+a=1p1(1p)p+(p1)(1)p12

公式足够简洁,可以 O(1) 计算。时间复杂度 O(p+nω(p))代码

*P8457 「SWTR-8」幂塔方程

题解

CHANGE LOG

  • 2021.12.6:重构文章,删去线性筛部分,修改部分表述。
  • 2022.3.15:重构文章。
  • 2022.3.24:新增 Wilson 定理,素数在阶乘和组合数中的幂次,阶与原根,二次剩余和 Lucas 定理。
  • 2022.6.12:勘误,修改表述。
  • 2022.6.22:将数论分块移出本文。
  • 2024.2.22:增加在线 O(1) 逆元。
  • 2025.1.22:重构文章,移除欧拉函数。
  • 2025.1.27:补充任意模数开平方根,开 N 次方根和 (Z/nZ)× 的结构性质。

参考资料

第一章

第二章

第六章

第七章

全文

posted @   qAlex_Weiq  阅读(10129)  评论(11编辑  收藏  举报
相关博文:
阅读排行:
· 阿里最新开源QwQ-32B,效果媲美deepseek-r1满血版,部署成本又又又降低了!
· 单线程的Redis速度为什么快?
· SQL Server 2025 AI相关能力初探
· AI编程工具终极对决:字节Trae VS Cursor,谁才是开发者新宠?
· 展开说说关于C#中ORM框架的用法!
点击右上角即可分享
微信分享提示