初等数学瞎扯Ⅰ：同余相关

0. 基本定义

• 整除：若对于 \(a,b\in \mathbb{Z}\) 且 \(b \neq 0\)，有 \(a\bmod b=0\)，则我们称 \(b\) 整除 \(a\) 或 \(a\) 被 \(b\) 整除，记作 \(b\mid a\)，否则我们称 \(b\) 不整除 \(a\) 或 \(a\) 不被 \(b\) 整除，记作 \(b \nmid a\)。

• 同余：若对于 \(a,b,p\in \mathbb{Z}\) 且 \(p \neq 0\)，有 \(a\bmod p=b\bmod p\)，则称 \(a\) 与 \(b\) 在模 \(p\) 意义下同余，记作 \(a\equiv b\pmod p\)，否则称 \(a\) 与 \(b\) 在模 \(p\) 意义下不同余，记作 \(a\not\equiv b\pmod p\)。本文中不会特别区分 \(=\) 和 \(\equiv\)，故可能存在混用情况，望见谅。

• 质数集：质数的集合，记作 \(\mathbb{P}\)

• 最大公约数：对于两个数 \(a,b\)，我们称数字 \(d\) 是 \(a,b\) 的最大公约数当 \(d\) 是整除 \(a,b\) 的最大整数，记作 \(\gcd(a,b)\)，简记为 \((a,b)\)。

• 互质：两个数 \(a,b\) 互质当且仅当 \((a,b)=1\)，记作 \(a\bot b\)。

• 质因子次数符号：把 \(n\) 分解质因数后数字 \(p\) 的幂次记作 \(\nu_p(n)\)

• 各位数之和符号：在 \(p\) 进制下的数 \(n\) 各位数字的和记作 \(s_p(n)\)，例如，\(s(123)=6\)。

• 完系与缩系：对于一个数 \(m\)，我们称一个集合 \(S\) 是完系，当且仅当 \(\forall n\in\mathbb{Z},T=\{x|x\in S,n\equiv x\pmod m\},|T|=1\)。我们称一个集合 \(S\) 是缩系，当且仅当 \(\forall n\in\mathbb{Z},n\bot m,T=\{x|x\in S,n\equiv x\pmod m\},|T|=1\)。

1. 常见定理

1-1. 费马小定理

先给出两个引理。

引理 1：当 \(p\) 是质数时，其因子只有 \(1\) 和 \(p\) 两个。因此，若两个数相乘是 \(p\) 的倍数，其中必然至少有一个是 \(p\) 的倍数。

引理 2：对于 \(a\in \mathbb{Z},p\in \mathbb{P},p\nmid a\)，则 \(\forall i,j\in [1,p),i\neq j,ai\not\equiv aj\pmod p\)。

考虑反证，若 \(ai\equiv aj\pmod p\)，则 \(a(i-j)\equiv 0\pmod p\)，与引理 1 矛盾。

引理二的另一表述为，对于 \(a\in \mathbb{Z},p\in \mathbb{P},p\nmid a\)，\(\forall i\in [1,p),ai\) 的取值互不相同且仍为 \([1,p)\).

所以，对于 \(a\in \mathbb{Z},p\in \mathbb{P}\)，有 \(\prod ai\equiv\prod i\pmod p\)，所以 \(a^{p-1}\equiv1\pmod p\)。这就是费马小定理。

1-2. 欧拉定理

考虑将费马小定理推广到非素数的形式。首先我们试图将引理 1 推广到更一般的形式。

引理 1 拓展：若 \(a,b,p\in \mathbb{Z}\) 且 \(p \mid ab\)，则 \(a,b\) 至少有一个与 \(p\) 不互质。

我们设在 \([1,p]\) 中有 \(x\) 个数与 \(p\) 互质，且构成的数字集合为 \(A=\{a_1,a_2,...,a_x\}\)。此时引理 2 可以被推广为

引理 2 拓展：对于 \(a,p\in \mathbb{Z},p\nmid a\)，则 \(\forall i,j\in A,i\neq j,ai\not\equiv aj\pmod p\)。

证明同上。对应的另一表述为，对于 \(a,p\in \mathbb{Z},p\nmid a\)，\(\forall i\in A,ai\) 的取值互不相同且仍属于 \(A\)。

那么此时仿照上文，有 \(\prod aa_i\equiv\prod a_i\pmod p\)，所以 \(a^{|A|}\equiv1\pmod p\)。其中 \(|A|\) 是 \(A\) 集合的大小，我们用 \(\varphi(x)\) 表示对于数字 \(x\) 对应 \(A\) 集合的大小。\(\varphi(x)\) 为欧拉函数。代入验证可知当 \(p\) 为质数时，此定理依然成立。

关于欧拉函数，详见 初等数学瞎扯Ⅲ：数论函数与筛法，这里不再赘述。

1-3. 威尔逊定理

下文中，若无特殊说明，默认 \(p\in \mathbb{P}\)。

1-3-1. 基本威尔逊定理

\(\prod_{i\in[1,p)} i\equiv-1\pmod p\Leftrightarrow p\in \mathbb{P}\)。

说人话，\((p-1)!\equiv -1\pmod p\)。

为证明上式，先证明一个引理。

引理 1：若 \(x\) 在模 \(p\) 意义下逆元为 \(y\)，则 \(y\) 在模 \(p\) 意义下逆元为 \(x\)。

证明：由已知，\(xy=1\)，则 \(x^{-1}\equiv y,y^{-1}\equiv x\)。

另外可以注意到，若 \(x\neq 1\)，则 \(x\neq y\)。

