数论学习笔记 Part 2
luuia
·
2024-03-09 21:28:08
·
算法·理论
最近一个月学习了一下数论的相关知识，写一篇笔记记录一下。
目录
Part 1
Part 2
Part 3
Part 4
Part 5
数论分块
数论分块可以在 O(\sqrt{n}) 的复杂度内快速计算含有除法下取整的式子，例如 \sum_{i=1}^n\lfloor\frac{n}{i}\rfloor 的值。
记 f(x) = \lfloor \dfrac{n}{x} \rfloor ，令 n = 10 ，观察 f 前几项的规律。
f(1) = 10,f(2) = 5,f(3) = 3,f(4) = f(5)=2,f(6) = f(7) = f(8) = f(9) = 1
发现有很多项结果都是一样的，于是我们可以打包计算这些块。
下面给一个更普遍的例题。
求 \sum_{i=1}^n \varphi(i)\lfloor\frac{n}{i}\rfloor,n \le 10^{6} 。
不加证明的给出一个定理：一个块的左端为 $l$，则右端点 $r$ 为 $\lfloor \dfrac{n}{\lfloor \frac{n}{l}\rfloor}\rfloor$。
因为一个块内的所有 $\lfloor\dfrac{n}{i}\rfloor$ 都是相同的，因此这个块对最终答案的贡献为：
$$(f(r) - f(l-1)) \times \lfloor\frac{n}{l}\rfloor$$
然后我们就可以 $O(\sqrt{n})$ 解决这个问题了。
### 例题
[UVA15526 H(n)](https://www.luogu.com.cn/problem/UVA11526)
#### 题意
求
$$\sum_{i=1}^n\,\lfloor\frac{n}{i}\rfloor$$
#### 思路
板子。对于每一个块 $[l,r]$，对答案的贡献为 $(r-l+1)\times \lfloor \frac{n}{l}\rfloor$，$O(\sqrt{n})$ 计算即可。
#### 代码
```cpp
#include<bits/stdc++.h>
using namespace std;
long long read()
{
long long s = 0,w = 1;
char ch = getchar();
while(ch < '0' || ch > '9')
{
if(ch == '-') w = -1;
ch = getchar();
}
while(ch >= '0' && ch <= '9')
{
s = (s << 1) + (s << 3) + (ch ^ 48);
ch = getchar();
}
return s * w;
}
long long T,n,ans;
int main()
{
//freopen("input.in","r",stdin);
T = read();
while(T--)
{
n = read(),ans = 0;
if(n == 0)
{
cout << 0 << endl;
continue;
}
for(long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
ans += (r - l + 1) * (n / l);
}
cout << ans << endl;
}
return 0;
}
```
## 乘法逆元
如果一个线性同余方程 $ax \equiv 1 \pmod b$，则 $x$ 称为 $a \bmod b$ 的逆元，记作 $a^{-1}$。
### 扩展欧几里得法
定睛一看，这不还是扩欧的板子吗？
设 $ax - 1 = by$，变换一下得到 $ax - by = 1$，扩展欧几里得求解即可，要求 $a,b$ 互质。
```cpp
void exgcd(int a,int b,int& x,int& y)
{
if (b == 0) {x = 1, y = 0;return;}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
```
### 快速幂法
只适用于 $b$ 为素数的情况。
因为 $ax \equiv 1 \pmod b$；
所以 $ax \equiv a^{b-1} \pmod b$（根据费马小定理）
所以 $x \equiv a^{b-2} \pmod b$。
然后我们就可以用快速幂来求了。
### 线性求逆元
转自 [OI-Wiki](https://oi-wiki.org/math/number-theory/inverse/) 。
#### 求出 $1,2,\dots,n$ 中每个数关于 $p$ 的逆元。
如果对于每个数进行单次求解，以上两种方法就显得慢了，很有可能超时，所以下面来讲一下如何 $O(n)$ 求逆元。
首先，很显然的 $1^{-1} \equiv 1 \pmod p$；
证明：对于 $\forall \,p \in \mathbf{Z}$，有 $1 \times 1 \equiv 1 \pmod p$ 恒成立，故在 $p$ 下 $1$ 的逆元是 $1$，而这是推算出其他情况的基础。
其次对于递归情况 $i^{-1}$，我们令
$k = \lfloor \frac{p}{i} \rfloor$，$j = p \bmod i$，有 $p = ki + j$。再放到 $\bmod \,p$ 意义下就会得到：$ki+j \equiv 0 \pmod p$；
两边同时乘 $i^{-1} \times j^{-1}$：
$$kj^{-1}+i^{-1}\equiv 0 \pmod p$$
$$i^{-1} \equiv -kj^{-1} \pmod p$$
再带入 $j = p \bmod i$，有 $p = ki + j$，有：
$$i^{-1} \equiv -\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1} \pmod p$$
我们注意到 $p \bmod i < i$，而在迭代中我们完全可以假设我们已经知道了所有的模 $p$ 下的逆元 $j^{-1}$, $j < i$。
故我们就可以推出逆元，利用递归的形式，而使用迭代实现：
$$
i^{-1} \equiv \begin{cases}
1, & \text{if } i = 1, \\
-\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1}, & \text{otherwise}.
\end{cases} \pmod p
$$
例题：[P3811【模板】模意义下的乘法逆元](https://www.luogu.com.cn/problem/P3811)
无脑直接套板子即可。
```cpp
#include<bits/stdc++.h>
using namespace std;
int read()
{
int s = 0,w = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') w = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') {s = (s << 1) + (s << 3) + (ch ^ 48);ch = getchar();}
return s * w;
}
const int N = 3e6 + 10;
long long n,p,inv[N];
int main()
{
//freopen("input.in","r",stdin);
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);n = read(),p = read();inv[1] = 1;
for(int i = 2;i <= n;i++) inv[i] = (long long)(p - p / i) * inv[p % i] % p;
for(int i = 1;i <= n;i++) cout << inv[i] << '\n';return 0;
}
```
#### 求出任意 $n$ 个数的逆元
上面的方法只能求 $1$ 到 $n$ 的逆元，如果需要求任意给定 $n$ 个数 $(1 \le a_i < p)$ 的逆元，就需要下面的方法：
首先计算 $n$ 个数的前缀积，记为 $s_i$，然后使用快速幂或扩展欧几里得法计算 $s_n$ 的逆元，记为 $sv_n$。
因为 $sv_n$ 是 $n$ 个数的积的逆元，所以当我们把它乘上 $a_n$ 时，就会和 $a_n$ 的逆元抵消，于是就得到了 $a_1$ 到 $a_{n-1}$ 的积逆元，记为 $sv_{n-1}$。
同理我们可以依次计算出所有的 $sv_i$，于是
$a_i^{-1}$ 就可以用 $s_{i-1} \times sv_i$ 求得。
所以我们就在 $O(n + \log p)$ 的时间内计算出了 $n$ 个数的逆元。
代码实现：
```cpp
s[0] = 1;
for (int i = 1;i <= n;i++) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2);
for (int i = n;i >= 1;i--) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1;i <= n;i++) inv[i] = sv[i] * s[i - 1] % p;
```
## 线性同余方程
懒得写了。转自 [OI-Wiki](https://oi-wiki.org/math/number-theory/linear-equation/) 。
### 定义
形如 $ax\equiv b\pmod n
的方程称为 线性同余方程。其中，a、b 和 n 为给定整数，x 为未知数。需要从区间 [0, n-1] 中求解 x ，当解不唯一时需要求出全体解。
用逆元求解
首先考虑简单的情况，当 a 和 n 互素时，即 \gcd(a, n) = 1 。
此时可以计算 a 的逆元，并将方程的两边乘以 a 的逆元，可以得到唯一解。
证明
x\equiv ba ^ {- 1} \pmod n
接下来考虑 a 和 n 不互素，即 \gcd(a, n) \ne 1 的情况。此时不一定有解。例如，2x\equiv 1\pmod 4 没有解。
设 g = \gcd(a, n) ，即 a 和 n 的最大公约数，其中 a 和 n 在本例中大于 1。
当 b 不能被 g 整除时无解。此时，对于任意的 x ，方程 ax\equiv b\pmod n 的左侧始终可被 g 整除，而右侧不可被 g 整除，因此无解。
如果 g 整除 b ，则通过将方程两边 a、b 和 n 除以 g ，得到一个新的方程：
a'x\equiv b' \pmod{n'}
其中 a' 和 n' 已经互素，这种情形已经解决，于是得到 x' 作为 x 的解。
很明显，x' 也将是原始方程的解。这不是唯一的解。可以看出，原始方程有如下 g 个解：
x_i\equiv (x'+ i\cdot n') \pmod n \quad \text{for } i = 0 \ldots g-1
总之，线性同余方程的 解的数量 等于 g = \gcd(a, n) 或等于 0 。
用扩欧求解
一眼板子。懒得说了。
int exgcd(int a, int b, int& x, int& y)
{
if (b == 0) {x = 1;y = 0;return a;}
int d = ex_gcd(b, a % b, x, y),temp = x;
x = y;y = temp - a / b * y;return d;
}
bool lieu(int a, int b, int c, int& x, int& y)
{
int d = exgcd(a, b, x, y);
if (c % d != 0) return 0;
int k = c / d;x *= k;y *= k;return 1;
}中国剩余定理
内容
中国剩余定理可求解如下形式的一元线性同余方程组（其中 n_1, n_2, \cdots, n_k 两两互质）：
\begin{cases}
x\equiv a_1 \pmod {n_1} \\
x\equiv a_2 \pmod {n_2} \\
\,\,\,\,\,\,\vdots \\
x\equiv a_k \pmod {n_k} \\
\end{cases}
过程
$2$. 对于第 $i$ 个方程：
- 计算 $m_i=\frac{n}{n_i}$，
- 计算 $m_i$ 在模 $n_i$ 意义下的逆元 $m_i^{-1}$；
计算 $c_i=m_im_i^{-1}$（不要对 $n_i$ 取模）。
$3$. 方程组在模 $n$ 意义下的唯一解为：
$x=\sum_{i=1}^k a_ic_i \pmod n$。
### 代码
```cpp
void exgcd(long long a,long long b,long long &x,long long &y)
{
if(b == 0)
{
x = 1;
y = 0;
return;
}
exgcd(b,a % b,x,y);
long long t = x;
x = y;
y = t - (a / b) * y;
return;
}
long long CRT(int k, long long* a, long long* r)
{
long long n = 1, ans = 0;
for (int i = 1; i <= k; i++) n = n * r[i];
for (int i = 1; i <= k; i++)
{
long long m = n / r[i], b, y;
exgcd(m,r[i],b,y);
ans = (ans + a[i] * m * b % n) % n;
}
return (ans % n + n) % n;
}
```
例题：[P2480【SDOI2010】古代猪文](https://www.luogu.com.cn/problem/P2480)
**数论全家桶**：扩展欧几里得算法 + 欧拉定理/费马小定理 + 中国剩余定理 + 卢卡斯定理。
观察题目，发现题目让我们求的式子是下面这个：
$$g^{\sum_{k \mid n} C_{n}^{\frac{n}{k}}} \bmod 999911659$$
对称的，我们将 $k$ 换为 $\dfrac{n}{k}$，式子就变成了一下这样：
$$g^{\sum_{k \mid n} C_{n}^{k}} \bmod 999911659$$
查一下发现 $999911659$ 是一个质数，所以我们可以利用费马小定理(欧拉定理)换一下得到这个式子：
$$g^{\sum_{k \mid n} C_{n}^{k} \bmod 999911658} \bmod 999911659$$
考虑计算 $\sum_{k \mid n} C_{n}^{k} \bmod 999911658$，剩余部分快速幂即可。
注意到 $n$ 的因数不会太多，考虑枚举每一个 $k$，计算对应的 $C_{n}^{k}$ 再求和。
对于组合数模一个数的计算，我们很容易想到卢卡斯定理，但是卢卡斯定理要求模数为质数，$999911658$ 显然不是质数，考虑两种方法：
- 扩展卢卡斯定理：~~太难了不会用~~ 杀鸡焉用牛刀
- 中国剩余定理
我们将 $999911658$ 分解为 $2 \times 3 \times 4679 \times 35617$，对于这 $4$ 个数分别用卢卡斯定理，假设得出的解分别为 $a_1,a_2,a_3,a_4$，用中国剩余定理合并答案即可。
$$
\begin{cases}
x\equiv a_1 \pmod {2} \\
x\equiv a_2 \pmod {3} \\
x\equiv a_3 \pmod {4679} \\
x\equiv a_4 \pmod {35617} \\
\end{cases}
$$
得出的唯一解是模 $2 \times 3 \times 4679 \times 35617$ 即 $999911658$ 的值，就是要求的答案，最后快速幂求答案即可。
代码：
```cpp
#include <bits/stdc++.h>
using namespace std;
long long read()
{
long long s = 0,w = 1;
char ch = getchar();
while(ch < '0' || ch > '9')
{
if(ch == '-') w = -1;
ch = getchar();
}
while(ch >= '0' && ch <= '9')
{
s = (s << 1) + (s << 3) + (ch ^ 48);
ch = getchar();
}
return s * w;
}
const int N = 5e4 + 10,p = 999911658;
long long q[N],r[N],x[N],b[5];
long long T,n,m,g,f[N];
long long multiple = 1,ans;
void init(long long p)
{
f[0] = 1;
for(int i = 1;i <= N;i++) f[i] = (f[i - 1] * i) % p;
}
long long exgcd(long long a,long long b,long long &x,long long &y)
{
long long t;
if (b == 0)
{
x = 1;
y = 0;
return a;
}
else
{
t = exgcd(b,a % b,x,y);
int xx = x;
x = y;
y = xx - a / b * y;
return t;
}
}
long long quickpow(long long c,long long n,long long m)
{
long long res = 1;
while(n)
{
if(n & 1) res = res * c % m;
c = c * c % m;
n >>= 1;
}
return res;
}
long long C(long long n,long long m,long long p)
{
if(m > n)return 0;
return((f[n] * quickpow(f[m],p - 2,p)) % p * quickpow(f[n - m],p - 2,p) % p);
}
long long Lucas(long long n, long long m, long long p)
{
if (m == 0) return 1;
return (C(n % p, m % p, p) * Lucas(n / p, m / p, p)) % p;
}
int a[5] = {0,2,3,4679,35617};
long long CRT()
{
for(int i = 1;i <= 4;i++)
{
q[i] = a[i];
multiple = multiple * a[i];
}
for(int i = 1;i <= 4;i++)
{
r[i] = multiple / q[i];
long long x = 0,y = 0;
exgcd(r[i],q[i],x,y);
if (x < 0) ans += b[i] * r[i] * (x + q[i]);
else ans += b[i] * r[i] * x;
}
return ans % multiple;
}
int main()
{
//freopen("input.in","r",stdin);
n = read(),g = read();
long long gg = g;
if(!(g % 999911659)){cout << 0 << endl;return 0;}
for(int i = 1;i <= 4;i++)
{
init(a[i]);
for(int j = 1;j * j <= n;j++)
{
if(n % j == 0)
{
b[i] = (b[i] + Lucas(n,j,a[i])) % a[i];
if(j * j == n) continue;
b[i] = (b[i] + Lucas(n,n / j,a[i])) % a[i];
}
}
}
cout << quickpow(gg,CRT(),p + 1) << endl;
return 0;
}
```
## 扩展中国剩余定理
### 内容
中国剩余定理可求解如下形式的一元线性同余方程组：
$$
\begin{cases}
x\equiv a_1 \pmod {n_1} \\
x\equiv a_2 \pmod {n_2} \\
\,\,\,\,\,\,\vdots \\
x\equiv a_k \pmod {n_k} \\
\end{cases}
$$
与中国剩余定理不同的是，这里不保证 $n_1,n_2,\cdots,n_k$ 两两互质。
### 过程
$1$. 取出方程组的前两个方程。
设两个方程分别是 $x\equiv a_1 \pmod {m_1}$、$x\equiv a_2 \pmod {m_2}$；
将它们转化为不定方程：$x=m_1p+a_1=m_2q+a_2$，其中 $p, q$ 是整数，则有 $m_1p-m_2q=a_2-a_1$。
由裴蜀定理，当 $a_2-a_1$ 不能被 $\gcd(m_1,m_2)$ 整除时，无解；
其他情况下，可以通过扩展欧几里得算法解出来一组可行解 $(p, q)$；
则原来的两方程组成的模方程组的解为 $x\equiv b\pmod M$，其中 $b=m_1p+a_1$，$M=\text{lcm}(m_1, m_2)$。
$2$. 将新方程放入方程组的第一个位置，删去原先的两个方程。
$3$. 重复以上操作，直到只剩下一个方程。
$4$. 记这个方程为 $x \equiv a \pmod n$，输出 $(a \bmod n + n) \bmod n$ 即可。
### 代码
```cpp
#include<bits/stdc++.h>
using namespace std;
__int128 read()
{
__int128 s = 0,w = 1;
char ch = getchar();
while(ch < '0' || ch > '9')
{
if(ch == '-') w = -1;
ch = getchar();
}
while(ch >= '0' && ch <= '9')
{
s = (s << 1) + (s << 3) + (ch ^ 48);
ch = getchar();
}
return s * w;
}
void write(__int128 a)
{
if(a > 9) write(a / 10);
putchar(a % 10 + '0');
}
void exgcd(__int128 a,__int128 b,__int128 &x,__int128 &y)
{
if(b == 0)
{
x = 1;
y = 0;
return;
}
exgcd(b,a % b,x,y);
__int128 t = x;
x = y;
y = t - (a / b) * y;
return;
}
__int128 gcd(__int128 a,__int128 b){return __gcd(a,b);}
__int128 lcm(__int128 a,__int128 b){return a * b / gcd(a,b);}
struct EQUATION
{
__int128 a,b;
}e[100100],ans;
int main()
{
//freopen("input.in","r",stdin);
__int128 n = read();
for(__int128 i = 1;i <= n;i++) e[i].a = read(),e[i].b = read();ans = e[1];
for(__int128 i = 2;i <= n;i++)
{
__int128 g = gcd(ans.a,e[i].a),x,y;
exgcd(ans.a / g,e[i].a / g,x,y);
EQUATION newans;
newans.b = (ans.b + (e[i].b - ans.b) / g * x * ans.a) % lcm(ans.a,e[i].a);
newans.a = lcm(ans.a,e[i].a);
ans = newans;
}
write((ans.b % ans.a + ans.a) % ans.a);
return 0;
}
```
例题：[P4777【模板】扩展中国剩余定理（EXCRT）](https://www.luogu.com.cn/problem/P4777)
板子题，直接写板子即可。注意需要 $\texttt{int128}$。
## 升幂引理
见 [OI-Wiki-升幂引理](https://oi-wiki.org/math/number-theory/lift-the-exponent/) 。
## 威尔逊定理
见 [OI-Wiki-威尔逊定理](https://oi-wiki.org/math/number-theory/wilson/) 。
对于素数 $p$ 有 $(p-1)!\equiv -1\pmod p$。
## 卢卡斯定理
部分转自 [OI-Wiki](https://oi-wiki.org/math/number-theory/lucas/) 。
Lucas 定理用于求解大组合数取模的问题，其中模数必须为素数。正常的组合数运算可以通过递推公式求解但当问题规模很大，而模数是一个不大的质数的时候，就不能简单地通过递推求解来得到答案，需要用到 Lucas 定理。
### 内容
Lucas 定理内容如下：对于质数 $p$，有
$$\binom{n}{m}\bmod p = \binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}\cdot\binom{n\bmod p}{m\bmod p}\bmod p$$
观察上述表达式，可知 $n\bmod p 和 m\bmod p$ 一定是小于 $p$ 的数，可以直接求解，
$\displaystyle\binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}$ 可以继续用 Lucas 定理求解。这也就要求 $p$ 的范围不能够太大，一般在 $10^5$ 左右。边界条件：当 $m=0$ 的时候，返回 $1$。
时间复杂度为 $O(f(p) + g(n)\log n)$，其中 $f(n)$ 为预处理组合数的复杂度，$g(n)$ 为单次求组合数的复杂度。
### 代码实现
```cpp
long long Lucas(long long n, long long m, long long p)
{
if (m == 0) return 1;
return (C(n % p, m % p, p) * Lucas(n / p, m / p, p)) % p;
}
```
### 证明
不会证。转自 [OI-Wiki](https://oi-wiki.org/math/number-theory/lucas/) 。
考虑
$\displaystyle\binom{p}{n} \bmod p$ 的取值，
注意到
$\displaystyle\binom{p}{n} = \frac{p!}{n!(p-n)!}$，分子的质因子分解中 $p$ 的次数恰好为 $1$，因此只有当 $n = 0$ 或 $n = p$ 的时候 $n!(p-n)!$ 的质因子分解中含有 $p$，因此
$\displaystyle\binom{p}{n} \bmod p = [n = 0 \vee n = p]$。
进而我们可以得出
$$
\begin{aligned}
(a+b)^p &= \sum_{n=0}^p \binom pn a^n b^{p-n}\\
&\equiv \sum_{n=0}^p [n=0\vee n=p] a^n b^{p-n}\\
&\equiv a^p + b^p \pmod p
\end{aligned}
$$
注意过程中没有用到费马小定理，因此这一推导不仅适用于整数，亦适用于多项式。因此我们可以考虑二项式 $f^p(x)=(ax^n + bx^m)^p \bmod p$ 的结果
$$
\begin{aligned}
(ax^n + bx^m)^p &\equiv a^p x^{pn} + b^p x^{pm} \\
&\equiv ax^{pn} + bx^{pm}\\
&\equiv f(x^p)
\end{aligned}
$$
考虑二项式 $(1+x)^n \bmod p$，那么
$\displaystyle\binom n m$ 就是求其在 $x^m$ 次项的取值。使用上述引理，我们可以得到
$$
\begin{aligned}
(1+x)^n &\equiv (1+x)^{p\lfloor n/p \rfloor} (1+x)^{n\bmod p}\\
&\equiv (1+x^p)^{\lfloor n/p \rfloor} (1+x)^{n\bmod p}
\end{aligned}
$$
注意前者只有在 $p$ 的倍数位置才有取值，而后者最高次项为 $n\bmod p \le p-1$，因此这两部分的卷积在任何一个位置只有最多一种方式贡献取值，即在前者部分取 $p$ 的倍数次项，后者部分取剩余项，即
$\displaystyle\binom{n}{m}\bmod p = \binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}\cdot\binom{n\bmod p}{m\bmod p}\bmod p。
扩展卢卡斯定理
还是转自 OI-Wiki 。算法能理解，证明看不懂。菜死了。
Lucas 定理中对于模数 p 要求必须为素数，那么对于 p 不是素数的情况，就需要用到 exLucas 定理。
实现
第一部分：中国剩余定理
要求计算二项式系数
考虑利用中国剩余定理合并答案，这种情况下我们只需求出
$\binom{n}{m}\bmod p^\alpha$ 的值即可（其中 $p$ 为素数且 $\alpha$ 为正整数）。
根据 唯一分解定理，将 $M$ 质因数分解：
$$M{p_1}^{\alpha_1}\cdot{p_2}^{\alpha_2}\cdots{p_r}^{\alpha_r}=\prod_{i=1}^{r}{p_i}^{\alpha_i}$$
对于任意 $i,j$，有 ${p_i}^{\alpha_i}$ 与 ${p_j}^{\alpha_j}$ 互质，所以可以构造如下 $r$ 个同余方程：
$$
\left\{
\begin{aligned}
a_1\equiv \displaystyle\binom{n}{m}&\pmod {{p_1}^{\alpha_1}}\\
a_2\equiv \displaystyle\binom{n}{m}&\pmod {{p_2}^{\alpha_2}}\\
&\cdots\\
a_r\equiv \displaystyle\binom{n}{m}&\pmod {{p_r}^{\alpha_r}}\\
\end{aligned}
\right.
$$
我们发现，在求出 $a_i$ 后，就可以用中国剩余定理求解出
$\displaystyle\binom{n}{m}$。
第二部分：移除分子分母中的素数
根据同余的定义，
$\displaystyle a_i=\binom{n}{m}\bmod {p_i}^{\alpha_i}$，问题转化成，求
$\displaystyle \binom{n}{m} \bmod p^\alpha$（$p$ 为质数）的值。
根据组合数定义
$\displaystyle \binom{n}{m} = \frac{n!}{m! (n-m)!}$，
$\displaystyle \binom{n}{m} \bmod p^\alpha = \frac{n!}{m! (n-m)!} \bmod p^\alpha$。
由于式子是在模 $p^\alpha$ 意义下，所以分母要算乘法逆元。
同余方程 $ax \equiv 1 \pmod p$（即乘法逆元）有解 的充要条件为 $\gcd(a,p)=1$（裴蜀定理），
然而 无法保证有解，发现无法直接求 $\operatorname{inv}_{m!}$ 和 $\operatorname{inv}_{(n-m)!}$，
所以将原式转化为：
$$\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z} \bmod p^\alpha$$
$x$ 表示 $n!$ 中包含多少个 $p$ 因子，$y$, $z$ 同理。
第三部分：Wilson 定理的推论
问题转化成，求形如：
$\frac{n!}{q^x}\bmod q^k
的值。这时可以利用 Wilson 定理的推论。
即：
$\displaystyle \prod_{i,(i,q)=1}^{q^k}i$ 一共循环了
$\displaystyle \lfloor\frac{n}{q^k}\rfloor$ 次，暴力求出
$\displaystyle \prod_{i,(i,q)=1}^{q^k}i$，然后用快速幂求
$\displaystyle \lfloor\frac{n}{q^k}\rfloor$ 次幂。
最后要乘上
$\displaystyle \prod_{i,(i,q)=1}^{n \bmod q^k}i$，显然长度小于 $q^k$，暴力乘上去。
上述三部分乘积为 $n!$。最终要求的是
$$\frac{n!}{q^x}\bmod{q^k}$$
所以有：
$$
n! = q^{\left\lfloor\frac{n}{q}\right\rfloor} \cdot \left(\left\lfloor\frac{n}{q}\right\rfloor\right)! \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right)$$
于是：
$$
\frac{n!}{q^{\left\lfloor\frac{n}{q}\right\rfloor}} = \left(\left\lfloor\frac{n}{q}\right\rfloor\right)! \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right)
$$
$\displaystyle \left(\left\lfloor\frac{n}{q}\right\rfloor\right)!$ 同样是一个数的阶乘，所以也可以分为上述三个部分，于是可以递归求解。
等式的右边两项不含素数 $q$。事实上，如果直接把 $n$ 的阶乘中所有 $q$ 的幂都拿出来，等式右边的阶乘也照做，这个等式可以直接写成：
$$
\frac{n!}{q^{x}} = \frac{\left(\left\lfloor\frac{n}{q}\right\rfloor\right)!}{q^{x'}} \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right)
$$
式中的 $x$ 和 $x'$ 都表示把分子中所有的素数 $q$ 都拿出来。改写成这样，每一项就完全不含 $q$ 了。
递归的结果，三个部分中，左边部分随着递归结束而自然消失，中间部分可以利用 Wilson 定理的推论 $0$，右边部分就是推论 2 中的 $\prod_{j\geq 0}(N_j!)_p$。
下面这种写法，拥有单次询问 $O(p\log p)$ 的时间复杂度。
### 代码
代码实现：
```cpp
long long calc(long long n, long long x, long long P)
{
if (!n) return 1;
long long s = 1;
for (long long i = 1; i <= P; i++) if (i % x) s = s * i % P;
s = Pow(s, n / P, P);
for (long long i = n / P * P + 1; i <= n; i++) if (i % x) s = i % P * s % P;
return s * calc(n / x, x, P) % P;
}
long long multilucas(long long m, long long n, long long x, long long P)
{
int cnt = 0;
for (long long i = m; i; i /= x) cnt += i / x;
for (long long i = n; i; i /= x) cnt -= i / x;
for (long long i = m - n; i; i /= x) cnt -= i / x;
return Pow(x, cnt, P) % P * calc(m, x, P) % P * inverse(calc(n, x, P), P) %P * inverse(calc(m - n, x, P), P) % P;
}
long long exlucas(long long m, long long n, long long P)
{
int cnt = 0;
long long p[20], a[20];
for (long long i = 2; i * i <= P; i++)
{
if (P % i == 0)
{
p[++cnt] = 1;
while (P % i == 0) p[cnt] = p[cnt] * i, P /= i;
a[cnt] = multilucas(m, n, i, p[cnt]);
}
}
if (P > 1) p[++cnt] = P, a[cnt] = multilucas(m, n, P, P);
return CRT(cnt, a, p);
}
```
上一篇笔记：[Part 1](https://www.luogu.com.cn/article/jmt0skay) 。
下一篇笔记：[Part 3](https://www.luogu.com.cn/article/lycvbvlm) 。