你可能会说，这么看来，\((p-1)!\equiv1\pmod p\)，但是对于 \(i=1\) 其逆元为 \(1\)，也就是对于 \(p\neq 2\)，逆元两两配对后会剩下一个元素，考虑这个元素是谁，注意到逆元也是属于 \([1,p)\) 的，也就是说，设元素为 \(x\)，有 \(x^2\equiv 1\)，解得 \(x=\pm 1=\{1,p-1\}\)。所以 \((p-1)!\equiv p-1\equiv -1\pmod p\) 在 \(p\neq 2\) 下成立。特别检验 \(p=2\) 依然成立。

但这还没完，上面只说明了必要性，下面来证明充分性。考虑大力分讨。

• 若 \(x\) 为完全平方数且 \(x>4\)，设 \(x=a^2\)，则 \(a\) 与 \(2a\) 存在与 \((x-1)!\) 中，故 \((x-1)!\equiv 0\pmod x\)。
• 若 \(x\) 为合数且不是完全平方数，设 \(x=pq(p,q>1)\)，则 \(p\) 与 \(q\) 存在与 \((x-1)!\) 中，故 \((x-1)!\equiv 0\pmod x\)。
• 若 \(x=4\)， \(3!\equiv 2\pmod 4\)。

得证。

所以你如果想求一个数是不是 4，你可以大力计算阶乘然后判断取模后是不是 2。

1-3-2. 威尔逊定理扩展

我们考虑一种修正阶乘，记作 \((n!)_p\)，这表示在 \(n\) 的阶乘中去掉所有含 \(p\) 的因子后的值。例如，\((7!)_2=1\times1\times3\times1\times5\times3\times5\bmod2\)，我们考虑该式子在模 \(p\) 意义下的值。我们先考虑素数幂次的取值，即 \((p^k!)_p\)。

首先试图把原式展开

\[\begin{aligned} &(p^k!)\\ =&1\times 2\times...\times(p-1)\times1\\ \times&(p+1)\times (p+2)\times...\times(2p-1)\times2\\ \times& ...\times1\\ =&1\times 2\times...\times(p-1)\times1\\ \times&1\times 2\times...\times (p-1)\times2\\ \times&1\times 2\times...\times (p-1)\times1 \end{aligned} \]

然后用威尔逊定理的普通形式代入，得到原式变为

\[(-1)^{p^{k-1}}\times1\times2\times... \]

注意到剩下的东西也是修正阶乘的形式，即为 \((p^{k-1}!)_p\)，答案可以递归下去，得到如下结果

\[(-1)^{p^{\dfrac{k(k-1)}{2}}} \]

那么对于任意正整数 \(n\)，只需要计算遗留下来的其他东西即可，注意到个数不会超过 \(p\)，所以时间复杂度是 \(O(p)\) 的，尽管这看起来没啥大用。

1-3-3. 威尔逊定理的其它拓展

求 \((p^k!)_p\)，在这里模数是 \(p^k\)。

首先考虑 \(p>2\) 的情况。仍去解方程 \(x^2\equiv 1\pmod {p^k}\)，得到 \(x=\pm 1\)，结果仍为 \(-1\)。

考虑 \(p=2\)，移项得 \(x^2-1=(x+1)(x-1)\equiv 0\pmod {2^k}\)，\(x=\pm 1\) 仍成立，此外有解 \(x=2^{k-1}\pm1\)，故答案为 \(1\)，但应当注意到 \(k<3\) 时，有重根存在，故答案为 \(-1\)。

综上，可得出如下结论

\[(p^k!)_p \equiv \begin{cases} 1 & p = 2 \land k \geq 3 \\ -1 & {\rm otherwise}\end{cases} \]

1-4. 素数在阶乘中的幂次

求 \(\nu_p(n!)\)

首先考虑朴素算法。每次从 \(n!\) 中提取出所有模 \(p\) 为 \(0\) 的数，然后将其除 \(p\)。重复此流程，知道答案不再改变。第 \(i\) 轮中提取出 \(\lfloor \dfrac{n}{p^i} \rfloor\) 个数，故答案为

\[\sum_{i=1}^{p^i\leq n}\lfloor \dfrac{n}{p^i} \rfloor \]

时间复杂度 \(O(\log n)\)。

考虑另外一种表达方式，我们考虑在 \(p\) 进制意义下，每一位的贡献，令 \(l\) 为此时 \(n\) 的长度（姑且这么称呼），令 \(n=\sum_{i=0}^p a_ip^i\)，此时有

\[\begin{aligned} ans&=\sum_{i=1}^c\sum_{j=0}^c\lfloor \dfrac{a_jp^j}{p^i} \rfloor\\ &=\sum_{i=1}^c\sum_{j=i}^ca_jp^{j-i}\\ &=\sum_{j=0}^ca_j\sum_{i=0}^{j-1}p^j\\ &=\dfrac{\sum_{j=0}^ca_jp^j-a_j}{p-1}\\ &=\dfrac{n-s_p(n)}{p-1} \end{aligned} \]

注意到 \(p=2\) 时，原式等于 \(n-\text{popcount}(n)\)，这使得求解复杂度可以降为常数级。

1-5. 库默尔定理

由 1-4 提到的幂次可知，对于组合数 \(\dbinom{n}{m}\)，\(\nu_p(\dbinom{n}{m})\) 可以用如下方式计算。

\[\begin{aligned} \nu_p(\dbinom{n}{m})&=\nu_p(\dfrac{n!}{m!(n-m)!})\\ &=\nu_p(n!)-\nu_p(m!)-\nu_p(n-m)!\\ &=\dfrac{s_p(m)+s_p(m-n)-s_p(n)}{p-1} \end{aligned} \]

考虑 \(n-m\) 和 \(m\) 在 \(p\) 进制下相加得到 \(n\) 的过程，其中 \(p\) 是质数。每产生一次进位，各位数字之和就会减少 \(p-1\)，因此 \(p\) 在组合数 \(\dbinom{n}{m}\) 中的幂次 \(\nu_p(\dbinom{n}{m})\) 即 \(p\) 进制下 \(n-m\) 加上 \(m\) 需要进位的次数，也等于 \(n\) 减去 \(m\) 需要借位的次数。

这个定理将会在 Lucas 定理中用到。

1-6. 步数-子环长定理

对于长度为 \(n\) 的环，从一个点走每次走 \(k\) 步能走到的点数恰为 \(\dfrac{n}{\gcd(n,k)}\)，且形成的本质不同的环有 \(\gcd(n,k)\) 个。

直觉理解是很对的，但是我们还是严谨证明一下。

不妨考虑 \(n\perp k\) 的情况，若走 \(n-1\) 步经过的节点均不同，即可说明可以走到所有节点。观察可知原命题等价于 \(r+ik\not\equiv r+jk\pmod n\) 两边消掉 \(r\) 后即可转化到费马小定理的引理 2，这里不再赘述。

考虑一般情况，我们对模 \(k\) 的同余类进行重编号，假设当前考察的同余类为 \(r\)，则共有 \(\dfrac{n}{\gcd(n,k)}\) 个数。我们不妨将其编号为 \(1,2,...,\dfrac{n}{\gcd(n,k)}\)，令 \(n'=\dfrac{n}{\gcd(n,k)}\)，则此时在原环上走 \(k\) 步等价于在子环上走 \(\dfrac{k}{\gcd(n,k)}\) 步。令 \(k'=\dfrac{k}{\gcd(n,k)}\)，则有 \(k'\perp n'\)，此时由上文可知定理成立。

2. 取模意义下的基本运算

2-1. 加法，减法与乘法

在普通四则运算下的运算方式仍然成立，不过多提及。

2-2. 除法

注意到 \(\dfrac{a}{b}\not\equiv \dfrac{a\bmod p}{b\bmod p}\pmod p\)，所以我们需要考虑将其变化为乘法运算。

\(\dfrac{a}{b}=a\dfrac{1}{b}=ab^{-1}\)，我们只需要计算 \(b^{-1}\) 即可。具体计算见 2-3 乘法逆元。

2-3. 乘法逆元

2-3-1. 欧拉定理形式

由 \(a^{\varphi(p)}\equiv1\pmod p\)，得 \(a^{-1}=a^{\varphi(p)-1}\)，可用快速幂解决，时间复杂度 \(O(\log p)\)。

关于快速幂，详见 初等数学瞎扯Ⅱ：辅助工具，这里不再赘述。

2-3-2. 线性求逆元

假设 \([1,i)\) 的逆元已知。设 \(p=ki+r(0\leq r<i)\) \(ki+r\equiv 0\pmod p\)，两边同时除以 \(ir\)，得 \(i^{-1}\equiv -\dfrac{k}{r}\equiv \lfloor \dfrac{p}{i}\rfloor r^{-1}\pmod p\)。时间复杂度 \(O(n)\)。

2-3-3. 线性任意 \(n\) 个数逆元

令 \(a_1,a_2,...,a_n\) 为所求数，对 \(a_1,a_2,...,a_n\) 做前缀积，得到 \(p_1,p_2,...,p_n\)，对 \(p_n\) 做 1-3-1 中提及的求法，设得到 \(inv_n\)，即可得到 \(a_n\) 的逆元为 \(inv_n\times p_{n-1}\)，以及 \(inv_{n-1}=inv_n\times a_n\)，递归求解即可，时间复杂度 \(O(n+\log p)\)。

2-4. 开根

开根与除法同理，也不能直接取模，为此我们也要对开根进行转化。注意到 \((\sqrt a)^2=a\)，也就是说，若 \(x^2\equiv a\pmod p\)，则 \(x\equiv \sqrt a\pmod p\)，其中，我们称 \(x\) 为模 \(p\) 意义下的二次剩余。

具体计算见二次剩余。

3. 二元一次不定方程组问题

3-0. 问题简述

解二元一次不定方程。

3-1. 前置知识—裴蜀定理

对于一个二元一次方程，我们首先要判定解的存在性。

一个直观的结论是，对于方程 \(ax+by=c\)，若 \(c\nmid (a,b)\)，则该方程一定无整数解，原因是 \(\forall a,b\in \mathbb Z,(a,b)\mid ax+by\)，显然此时 \(ax+by\neq c\)。

我们试图对这个结论做进一步推广，即若 \(c\mid (a,b)\)，则原方程一定有解。不妨令 \(d=(a,b),a'=\dfrac{a}{d},b'=\dfrac{b}{d},c'=\dfrac{c}{d}\)，此时有 \(a'\bot b'\)，则我们若能证明 \(a'x+b'y\) 可以取遍任何数即可。

不难发现上文等价于存在一组 \(x\) 满足 \(a'x_i\bmod b'\) 取遍 \([0,b'-1]\)（这意味着对于每个数 \(n\)，存在一个解为 \(x_0\in x,y_0=\dfrac{n-a'x_0}{b'}\)，其中 \(x_0\) 满足 \(a'x_0\equiv n\pmod {b'}\)），此时注意到 \(a\bot b\) 是一个很强的性质，我们令 \(x=\{0,1,...,b'-1\}\)，可以用和费马小定理/欧拉定理一样的方式反证集合 \(x\) 中的数模 \(b'\) 互不相同，这里不再赘述。

综上，我们得到了二元一次不定方程解的存在性的充要条件，\(c\mid (a,b)\)。

3-2. 扩展欧几里得算法

我们再次审视一下上面的方程，不难发现我们可以视为解方程 \(ax+by=(a,b)\)，然后将结果扩大 \(\dfrac{c}{(a,b)}\) 倍。

扩展欧几里得算法是递归构造的算法，我们不妨回忆一下我们是如何辗转相除求 \(\gcd\) 的，我们每次将 \((a,b)\) 改为 \((b,a\bmod b)\)，然后递归求解，注意到此时二者的 \(\gcd\) 不变，这意味着我们如果得知了方程 \(bx+(a\bmod b)y=(a,b)\) 的解 \((x_0,y_0)\)，一定有一组满足 \(ax+by=(a,b)\) 的解为 \(x_0'=y_0,y'_0=x_0-\lfloor\dfrac{a}{b}\rfloor y_0\)，具体证明代入即可。这里用到了一步比较常见和重要的变换是 \(a\bmod b=a-b\lfloor\dfrac{a}{b}\rfloor\)。

3-3. 代码

先是容易理解的。

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

然后是短一点的

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

3-4. 通解的形式与特解范围

通解可表示为如下的形式

\[\begin{cases} x=x_0+\dfrac{b}{d} k\\ y=y_0-\dfrac{a}{d} k \end{cases} \]

其中 \(x_0,y_0\) 是一组通解。

关于特解的大小具体证明可以参见 ycx 的博客，有结论为其通解与 \(a,b\) 最大同阶。

3-5. 关于扩欧在逆元上的应用

• 对于求解逆元来说，费马定理只在 \(p\) 为质数的时候成立。而对于一类模数不为质数，且在 \(10^8\) 范围内的模数来说，扩展欧拉定理需要处理 \(\varphi\)，而常见的处理 \(\varphi\) 的复杂度是 \(O(\sqrt n)\)（处理所有约数），当然我们也许可以使用 pollard_pho 算法做到 \(O(n^{\dfrac{1}{4}})\) 分解质因数然后得到每个质因子幂次后根据约数个数和求 \(\varphi\)，但这样让我们看起来很呆。考虑问题的实质是什么，求解在模 \(p\) 意义下 \(d\) 的逆元，即需要找到一个数 \(x\)，使得 \(xd\equiv 1\pmod p\)，而同余方程可以转二元一次不定方程求解，即 \(dx+py=1\)，直接扩欧上，求出一组通解即可。

3-6. 一些奇怪的结论

你说的对，但是只有通解是没有用的。

一些题目中我们要求解需要满足一些性质，比如最小正整数解啥的，这里直接给出推导和结论。

• \(x\) 的最小正整数解 \(x_{min}\)

我们试图对通解进行调整，具体而言，在
\( \begin{cases} x=x_0+\dfrac{b}{d} k\\ y=y_0-\dfrac{a}{d} k \end{cases} \) 这个式子的基础上，找到一个最小的 \(k\)，满足 \(x=x_0+\dfrac{b}{d}k\ge 1\)。

求解得到，\(k>\lceil\dfrac{d(1-x_0)}{b}\rceil\)，因此 \(k\) 的最小值为\(\lceil\dfrac{d(1-x_0)}{b}\rceil+1\)

所以 \(x\) 的最小正整数解为
\( x=x_0+\dfrac{b}{d}(\lceil\dfrac{d(1-x_0)}{b}\rceil+1) \)

• \(y\) 的最小正整数解 \(y_{min}\)

与 \(x\) 的推导同理。毛毛估一下是 \(y=y_0+\dfrac{a}{d}(\lceil\dfrac{d(1-y_0)}{a}\rceil+1)\)。

• \(x\) 的最大正整数解 \(x_{max}\)

即 \(y\) 取到最小时的 \(x\)。

• \(y\) 的最大正整数解 \(y_{max}\)

即 \(x\) 取到最小时的 \(y\)。

• 正整数解的个数

我们先将解调整为 \(x_{min},y_{max}\)，则答案为把 \(y\) 从 \(y_{max}\) 调整到 \(y_{min}\) 的次数，大概是 \(\dfrac{d(y_{max}-y_{min})}{a}\)。

4. 离散对数问题

4-0. 问题简述

求解 \(a^x\equiv b\pmod p\) 中 \(x\) 的最小取值。

或者说，我们要找到一个 \(x\)，满足 \(x\equiv \log_ab\pmod p\)

4-1. 大步小步算法

大步小步算法又叫做 BSGS(Baby Step Giant Step)，他的普通形式适用于 \(a\bot p\) 的情况。

BSGS 算法运用了根号平衡的思想，我们令 \(x=kT-c\)，且 \(c\in[0,T-1]\)，则原式即为 \(a^{kT-c}\equiv b\pmod p\)，移项后为 \(a^{kT}=a^cb\pmod p\)，直接枚举 \(a^{kT}\) 的取值，由欧拉定理可知最多有 \(\dfrac{p}{T}\) 种，用哈希表存下这些值，然后枚举 \(c\)，判断其是否在哈希表内出现。而 \(c\) 的取值有 \(T\) 种，总复杂度为 \(O(\dfrac{p}{T}+T)\)，取 \(T=\sqrt p\) 得到最优复杂度 \(O(\sqrt p)\)。

unordered_map<int,int>Map;
int BSGS(int a,int b,int p){
	Map.clear();a%=p,b%=p;
	int t=sqrt(p)+1,v=b,s=1;
	for(int i=1;i<=t;i++)Map[v=1ll*v*a%p]=i,s=1ll*s*a%p;
	for(int i=1,v=s;i<=t;i++,v=1ll*v*s%p)
		if(Map.count(v))return i*t-Map[v];
	return -1;
}

4-2. 扩展大步小步算法

现在如果没有 \(a\bot p\) 的限制，考虑如何去做。

考虑我们在 gcd 里干的事情，把不互质的数变成互质的。

令 \(d=(a,p)\)，那么原方程变为 \(\dfrac{a}{d}a^{x-1}\equiv \dfrac{b}{d}\pmod{\dfrac{p}{d}}\)，而此时 \((\dfrac{a}{d},\dfrac{p}{d})=1\)，故有逆元存在，将两边同时除以逆元，得到 \(a^{x-1}=\dfrac{b}{d}(\dfrac{a}{d})^{-1}\pmod {\dfrac{p}{d}}\)，不妨令 \(b'=\dfrac{b}{d}(\dfrac{a}{d})^{-1},p'=\dfrac{p}{d}\)，则原式变为 \(a^{x-1}\equiv b'\pmod{p'}\)，此时又变成了原来的形式，直到 \(p\) 和 \(a\) 互质。但 \(p'\) 至少是原来的一半，故这个操作将会递归不超过 \(\log p\) 轮，最后用一遍 BSGS 求解即可，复杂度是求解逆元的 \(O(\log^2 p)\) 和一遍 BSGS 的 \(O(\sqrt p)\)，由于哈希表的大常数和实际增长速度我一般喜欢直接当成根号算法。

在做上述的过程中需要特判，比如 \(b=1\) 和 \(b\nmid d\) 的情况。另外就是假设上述过程进行了 \(k\) 轮，则答案需要 \(+k\)。

int exgcd(int a,int b,int &x,int &y){
	if(!b){x=1,y=0;return a;}
	int d=exgcd(b,a%b,y,x);y-=a/b*x;return d;
}
int getinv(int x,int p){
	int y;exgcd(x,p,x,y);
	return (x%p+p)%p;
}
unordered_map<int,int>Map;
int BSGS(int a,int b,int p){
	Map.clear();a%=p,b%=p;
	int t=sqrt(p)+1,v=b,s=1;
	for(int i=1;i<=t;i++)Map[v=1ll*v*a%p]=i,s=1ll*s*a%p;
	for(int i=1,v=s;i<=t;i++,v=1ll*v*s%p)
		if(Map.count(v))return i*t-Map[v];
	return -1;
}
int exBSGS(int a,int b,int p){
	int d=__gcd(a,p);a%=p,b%=p;int k=0;
	for(;d>1;d=__gcd(p,a%=p),k++){
		if(b==1)return k;if(b%d!=0)return -1;
		p/=d,b=1ll*b/d*getinv(a/d,p)%p;
	}int ans=BSGS(a,b,p);return ans>0?ans+k:-1; 
}

一个有趣的事实：我因为No Solution写成No solution 查了 20min BSGS。

5. 线性同余方程组

5-0. 问题简述

求解形如

\[\begin{cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \cdots \\ x \equiv a_k \pmod {m_k} \end{cases} \]

的线性同余方程组。

5-1. 中国剩余定理

中国剩余定理，简称 CRT(Chinese Remainder Theorem)。它用于求解模数两两互质的线性同余方程组，即对于任意 \(i\neq j\)，均有 \(m_i\bot m_j\)。

在CRT中，我们试图求出一个解集 \(x\)，对于 \(x\) 中的元素 \(x_i\)，均满足 \(x_i\) 是第 \(i\) 个方程的解，且对于其它模数 \(x_i\) 对其取模的值均为 \(0\)。不难发现,\(\sum\limits_{v\in x} v\) 即为原式的一个可行解。

我们试图对于单个方程求出一个解，令 \(M=\prod\limits_{i=1}^km_i\)，那么对于第 \(i\) 个方程而言，我们求出的答案必须满足下述性质。

• \(\dfrac{M}{m_i}\mid x_i\)
• \(x_i\equiv a_i\pmod {m_i}\)

首先试图满足第一个性质，直接令 \(x_i=p\dfrac{M}{m_i}\) 即可。

然后是第二个性质，我们不妨令 \(d_i\) 为 \(\dfrac{M}{m_i}\) 在模 \(m_i\) 意义下的逆元，则有 \(x_id_i\equiv p\pmod {m_i}\)，若令 \(p=a_i\) 则正好满足上述的第二条性质，得到答案 \(x_i=a_id_i\dfrac{M}{m_i}\) 对于每个方程跑一遍上述流程即可。复杂度 \(O(n\log n)\)。

int n,a[N],m[N],c[N],d[N];ll M=1,ans;
int exgcd(int a,int b,int &x,int &y){
	if(!b){x=1,y=0;return a;}
	int d=exgcd(b,a%b,y,x);y-=a/b*x;return d;
}
int getinv(int x,int p){
	int y;exgcd(x,p,x,y);
	return (x%p+p)%p;
}
ll CRT(){
	for(int i=1;i<=n;i++)M=M*m[i];
	for(int i=1;i<=n;i++){
		c[i]=M/m[i];
		ans+=a[i]*(M/m[i])%M*getinv(M/m[i]%m[i],m[i])%M;
	}return ans;
}

5-2. 扩展中国剩余定理

考虑对于模数不两两互质的情况，在这个算法里，我们每次试图合并两个同余方程。

考虑对于如下同余方程

\[\begin{cases} x\equiv a_1\pmod {m_1}\\ x\equiv a_2\pmod {m_2}\\ \end{cases} \]

第一个同余方程的通解可以表示为 \(x=a_1+km_1\)，将其代入第二个同余方程得到 \(a_1+km_1\equiv a_2\pmod{m_2}\)，即 \(k\equiv \dfrac{a_2-a_1}{m_1}\pmod{m_2}\)，\(k\) 直接 exgcd 求逆元即可。

合并后的方程即为 \(x\equiv a_1+km_1\pmod{\text{lcm}(m_1,m_2)}\)。

ll n,a[N],m[N],ans;
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1,y=0;return a;}
	int d=exgcd(b,a%b,y,x);y-=a/b*x;return d;
}
ll getinv(ll x,ll p){
	ll y;exgcd(x,p,x,y);
	return (x%p+p)%p;
}
ll exCRT(){
	ll A=a[1],M=m[1];
	for(ll i=2,d,c;i<=n;i++){
		d=__gcd(M,m[i]);c=(__int128)(a[i]-A%m[i]+m[i])%m[i];c/=d;
		A+=(__int128)c*getinv(M/d,m[i]/d)%(m[i]/d)*M;M*=m[i]/d,A%=M; 
	}return A;
}

6. 阶与原根

6-1. 阶

\(a\) 在模 \(m\) 意义下的阶定义为使得 \(a^x\equiv 1\pmod m\) 的最小正整数 \(x\)，记为 \(\delta_m(a)=x\)。

对于一个数 \(a\)，其在模 \(m\) 意义下存在阶的充要条件是 \(a\bot m\)，证明如下。

考虑一个数字 \(a\) 和 \(m\) 互质，若 \(a^x\equiv 1\pmod m\)，则 \(a^x\bot m\)（扩展欧几里得），而 \(a\)，\(m\) 不互质，故原命题不成立。

再考虑一个数字 \(a\bot m\)，则由扩展欧拉定理有一解为 \(\varphi(m)\)。

接下来给出阶的性质。设存在两个数字 \(a\bot m\)。

1. \(\delta_m(a)\mid\varphi(m)\)，证明如下。

若存在一个 \(\delta_m(a)\mid\varphi(m)\)，令 \(k=\varphi(m)\bmod \delta_m(a)\)，则 \(k\) 也是一个使得 \(a^x\equiv 1\pmod m\) 的值，且 \(k<\delta_m(a)\)，与上文矛盾。

2. 若 \(a^x\equiv 1\pmod m\) ，则一定有 \(\delta_m(a)\mid x\)，证明同上。

3. 若 \(a\) 在模 \(m\) 意义下的阶为 \(\delta_m(a)\)，则 \(\delta_m(a^k)=\dfrac {\delta_m(a)}{\gcd(\delta_m(a),k)}\)，证明如下。

按照原根的定义写出下式。

\[\begin{aligned} &(a^k)^{\delta_m(a^k)}\equiv a^{k\delta_m(a^k)}\equiv 1\pmod m\\ \Rightarrow &\delta_m(a)\mid k\delta_m(a^k)\\ \Rightarrow &\dfrac{\delta_m(a)}{\gcd(k,\delta_m(a))}\mid \delta_m(a^k)\\ \end{aligned} \]

再根据 \(a^{\delta_m(a)}\equiv 1\pmod m\)，可以得到下式。

\[\begin{aligned} &a^\dfrac{k\delta_m(a)}{\gcd(k,\delta_m(a))}\equiv (a^{\delta_m(a)})^\dfrac{k}{\gcd(k,\delta_m(a))}\equiv 1\pmod m\\ \Rightarrow &\delta_m(a^k)\mid \dfrac{\delta_m(a)}{\gcd(k,\delta_m(a))} \end{aligned} \]

综上，\(\delta_m(a^k)=\dfrac {\delta_m(a)}{\gcd(\delta_m(a),k)}\)。

4. \(a^1.a^2,...,a^{\delta_m(a)}\) 在模 \(m\) 意义下两两不同，证明如下。

若存在 \(i<j\) 满足 \(a^j\equiv a^i\pmod m\)，则 \(a^{j-i}\equiv 1\pmod m\)，而 \(j-i<\delta_m(a)\)，与阶的最小性矛盾，故不成立。

6-2. 原根

对于 \(m\in \mathbb{N}^*,g\in \mathbb{Z},(m,g)=1\)，我们称 \(g\) 是模 \(m\) 意义下的原根当且仅当 \(\delta_m(a)=\varphi(m)\)。

原根的若干次幂构成了一个模 \(m\) 意义下的缩系，这会在部分题目内有用。

原根判定定理：对于 \(m\ge 3,a\bot m\)，\(a\) 是 \(m\) 的原根的充要条件是 \(a^{\dfrac{\varphi(m)}{p}}\not\equiv 1\pmod m\)。证明如下。

• 必要性：若 \(a^{\dfrac{\varphi(m)}{p}}\equiv 1\pmod m\)，则 \(\delta_m(a)\leq\dfrac{\varphi(m)}{p}\)，与原根的定义矛盾。

• 充分性：不妨假设存在一个 \(a\) 满足 \(a^{\dfrac{\varphi(m)}{p}}\not\equiv 1\pmod m\)，则此时若 \(a\) 不是原根，则令 \(a^t\equiv 1\pmod m\)，此时有 \(t\mid \varphi(m)\)，不妨假设 \(\dfrac{\varphi(m)}{t}=x\)，令 \(x\) 的最小质因子为 \(p'\)，则有 \(a^\dfrac{\varphi(m)}{p'}\equiv 1\pmod m\)，与原假设不符，故舍去。

综上，我们证明了原命题。

原根存在定理：对于一个数 \(m\)，其存在原根的充要条件是 \(m=2,4,p^{\alpha},2p^{\alpha}\)，具体证明可见 oi-wiki - 原根 - 原根存在定理，感兴趣的读者可自行了解，这里由于证明复杂略去。

6-3. 阶与原根的推论

对于一个存在原根的数字 \(m\)，使得 \(\delta_m(a)=x\) 的数字 \(a\) 的个数 \(n(x)\) 为

\[n(x)= \begin{cases} 0&x\nmid \varphi(m)\\ \varphi(x)&x\mid \varphi(m) \end{cases} \]

第一部分在阶的性质中已经有提及，在这里略去，考虑性质二的证明，令 \(g\) 为 \(m\) 的一个原根，则任意一个满足 \(a\bot m\) 的 \(a\) 均可以表示为 \(g^k\)，此时 \(\delta_m(a)=\delta_m(g^k)=\dfrac{\delta_m(g)}{\gcd(\delta_m(g),k)}=\dfrac{\varphi(m)}{\gcd(\varphi(m),k)}=x\)，移项处理可得 \(\dfrac{\varphi(m)}{x}=\gcd(\varphi(m),k)\)，由欧拉函数的性质，可知答案为 \(\varphi(\dfrac{\varphi(m)}{\dfrac{\varphi(m)}{x}})=\varphi(x)\)。

同时可以得到原根个数定理：原根个数为 \(\varphi(\varphi(m))\)。最小原根的数量级约为 \(m^{0.25+\epsilon}\)，由王元院士证明，感兴趣的读者可以自行寻找资料。

因此我们得到一个求原根的算法，暴力枚举 \(\varphi\) 的所有质因子并试除，复杂度为 \(O(m+\log m+\sqrt[4]{m}\log m\omega(m))\)。

同时，不难发现令 \(n=\varphi(m)\)，则 \(g\) 与 \(m\) 的 \(n\) 次单位根 \(\omega\) 同构。因为 \(g^k\equiv 1\pmod m,\omega_1^k=1\)，且 \(g^i\times g^j=g^{i+j},\omega_i\times\omega_j=\omega_{i+j}\) 因此，在某种条件下可以用原根代替单位根做快速数论变换，即 NTT。

在这里给出求原根的代码。

const int M=1e6+7;
int minp[M],p;
void init(int lim){
	for(int i=2;i<=lim;i++)if(!minp[i])
		for(int j=i;j<=lim;j+=i)minp[j]=i;
}
int qp(int b,int m,int p){
	int res=1;
	for(;m;m>>=1,b=1ll*b*b%p)if(m&1)res=1ll*res*b%p;
	return res;
}
int getphi(int m){
	int res=m; 
	while(m!=1){
		int p=minp[m];while(m%p==0)m/=p;
		res=1ll*res/p*(p-1);
	}return res;
}
int getg(int m){
	int g=1;init(m);
	int phi=getphi(m),q=phi;
	while(1){
		g++;phi=q;
		while(phi!=1){
			int p=minp[phi];while(phi%p==0)phi/=p;
			if(qp(g,q/p,m)==1)goto there;
		}return g;there:;
	}
}

我们将会在后文的高次剩余（也许会有？）中进一步探讨原根和阶的使用，在这里告一段落。

7. 卢卡斯定理

7-0. 问题简述

求 \(\dbinom{n}{m}\bmod p\)，其中 \(n,m\) 较大，\(p\) 较小。

7-1. 普通卢卡斯定理

设两个数字 \(a,b\) 和质数 \(p\)，由 1-5 的库默尔定理可知 \(\dbinom{a}{b}\equiv 0\pmod p\) 当且仅当 \(a-b\) 在 \(p\) 进制意义下需要借位，等价于令 \(c=\lfloor\log_p a\rfloor,a=\sum\limits_{i=0}^cx_ip^i,b=\sum\limits_{i=0}^cy_ip^i\)，满足 \(\forall i\in[0,c],x_i>y_i\)。进一步的，由于当 \(i<j\) 时，\(\dbinom{i}{j}=0\)，我们考虑是否 \(\dbinom{a}{b}\equiv\prod\limits_{i=0}^c\dbinom{x_i}{y_i}\pmod p\)？

为证明上述问题，优先证明一个引理。

引理：\((x+y)^p\equiv x^p+y^p\pmod p\)，证明如下。

二项式定理拆开后，易证中间几项系数均为 \(p\) 的倍数，故得证。

同理可推出，\((1+x)^p\equiv 1+x^p\pmod p\)。

考虑 Lucas 定理的另一个表述，\(\dbinom{a}{b}\equiv\dbinom{\lfloor\dfrac{a}{p}\rfloor}{\lfloor \dfrac{b}{p}\rfloor}\dbinom{a\bmod p}{b\bmod p}\pmod p\)，注意到这是进制转换的形式，故两等式等价，我们只需要证明上式即可。

不妨用二项式定理将组合数变为多项式系数（有一点生成函数的思想），即 \((1+x)^a\) 的 \(b\) 次项系数，令 \(r=a\bmod p,d=\lfloor\dfrac{a}{p}\rfloor\)，则原组合数可写为 \([x^b](1+x)^{dp}(1+x)^r\equiv[x^b](1+x^p)^d(1+x)^r\pmod p\)。不难发现，此时的 \(b\) 次方项只能由 \(x^{\dfrac{b}{p}\times p}\cdot x^{b\bmod p}\) 得到，此时项前的系数为两项的乘积，即 \(\dbinom{d}{\lfloor \dfrac{b}{p}\rfloor}\dbinom{r}{b\bmod p}\)。将 \(r,d\) 回代，即可得到原等式。时间复杂度 \(O(\log n+p)\)。

const int N=1e5+7;
int T,n,m,mod,fac[N],ifac[N];
int qp(int b,int t=mod-2){
	int res=1;
	for(;t;t>>=1,b=(ll)b*b%mod)if(t&1)res=(ll)res*b%mod;
	return res;
}
int binom(int x,int y){return (x<0||y<0||x<y)?0:(ll)fac[x]*ifac[y]%mod*ifac[x-y]%mod;}
int lucas(int x,int y){
	if(x<mod)return binom(x,y);
	return (ll)lucas(x/mod,y/mod)*binom(x%mod,y%mod)%mod;
}
void solve(){
	read(n,m,mod);
	for(int i=fac[0]=1;i<mod;i++)fac[i]=(ll)fac[i-1]*i%mod;
	ifac[mod-1]=qp(fac[mod-1]);
	for(int i=mod-2;i>=0;i--)ifac[i]=(ll)ifac[i+1]*(i+1)%mod;
	printf("%d\n",lucas(n,m));
}
int main(){
	read(T);
	while(T--)solve();
	return 0;
}

7-2. 扩展卢卡斯定理

现在考虑模数 \(q\) 不是质数的情况。使用中国剩余定理，将 \(q\) 分解为若干质数的幂次，算出每一个的答案，然后拿中国剩余定理合出来在模 \(p\) 意义下的答案。

现在的问题转化为了求 \(\dbinom{n}{m}\) 模 \(p^k\) 意义下的结果。把式子展开，有 \(\dfrac{n!}{m!(n-m)!}\)。考虑首先在原式中去掉 \(p\)，变为 \(\dfrac{\dfrac{n!}{p^{\nu(n!)}}}{\dfrac{m!(n-m)!}{p^{\nu(m!)}p^{\nu((n-m)!)}}}p^{\nu(n!)-\nu(m!)-\nu(n-m)!}\)。后面的东西是易处理的，问题在前面那一堆。更进一步的，其可以被规约为求 $ \dfrac{n!}{p^{\nu(n!)}}$.考虑这个东西怎么算。

将所有数字分成 \(p\) 的倍数和非倍数讨论。令 \(d=\lfloor\dfrac{n}{p}\rfloor\)对于其倍数，我们每次除掉一个 \(p\)，变成了求 \(\dfrac{d!}{p^{\nu(d!)}}\)。这是一个子问题的形式，我们直接递归下去做即可，边界为 \(d=0\)。接下来考虑后者，引入威尔逊定理的修正阶乘的概念。注意到所有数在模 \(p^k\) 意义下构建了一个循环节，此时即为求 \((p^k!)_k^{\lfloor\dfrac{n}{p^k}\rfloor}(n\bmod p^k)!\) 模 \(p^k\) 的值。那么预处理出来 \([1,p^k]\) 所有阶乘即可。

最后用中国剩余定理合并即可。

const int N=2e5+7;
int fac[N];
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){int a,b;exgcd(x,p,a,b);return (a%p+p)%p;}
ll calv(ll x,int p){ll res=0;while(x)x/=p,res+=x;return res;}
int calc(ll x,int p,int pk){return !x?1:(ll)calc(x/p,p,pk)*(((pk&1)||(pk<=4))&&((x/pk)&1)?pk-1:1)%pk*fac[x%pk]%pk;}
int qp(int b,int m,int p){
	int res=1;
	for(;m;m>>=1,b=(ll)b*b%p)if(m&1)res=(ll)res*b%p;
	return res;
}
ll solve(int n,int m,int p,int pk){
	for(int i=fac[0]=1;i<pk;i++)fac[i]=(ll)fac[i-1]*(i%p?i:1)%pk;
	int cur=calv(n,p)-calv(m,p)-calv(n-m,p);
	int res=(ll)calv(n,p)*inv((ll)calv(m,p)*calv(n-m,p)%pk,pk)%pk;
	return (ll)res*qp(p,cur,pk)%pk;
}
int solve(int n,int m,int mod){
	int ans=0,all=1;
	for(int i=2;i*i<=mod;i++){
		if(mod%i)continue;
		int p=i,d=1;
		while(mod%p==0)mod/=p,d*=p;
		ans=ans+(ll)(solve(n,m,p,d)-ans%d+d)*inv(all,d)%d*all;
		all*=d;ans%=all;
	}
	if(mod!=1){
		int p=mod,d=1;
		while(mod%p==0)mod/=p,d*=p;
		ans=ans+(ll)(solve(n,m,p,d)-ans%d+d)*inv(all,d)%d*all;
		all*=d;ans%=all;
	}
	return ans;
}
ll n,m;int p;

int main(){
	read(n,m,p);
	printf("%d\n",solve(n,m,p));
	return 0;
}

end.本文引用&鸣谢

Alex_Wei's Blog

ycx's Blog

oi-wiki