初等数论学习笔记
A2_Zenith
·
2023-10-13 22:10:50
·
算法·理论
本文将整理一些提高组数论知识以及一些在提高组可能会用到的数论知识（如数论分块）。
\
\
一些基础知识点
这里我们明确一些基础概念和基础的知识点：
$\mathbb{N}$：全体自然数。
$\mathbb{Z^+},\mathbb{N}^+$：全体正整数。
$\mathbb{P}$：我不知道这么表述严不严谨，但是我一般拿它当全体素数来用。
+ $a$ 是 $b$ 的倍数，当且仅当存在 $k \in \mathbb{Z},a=kb$。同理，$b$ 是 $a$ 的因数（约数），当且仅当 存在$k \in \mathbb{Z},a=kb$。这两种情况都可记作 $b\mid a$。
+ $\lfloor x\rfloor$ 表示不大于 $x$ 的最大整数。同理，$\lceil x\rceil$ 表示不小于 $x$ 的最小整数。
+ $a \bmod b$ 表示 $a$ 除以 $b$ 的余数。
+ $a\equiv b \pmod p$，当且仅当 $a \bmod p=b \bmod p$。
+ $\gcd(a,b)$ 为 $a,b$ 的最大公约数，亦即最大的 $k$ 且 $k$ 同时为 $a,b$ 的因数。$\text{lcm}(a,b)$ 表示其最小公倍数，亦即最小的 $k$ 且 $k$ 同时为 $a,b$ 的倍数。
这里我们顺便再普及一下基础的定理和算法吧。
### 唯一分解定理
任意大于 $1$ 的自然数 $n$ 必然能被表示成有限个素数的幂次的乘积，并且这种表示是唯一的。即 $n = \prod\limits_{i=1}^k{p_i^{k_i}}$。这就是数论中的唯一分解定理。
### $\gcd$ 和辗转相除法
$\gcd(a,b) = \gcd(a,b-a)$。证明不难。由此可以推出 $\gcd(a,b)=\gcd(a,b\bmod a)=\gcd(b\bmod a,a)$。很简单，不断减，减到不能再减为止，其实就是取模了。这样算 $\gcd$ 是 $\log$ 的，因为每次要递归的子问题数据规模至少减少一半。（ $b\bmod x < \frac{b}{2}$ ）。
### 素数判定
对于单个数 $n$ 判定其是否为素数，我们目前涉及的算法都是 $\sqrt n$ 的。因为一个数若有因数，那我们可以将其两两配对为形如 $(x,\frac{n}{x})$ 的形式。不难看出我们只需检测其中的较小的一个是否能被 $n$ 整除。这样的话，我们肯定希望 $x \leq \frac{n}{x}$。解得 $x \leq \sqrt n$。也就是说我们只需要对 $2$ 到 $\sqrt n$ 进行检测。
### 筛法
#### 埃氏筛
有时对于 $1$ 到 $n$ 的所有整数，我们希望能够快速的找出其中的素数。对于每个数进行素性检测显然很慢，$n \sqrt n$。 一个朴素的想法是：对于所有的正整数 $i$，我们能够知道：$ki$ 一定为合数，其中 $k>2$。那么我们不妨对这些 $kn$ 打上标记。这样对于所有 $i$ ，都要进行 $\frac{n}{i}$ 轮打标记。这其实是 $n \log n$ 的。
当然这还是太慢，你发现只需对每个素数 $p$ 进行 $\frac{n}{p}$ 轮打标记。有人问：你怎么知道这个数是不是素数？很简单，如果你发现这个数在它前面的数打标记时没有被标记为合数，那这个数一定是素数，这很显然。这大概是 $n \log \log n$ 的。这其实就是埃氏筛了。这里我相信聪明的大家一定能想出更多优化，这里听神 Gapple 说埃氏筛配合 bitset 比欧拉筛快，这里 %。
#### 欧拉筛
当然埃氏筛还有优化空间，因为一个合数可能会被标记多次。
这里我们只让每个合数被它的最小的素因子筛。那这样就是 $\mathcal{O}(n)$ 的了。这很棒。
那我们直接考虑从枚举 $i$ 变为枚举 $k$。这里我们先放代码：
```cpp
bool tag[maxn];
int p[maxn];
int cnt=0;
void init(int n){
for(int k = 2; k <= n; k++){
if(!tag[k]){
p[++cnt] = k;
}
for(int j = 1; j <= cnt && k * p[j] <= n; j++){
tag[k * p[j]]=1;
if(k % p[j] == 0)break;
}
}
}
```
这里为了让大家看的清楚加了空格。
注意这一行：
```if(k % p[j] == 0)break;```
当 $k\bmod p_j=0$ 时，也就是说明 $p_j$ 是最小的能整除 $k$ 的素数（因为你存下 $p_j$ 显然是从小往大了存的）。也就是说，若是继续打标记的话，后面的 $k \times p_l$ 一定都有素因子 $p_j$，也就是说它们都被 $p_j$ 筛过了，在筛一遍没有意义，所以直接 break 就好。
## Euler 函数
- 欧拉函数
$$\LARGE\varphi$$
读作 $\operatorname{/fai/}$。
先给出一些性质：
- $\varphi(n)$ 表示小于等于 $n$ 的与 $n$ 互质的数的个数。
- $p \in \mathbb{P}\Rightarrow \varphi(p)=p-1$。
- 设 $n = \prod\limits_{i=1}^k{p_i^{k_i}}$，则 $\varphi(n) = n\prod\limits_{i=1}^k(1-\frac{1}{p_i})$ 。特别的，$\varphi(1)=1$。
每个 $1-\frac{1}{p_i}$ 都删掉了与 $n$ 共享 $p_i$ 这个素因子的所有数。
由此我们容易发现 $\varphi$ 的积性。
- $n = \sum\limits_{d|n}\varphi(d)$(即欧拉反演）。
证：因为 $\varphi$ 是积性函数，所以我们只需证 $\sum\limits_{d|p^k}\varphi(p^k)=p^k$。
$\sum\limits_{d|p^k}\varphi(d)=1+\sum\limits_{i=1}^k(p-1)p^{i-1}
=1+(p-1)\sum\limits_{i=0}^{k-1}p^i
=1+(p-1)\times \frac{p^k-1}{p-1}
而 $\forall n$，我们令 $n=\prod\limits_{i=1}^mp_i^{c_i}
=\prod\limits_{i=1}^m\sum\limits_{d|p_i^{c_i}}\varphi(p_i^{c_i})
到了这一步，你发现这实际上这和枚举 n 的因子的 \varphi 是等价的。
下面 CODETIME。
单点求 \varphi 代码实现：
int phi(int n){
if(n==1)return 1;
int ans=n;
for(int i=2;i*i<=n;i++){
if(n%i==0){
ans=ans*(i-1)/i;
while(n%i==0)n/=i;
}
}
if(n>1)ans=ans*(n-1)/n;
return ans;
}线性筛 \varphi ：
vector<int>a;
bool isp[10000007];
int pp[10000007];
void sieve(int n){
pp[1]=1;
for(int i=1;i<=n;i++){
isp[i]=1;
}
for(int i=2;i<=n;i++){
if(isp[i]){
a.push_back(i);
pp[i]=i-1;
}
for(int j=0;j<a.size()&&a[j]*i<=n;j++){
isp[a[j]*i]=0;
if(i%a[j]==0){
pp[i*a[j]]=pp[i]*a[j];
break;
}
pp[i*a[j]]=pp[i]*pp[a[j]];
}
}
}我们在 i,p_j 互质时使用积性函数的性质来筛 \varphi(i\times p_j) ，在 i,p_j 不互质时，你发现 i\times p_j 相较于 i ，其本质不同质因子的情况没有改变，也就是 n\prod\limits_{i=1}^k(1-\frac{1}{p_i}) 后面的那个 Product 没有改变，改变的只有前面的那个 n ，它乘上了一个 p_j 。
下面是题目。
$\
\
GCDSUM
求
{\sum\limits_{x=1}^n}{\sum\limits_{y=1}^n}\gcd(x,y)
$${\sum\limits_{x=1}^n}{\sum\limits_{y=1}^n}\gcd(x,y)$$
$$=\sum\limits_{d=1}^{n}{\sum\limits_{x=1}^{n}}{\sum\limits_{y=1}^{n}}d[\gcd(x,y)=d]$$
$$=\sum\limits_{d=1}^{n}d{\sum\limits_{x=1}^{n}}{\sum\limits_{y=1}^{n}}[\gcd(x,y)=d]$$
$$=\sum\limits_{d=1}^{n}d{\sum\limits_{x=1}^{\lfloor{n/d}\rfloor}}{\sum\limits_{y=1}^{\lfloor{n/d}\rfloor}}[\gcd(x,y)=1]$$
$$=\sum\limits_{d=1}^{n}d(2*{(\sum\limits_{i=1}^{\lfloor{n/d}\rfloor}\varphi(i))-1)}$$
这一步是因为：
你枚举所有 $i\leq j$，产生了 $\sum\limits_{i=1}^{n/d}\varphi(i)$ 的贡献。对于 $i\geq j$ 同理。但是 $i=j=1$ 的情况被计算了两次，所以要减一。
线性筛筛 $\varphi$，求其前缀和。枚举 $d$ ，后面的东西有前缀和可以 $\mathcal{O}(1)$ 求。
总复杂度 $\mathcal{O}(n)$。
```cpp
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<string>
#include<cstring>
#include<queue>
#include<stack>
#include<cstdlib>
#include<iomanip>
#include<map>
#define int long long
#define double long double
#define lc(p) p<<1
#define rc(p) p<<1|1
#define pii pair<int,int>
using namespace std;
int a[100007];
bool isp[100007];
int phi[100007]={0,1};
int sum[100007];
int cnt=0;
void sieve(int n){
for(int i=1;i<=n;i++){
isp[i]=1;
}
for(int i=2;i<=n;i++){
if(isp[i]){
a[cnt]=i;
phi[i]=i-1;
cnt++;
}
for(int j=0;j<cnt&&a[j]*i<=n;j++){
isp[a[j]*i]=0;
if(i%a[j]==0){
phi[i*a[j]]=phi[i]*a[j];
break;
}
phi[i*a[j]]=phi[i]*phi[a[j]];
}
}
}
signed main(){ios::sync_with_stdio(0);
int n;
cin>>n;
sieve(n);
for(int i=1;i<=n;i++){
//cout<<phi[i]<<endl;
sum[i]=sum[i-1]+phi[i];
}
int ans=0;
for(int i=1;i<=n;i++){
ans+=((2*sum[n/i]-1)*i);
}
cout<<ans;
}
```
当然，这里也可以用我们前文提到过的欧拉反演实现。
$$\sum\limits_{x=1}^n\sum\limits_{y=1}^n\gcd(x,y)$$
$$=\sum\limits_{x=1}^n\sum\limits_{y=1}^n\sum\limits_{d|\gcd(x,y)}\varphi(d)$$
$$=\sum\limits_{x=1}^n\sum\limits_{y=1}^n\sum\limits_{d|x,d|y}\varphi(d)$$
$$=\sum\limits_{d=1}^n\varphi(d)\lfloor\frac{n}{d}\rfloor^2$$
比上面那个式子好看多了是不是？
并且这个式子可以用**数论分块**进行优化而到达 $\mathcal{O}(\sqrt n)$ 的复杂度。
这版代码（$\mathcal{O}(n)$）留作练习。
$\
\
GCD
求
\sum\limits_{p\in \mathbb{P}}{\sum\limits_{x=1}^n}{\sum\limits_{y=1}^n}[\gcd(x,y)=p]
$$\sum\limits_{p\in prime}{\sum\limits_{x=1}^n}{\sum\limits_{y=1}^n}[\gcd(x,y)=p]$$
$$=\sum\limits_{p\in prime}{\sum\limits_{x=1}^{\lfloor{n/p}\rfloor}}{\sum\limits_{y=1}^{\lfloor{n/p}\rfloor}}[\gcd(x,y)=1]$$
$$=\sum\limits_{p\in prime}(2*{(\sum\limits_{i=1}^{\lfloor{n/p}\rfloor}\varphi(i)})-1)$$
之后同上。
其实这个题用莫反做柿子会好看一点？
```cpp
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<string>
#include<cstring>
#include<queue>
#include<stack>
#include<cstdlib>
#include<iomanip>
#include<map>
#define int long long
#define double long double
#define lc(p) p<<1
#define rc(p) p<<1
#define pii pair<int,int>
using namespace std;
int phi(int n){
if(n==1)return 1;
int ans=n;
for(int i=2;i*i<=n;i++){
if(n%i==0){
ans=ans*(i-1)/i;
while(n%i==0)n/=i;
}
}
if(n>1)ans=ans*(n-1)/n;
return ans;
}
vector<int>a;
bool isp[10000007];
int p[10000007];
int sum[10000007];
void sieve(int n){
for(int i=1;i<=n;i++){
isp[i]=1;
}
for(int i=2;i<=n;i++){
if(isp[i]){
a.push_back(i);
p[i]=i-1;
}
for(int j=0;j<a.size()&&a[j]*i<=n;j++){
isp[a[j]*i]=0;unsigned
if(i%a[j]==0){
p[i*a[j]]=p[i]*a[j];
break;
}
p[i*a[j]]=p[i]*p[a[j]];
}
}
}
signed main(){ios::sync_with_stdio(0);
int n;
cin>>n;
sieve(n);
p[1]=1;
for(int i=1;i<=n;i++){
sum[i]=sum[i-1]+p[i];
}
int ans=0;
for(int i=0;i<a.size();i++){
ans=ans+2*sum[n/a[i]]-1;
}
cout<<ans<<endl;
return 0;
}
```
$\
\
Euler 定理（数论）
然后我们对欧拉定理进行一个证明。
我们设 R=\{r_1,r_2,r_3\cdots r_{\varphi(n)}\} ，R 中 每一个数都与 n 互质且模 n 不同余。
设 S=\{ar_1,ar_2,ar_3\cdots ar_{\varphi(n)}\} 。
显然 S 中的每个数模 n 互不同余且都与 n 互质。
你发现 S=R 。
然后 \prod\limits_{i=1}^{\varphi(n)}r_i\equiv\prod\limits_{i=1}^{\varphi(n)}ar_i\equiv a^{\varphi(n)}\prod\limits_{i=1}^{\varphi(n)}r_i(\bmod m) 。
两边的 \prod\limits_{i=1}^{\varphi(n)}r_i 与 n 显然互质，所以我们可以将其约去。
最后得到 a^{\varphi(n)}\equiv 1(\bmod m) 。
扩展欧拉定理：
a^b\equiv\begin{cases}a^b,b < \varphi(m)\\a^{{b \operatorname{mod} \varphi(m) }+\varphi(m)},b\geq\varphi(m)\end{cases}(\operatorname{mod} m)
这两个定理主要用于降幂，将一个叠叠乐降成人类可以计算的程度。
上帝与集合的正确用法
2^{2^{2^{2^{\cdots}}}} \operatorname{mod} P
最终会保持在同一个值，求这个值。
叠叠乐就是像这样的东西。
多测，T \leq 10^3 ，P \leq 10^7 。
这东西乍一看肥肠的不可做啊。你考虑用欧拉定理降一手幂。
2^{2^{2^{2^{\cdots}}}} \operatorname{mod} P
=2^{(2^{2^{2^{\cdots}}}) \operatorname{mod} \varphi(P)+\varphi(P)} \operatorname{mod} P
前一部分又变成了一个类似于原式的柿子！
但是模数在变小！
我们定义 Slv(P)=2^{2^{2^{2^{\cdots}}}} \operatorname{mod} P ，则 Slv(P)=2^{Slv(\varphi(P))+\varphi(P)} 。
当然了，这里不需要对扩展欧拉定理进行分讨，因为 2^{2^{2^{2^{\cdots}}}} 远大于 P 。
对一个数取很多次 \varphi ，变成 1 大概只需要 \log 次。
于是：
```cpp
int slv(int p){
if(p==1)return 0;
else return qpow(2ll,slv(pp[p])+pp[p],p);
}
```
关于只会递归 $\log$ 层的证明：
若 $P \operatorname{mod} 2 ==0$，由上文欧拉函数的计算公式，$P$ 至少会被砍一半；
否则，$\varphi(P) = \prod\frac{\text{偶数}}{\text{奇数}}$，又变成偶数了。
$\mathcal{O}(T \log P)$。
```cpp
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<string>
#include<cstring>
#include<queue>
#include<stack>
#include<cstdlib>
#include<iomanip>
#include<map>
#define int long long
#define double long double
#define lc(p) p<<1
#define rc(p) p<<1|1
#define pii pair<int,int>
using namespace std;
int phi(int n){
if(n==1)return 1;
int ans=n;
for(int i=2;i*i<=n;i++){
if(n%i==0){
ans=ans*(i-1)/i;
while(n%i==0)n/=i;
}
}
if(n>1)ans=ans*(n-1)/n;
return ans;
}
vector<int>a;
bool isp[10000007];
int pp[10000007];
int sum[10000007];
void sieve(int n){
pp[1]=1;
for(int i=1;i<=n;i++){
isp[i]=1;
}
for(int i=2;i<=n;i++){
if(isp[i]){
a.push_back(i);
pp[i]=i-1;
}
for(int j=0;j<a.size()&&a[j]*i<=n;j++){
isp[a[j]*i]=0;
if(i%a[j]==0){
pp[i*a[j]]=pp[i]*a[j];
break;
}
pp[i*a[j]]=pp[i]*pp[a[j]];
}
}
}
int qpow(int a,int b,int p){
if(b==1)return a%p;
int ans=1;
if(b&1){
int u=qpow(a,b>>1,p);
ans=ans*u%p;
ans=ans*u%p;
ans=ans*a%p;
return ans;
}
else{
int u=qpow(a,b>>1,p);
ans=ans*u%p;
ans=ans*u%p;
return ans;
}
}
int slv(int p){
if(p==1)return 0;
else return qpow(2ll,slv(pp[p])+pp[p],p);
}
signed main(){
int t;
cin>>t;
sieve(10000000);
while(t--){
int p;
cin>>p;
// for(int i=1;i<=mod;i++){
// cout<<pp[i]<<" ";
// }
cout<<slv(p)<<endl;
}
}
```
### Power Tower
*```2700```
呜呜呜怎么掉 ```2700``` 了。
给定 $n$ 与 长为 $n$ 的序列 $a$ ，给定模数 $P$ ，$Q$ 次询问，每次询问给定 $l$，$r$，求
$$ a_l^{a_{l+1}^{a_{l+2}^{\cdots^{a_r}}}} \operatorname{mod} P$$
这里，我们令 $b={a_{l+1}^{a_{l+2}^{\cdots^{a_r}}}}$，则答案 $=a^b\bmod P$。
$n,q\leq 10^5
P \leq 10^9
和上面其实是同理的。
可写出如下递归函数：
int qpow(int a,int b,int p){
int ans=1;
while(b>0){
if(b&1){
ans=ans*a;
if(ans>=p)ans=ans%p+p;
}
a=a*a;
if(a>=p)a=a%p+p;
b>>=1;
}
return ans;
}
int slv(int l,int r,int p){
if(l==r)return (a[l]<p?a[l]:a[l]%p+p);
if(p==1)return 1;
return qpow(a[l],slv(l+1,r,phi(p)),p);
}注意快速幂与递归函数，回想扩展欧拉定理：
a^b\equiv\begin{cases}a^b,b < \varphi(m)\\a^{{b \operatorname{mod} \varphi(m) }+\varphi(m)},b\geq\varphi(m)\end{cases}(\operatorname{mod} m)
快速幂与递归返回行要分讨。
#include<bits/stdc++.h>
#define int long long
using namespace std;
//LONG LIVE BRUTE-FORCE ALGORITHM!!!
int n;
int p;
int a[100007];
int q;
map<int,int> mp;
int phi(int x){
if(mp.count(x))return mp[x];
int ans=x;
int t=x;
for(int i=2;i*i<=x;i++){
if(!(t%i)){
ans=ans/i*(i-1);
while(!(t%i))t/=i;
}
}
if(t>1)ans=ans/t*(t-1);
mp[x]=ans;
return ans;
}
int qpow(int a,int b,int p){
int ans=1;
while(b>0){
if(b&1){
ans=ans*a;
if(ans>=p)ans=ans%p+p;
}
a=a*a;
if(a>=p)a=a%p+p;
b>>=1;
}
return ans;
}
int slv(int l,int r,int p){
if(l==r)return (a[l]<p?a[l]:a[l]%p+p);
if(p==1)return 1;
return qpow(a[l],slv(l+1,r,phi(p)),p);
}
signed main()
{
cin>>n;
cin>>p;
for(int i=1;i<=n;i++){
cin>>a[i];
}
cin>>q;
while(q--){
int l,r;
cin>>l>>r;
cout<<slv(l,r,p)%p<<endl;;
}
return 0;
}[Ynoi2016] 炸脖龙 I
上个题套个树状数组即可。注意实现细节。
\
\
逆元
数 a 在 \operatorname{mod} p 意义下的逆元定义为
ax \equiv 1 (\operatorname{mod} p)
的解，记作 a^{-1} \pmod p 。
显然，当 \gcd(a,m) = 1 时 a 才在模 m 意义下有逆元。
当 a 在模意义下的逆元都可以表示成 km +r 的形式，其中 k 是整数。
简证：
取 a 任意两个逆元。
有
ax \equiv 1(\operatorname{mod} m),
ax' \equiv 1(\operatorname{mod} m).
则
a(x-x') \equiv 0(\operatorname{mod} m).
因为 \gcd(a,m) = 1 ，故 m|x-x' ，x \equiv x' ，证毕。
问题来了：我们怎么求呢？
我会费小！
注意到当模数为质数时，a^{P-1} \equiv 1(\operatorname{mod} P) ，故 a \times a^{P-2} \equiv 1(\operatorname{mod} P) ,所以我们知道：a^{P-2} \operatorname{mod} P 是 a 在模 P 意义下的一个逆元。
这是 \log 的。
我会欧拉定理！
这是根号的。
那如果我要让你求 $1$ 到 $n$ 在模 $P$ 意义下的逆元，$n \leq 10^7$，$P \leq 10^8$，阁下又将如何应对？
这时我们需要一种**线性**的算法。
我们要求 $i^{-1} \bmod P$。
显然$1^{-1} \equiv 1$。
令 $p$（模数） $ = k \times i + r (0 < r \leq i)$，则
$$ki + r \equiv 0(\operatorname{mod} p)$$
两边同时乘$i^{-1}r^{-1} \mod P$，则
$$ kr^{-1}+i^{-1} \equiv 0$$
$$i^{-1} \equiv -kr^{-1}$$
$$ i^{-1} \equiv -\lfloor{\frac{p}{i}}\rfloor (p \operatorname{mod} i)^{-1}$$
然后你发现 $i^{-1}\bmod P$ 可以从前面递推而来。
于是写下如下代码：
```cpp
inv[x] = (p-p/i) * inv[p%x];
```
逆元 1 就是这个递推（费小亦可），逆元 2 稍微推一下式子，做一些转化即可。
## EXgcd
用于求方程
$$ax+by=\gcd(a,b) { \ \ \ \ \ \ \ \ \ (1)}$$
的一组特解。
令 $\gcd(a,b)=c$。
考虑方程
$$bx'+(a \operatorname{mod} b)y'=c {\ \ \ \ \ \ \ \ \ (2)}$$
假设我们有 $x'$，$y'$，则
$$bx'+(a - \lfloor{\frac{a}{b}}\rfloor \times b)y'= c$$
$$b(x'-\lfloor{\frac{a}{b}}\rfloor y')+ay'=c$$
则我们有一组 $x$，$y$ 满足
$$\begin{cases}x = y'\\y = x' - \lfloor{\frac{a}{b}}\rfloor y'\end{cases}$$
$(2)$ 可以继续往下递归。
当 $y$ 系数为 $0$ 时，这时是
$$ax=\gcd(a,0)$$
，有特解 $x = 1$，$y = 0$。
所以有代码
```cpp
void exgcd(int a,int b,int &x,int &y){
if(!b){x=1;y=0;return;}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
```
应用：
$\log$ 求单个数逆元。不要求 $P$ 是质数。
比欧拉定理好多了不是？
$$ax \equiv 1(\operatorname{mod} p)$$
则存在 $k$ ，使得
$$ax+py=1$$
用 $\operatorname{exgcd}$ 算一下，求出来的 $x$ 就是 $a$ 的逆元。
$\rm exgcd$ 板题我不建议大家写，细节很多很麻烦。建议大家写 P1516【青蛙的约会】。
## CRT
亦可称作孙子定理。
用于求解一组形如
$$
\begin{cases}x \equiv a_1\pmod {m_1}\\
x \equiv a_2\pmod {m_2}\\
x \equiv a_3\pmod {m_3}\\
\cdots\cdots\\
x \equiv a_n\pmod {m_n}\end{cases}
$$
的线性同余方程组的最小非负整数解，其中 $m$ 两两互质。
~~这是个构造题~~
我们尝试构造一组特解：
令 $M = \prod\limits_{i=1}^nm_i$，$mk_i=\frac{M}{m_i},inv_i=mk_i^{-1}(\operatorname{mod} m_i)$。
则我们肯定有解 $x'=\sum\limits_{i=1}^na_i \times mk_i\times inv_i$。
对于第 $i$ 个同余方程，$\forall j \neq i$，$a_j\times mk_j\times inv_j \equiv 0 (\operatorname{mod} m_i)$。
原因很简单，$mk_j$ 显然有因数 $m_i$。
而 $a_i\times mk_i\times inv_i \operatorname{mod} m_i$ 显然等于 $a_i$。
所以 $\forall 1\leq i \leq n $，$x' \operatorname{mod} m_i = a_i$ 。
$x' \operatorname{mod} M$ 就是我们要的最小非负整数解。
原因很简单，$\forall 1\leq i \leq n$，$M \operatorname{mod} m_i = 0$。答案肯定是 $x'+kM$。最小非负整数解显然是 $x' \operatorname{mod} M$ 啊。
### 【模板】
```cpp
#include<bits/stdc++.h>
#define int long long
using namespace std;
void exgcd(int a,int b,int &x,int &y){
if(!b){x=1;y=0;return;}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
int a[15],b[15];
int M[15];
int MProd=1;
int n;
signed main(){
cin>>n;
for(int i=1;i<=n;i++){cin>>a[i]>>b[i];MProd*=a[i];}
for(int i=1;i<=n;i++){M[i]=MProd/a[i];}
int sum=0;
for(int i=1;i<=n;i++){
int inv,y;
exgcd(M[i],a[i],inv,y);
inv=(inv%a[i]+a[i])%a[i];
sum=(sum+b[i]*M[i]*inv)%MProd;
}
cout<<sum<<endl;
}
```
有些时候我们做题，很多题要你把答案对一个数取模。你求出了这个问题的解，但是毒瘤出题人给定的模数不是质数，怎么办呢，你把那个模数拆成好多个质数的幂次相乘，对于每个素因子幂次算出模它们意义下的答案，用 CRT 合并。
扩展 $\rm Lucas$ 定理就是一个很好的例子。~~再比如说任意模数 NTT~~
## Lucas
### 前情提要
相信大家一定知道组合数了，$C_n^m$ 就是在 $n$ 个球里选出 $m$ 个球的方案数，不考虑顺序。
由定义，我们知道 $C_n^m=\frac{A_n^m}{A_m^m}=\frac{n!}{m!(n-m)!}
然后是二项式定理：
(x+y)^n=\sum\limits_{i=0}^nC_n^ix^iy^{n-i}
原因还是肥肠的简单的，你考虑 x^iy^{n-i} ，这就意味着你要在 n 个 (x+y) 中选出 i 个来贡献 x ，剩下的贡献 y ，显然不考虑顺序。
正文
Lucas 定理用于快速求组合数取模，其中模数比较小。
你在看别人博客时有时会看到 \dbinom{n}{m} ，它 =C_n^m 。
好然后我们给出 Lucas 定理的结论。
结论：p 是质数，则 C_n^m \operatorname{mod} p=C_{n/p}^{m/p}C_{n \operatorname{mod} p}^{m\operatorname{mod}p} \operatorname{mod} p 。
证明：
(1+x)^p = \sum\limits_{i=0}^pC_p^ix^i
有 C_n^m=\frac{n!}{m!(n-m)!} ，当 n 为 p ，m \in (0,p) 时，分子显然没有能约掉上面的 p 的（因为 p 是质数！），所以 \forall 0 < i < p ，C_p^i \equiv 0(\operatorname{mod} p) 。
所以上面那一坨 \operatorname{mod} p 肯定只剩一个 (1 + x^p) 。
我们设 n = kp + r (0 \leq r < p ) ，m = tp + j(0 \leq j < p ) 则
(1+x)^n (\operatorname{mod }p)
=(1+x)^{kp}(1+x)^{r} (\operatorname{mod }p)
=((1+x)^p)^k(1+x)^r(\operatorname{mod }p)
=(1+x^p)^k(1+x)^r(\operatorname{mod }p)
观察 x^m 项系数。则：
C_n^mx^m \equiv C_k^tx^{tp}C_r^jx^j(\operatorname{mod }p)
C_n^mx^m \equiv C_k^tC_r^jx^{tp+j}(\operatorname{mod }p)
C_n^mx^m \equiv C_k^tC_r^jx^{m}(\operatorname{mod }p)
\therefore C_n^m \equiv C_k^tC_r^j (\operatorname{mod }p)
得证！
（这里其实是用到了生成函数的思想，x 只是一个形式符号，一个位置的标识符。）
后面的那一坨只能爆算了，但前面的那一坨可以通过 Lucas 加速！
于是，我们写出如下代码：
int Lucas(int n,int m,int p){
if(m==0)return 1;
return Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}我们设 n 在 p 进制下的表示是 n_kn_{k-1}n_{k-2}\cdots n_0 ，
$$ C_n^m \operatorname{mod} p$$
$$= C_{n_kn_{k-1}\cdots n_0}^{m_km_{k-1}\cdots m_0} \operatorname{mod} p$$
$$ =C_{n_kn_{k-1}\cdots n_1}^{m_km_{k-1}\cdots m_1} C_{n_0}^{m_0}\operatorname{mod} p$$
$$ =C_{n_kn_{k-1}\cdots n_2}^{m_km_{k-1}\cdots m_2} C_{n_1}^{m_1}C_{n_0}^{m_0}\operatorname{mod} p$$
$$ = \cdots $$
$${ = C_{n_k}^{m_k}C_{n_{k-1}}^{m_{k-1}} \cdots C_{n_0}^{m_0} \operatorname{mod} p}$$
由此我们得出一个结论：
$C_n^m \operatorname{mod }p$ 等于 $n$，$m$ 在 $p$ 进制下，每一位求一次组合数乘起来的值。
由此我们便有两个推论：
- $\exists 0 \leq i \leq k,m_i>n_i \iff C_n^m \equiv 0 (\operatorname{mod} p) $。
- $C_n^m$ 为奇数的充要条件是 $n \operatorname{bitand} m = m$。
第二个推论也能从上面推过来。
例题时间。
### 《瞿葩的数字游戏》T3-三角圣地
宝宝题。
你需要构造一个 $1 - n$ 的排列。
像这样：
$$1\ 3\ 4\ 2$$
相邻两个数向下合并，像这样：
$$1\ 3\ 4\ 2$$
$$4\ 7\ 6 $$
这样，每一层都向下合并，最后只剩一个数，最大化这个数，$\operatorname{mod} 10007$ 并输出。
$10007$ 是质数。
$n \leq 10^6$。
你把两个数合并想象成上面一层的每个数都向左下，右下做一次贡献，则最上面一层的一个数对答案的贡献就是在一个三角形图上走向顶点的路径数量。
最上面一层的第 $m$ 个位置对答案的贡献显然是 $C_n^m$ 。
贪心的想，肯定是把大的数放在中间更有，小的一边站着去。
注意到 $ C_n^m = C_n ^ {n-m}$，所以肯定 $1$，$2$ 配 $C_n^0(C_n^n)$，$3$，$4$ 配 $C_n^1(C_n^{n-1})$ 像这样。
所以答案其实是
$$\large\sum\limits_{i=1}^ni \times C_n^{(i-1)/2}$$
$n$ 很大，我们使用 Lucas 加速。
代码：
```cpp
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e4+7;
int fac[10100];
int ifac[10100];
int qpow(int a,int b){
a%=mod;
int ans=1;
for(int t=b;t;t>>=1,a=a*a%mod)if(t&1)ans=ans*a%mod;
return ans;
}
int inv(int x){
return qpow(x,mod-2);
}
int C(int n,int m){
if(m>n)return 0;
return fac[n]%mod*ifac[m]%mod*ifac[n-m]%mod;
}
int Lucas(int n,int m){
if(m==0)return 1;
return Lucas(n/mod,m/mod)*C(n%mod,m%mod)%mod;
}
int n;
int sum=0;
signed main(){
cin>>n;
ifac[0]=fac[0]=1;
for(int i=1;i<mod;i++){
fac[i]=fac[i-1]*i%mod;
ifac[i]=inv(fac[i]);
}
for(int i=1;i<=n;i++){
sum=(sum+i*Lucas(n-1,i-1>>1))%mod;
}
cout<<sum<<endl;
}
```
### [CTSC2017] 吉夫特
结论题，有结论就简单了。
输入一个长度为 $n$ 的数列 $a_1, a_2, \cdots , a_n$ 问有多少个长度大于等于 $2$ 的不上升的子序列满足：
$$\prod _{i=2}^{k} \binom{a_{b_{i-1}}}{a_{b_i}} \bmod 2 = \binom{a_{b_1}}{a_{b_2}} \times \binom{a_{b_2}}{a_{b_3}} \times \cdots \binom{a_{b_{k-1}}}{a_{b_k}} \bmod 2 > 0$$
输出这个个数对 $10^9+7$ 取模的结果。
$1\leq n\leq 211985$，$1\leq a_i\leq 233333$。所有的 $a_i$ 互不相同。
这不是 ```Lucas``` 定理第二个推论嘛。
肯定是 $\forall 1 \leq i < k$，$a_{b_i} \operatorname{bitand} a_{b_{i+1}} = a_{b_{i+1}}$。
你都是上一个数的子集了，那单调性还考虑啥呀！
这值域，小的很嘛！
直接设 $f_i$ 为以 $i$ **这个数** 为结尾的子序列个数刷表转移一下就行。
```cpp
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<string>
#include<cstring>
#include<queue>
#include<stack>
#include<cstdlib>
#include<iomanip>
#include<map>
#define int long long
#define db long double
#define pii pair<int,int>
#define up(i,l,r) for(int i=(l);i<=(r);i++)
#define down(i,l,r) for(int i=(l);i>=(r);--i)
#define p_b push_back
#define m_p make_pair
using namespace std;
//LONG LIVE BRUTE-FORCE ALGORITHM!!!
const int maxn=250007;
const int mod=1e9+7;
int f[maxn];
int n;
int a[maxn];
signed main(){
cin>>n;
for(int i=1;i<=n;i++)cin>>a[i];
int ans=0;
for(int i=1;i<=n;i++){
int t=a[i];
//cout<<f[t]<<endl;;
for(int m=t&(t-1);m;m=m-1&t)f[m]+=f[t]+1;
/*
t的每个子集（即n&m==m的每个m自己都是一个子序列），所以要+1。
*/
ans=(ans+f[t])%mod;
}
cout<<ans<<endl;
}
```
我们补充一下扩展 Lucas。
当模数不是质数时，$(1+x)^p \equiv 1+x^p ( \bmod p)$ 这东西就寄了。
怎么办呢。
我们先考虑模数是素数的幂次怎么办，这样我们就能够通过 CRT 合并答案了。
现在问题变成了 $\dbinom{n}{m} \bmod p^k$。
众所周知，$\dbinom{n}{m}=\large\frac{n!}{m!(n-m)!}
这时聪明的你或许就会想到：直接对分母求逆元就好了啊！
先别急嘛。
阶乘这个东西不一定和模数互质的。
所以，我们直接暴力拆开：
\Large\frac{n!}{m!(n-m)!}=\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}\normalsize p^{x-y-z}
我们扔掉了阶乘中所有的 p ，这样就可以求逆元了。
这个非常朴素吧。
然后就是大头了：
首先，我们推柿子：
\ n! \bmod p^k
=(p\times 2p\times3p\cdots\lfloor\frac{n}{p}\rfloor p)\times\prod_{i=1,p\nmid i}^ni\bmod p^k
先把所有的带 p 的项都拎出来。
=p^{\lfloor\frac{n}{p}\rfloor}\times\lfloor\frac{n}{p}\rfloor!\times\prod_{i=1,p\nmid i}^ni\bmod p^k
这里发现一堆 p 能提出来，剩下的是阶乘。
=p^{\lfloor\frac{n}{p}\rfloor}\times\lfloor\frac{n}{p}\rfloor!\times\large(\prod_{i=1,p\nmid i}^{p^k}i)^{\lfloor\frac{n}{p^k}\rfloor}\normalsize\times\large\prod\limits_{i=\lfloor\frac{n}{p^k}\rfloor p^k+1,p\nmid i}^ni
发现剩下的东西 \bmod p^k 以 p^k 为一个循环节，所以我们把循环节拎出来（就是第三段），剩下的余项是第四段。
所以我们如果要求消除 p 影响的阶乘（我们设其为f(n) ）的话，我们就把那个 p^{\lfloor\frac{n}{p}\rfloor} 扔掉。就是：
=f(\lfloor\frac{n}{p}\rfloor)\times\large(\prod_{i=1,p\nmid i}^{p^k}i)^{\lfloor\frac{n}{p^k}\rfloor}\normalsize\times\large\prod\limits_{i=\lfloor\frac{n}{p^k}\rfloor p^k+1,p\nmid i}^ni
第三，四段暴力计算，剩下的那个阶乘是一个规模更小的子问题，可以递归计算。
所以啊，每层要计算 p^k 次，要向下递归 \log 层，这是 p^k \log P 的。
边界是 f(0)=1 。
最后一个问题：我们最后还要把 p 的贡献加回来。
这东西我们其实可以由上面的柿子中直观的看出来。设 g(n) 为 n! 中所含 p 的幂次，我们很显然的有：
g(n)=\lfloor\frac{n}{p}\rfloor + g(\lfloor\frac{n}{p}\rfloor)
即
\sum\limits_{i=0}^\infty\lfloor\frac{n}{p^i}\rfloor
这其实就是大家所熟知的 Legendre 公式。即使你并没有听说过它的大名，也一定在小学奥数时就对其有所耳闻。
所有问题都解决了。
上代码。
#include<bits/stdc++.h>
#define int long long
using namespace std;
inline int qpow(int a,int b,int p){
int ans=1;
while(b){
if(b&1)ans=ans*a%p;
a=a*a%p;
b>>=1;
}
return ans;
}
void exgcd(int a,int b,int &x,int &y){
if(!b){x=1,y=0;return;}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
inline int inv(int a,int p){
a%=p;
int x=0,y=0;
exgcd(a,p,x,y);
return (x%p+p)%p;
}
int pcf(int n,int p){
if(n<p)return 0;
return (n/p)+pcf(n/p,p);
}
int fac(int n,int p,int pc){
if(n==0)return 1;
int BL=1;
int T=(n/pc);
int ans=(T&1?pc-1:1);
int rem=1;
for(int i=n/pc*pc+1;i<=n;i++)if(i%p)rem=rem*(i%pc)%pc;
return ans*rem%pc*fac(n/p,p,pc)%pc;
}
int C_prt(int n,int m,int p,int pc){
if(n<m)return 0;
int fz=fac(n,p,pc);
int fm=fac(n-m,p,pc)*fac(m,p,pc)%pc;
return fz*inv(fm,pc)%pc*qpow(p,pcf(n,p)-pcf(n-m,p)-pcf(m,p),pc)%pc;
}
int C(int n,int m,int mod){
static int p[30];
static int pc[30],cnt;
static int val[30];
static int M[30];
static int Mi[30];
int k=mod;
for(int i=2;i*i<=mod;i++){
if(k%i==0){
++cnt;
p[cnt]=i;
pc[cnt]=1;
while(k%i==0){
pc[cnt]*=i;
k/=i;
}
}
}
if(k>1)p[++cnt]=k,pc[cnt]=k;
for(int i=1;i<=cnt;i++)M[i]=mod/pc[i],Mi[i]=inv(M[i],pc[i]);
for(int i=1;i<=cnt;i++)val[i]=C_prt(n,m,p[i],pc[i]);
int ans=0;
for(int i=1;i<=cnt;i++){
ans=(ans+val[i]*M[i]%mod*Mi[i]%mod)%mod;
}
return ans;
}
int n,m,p;
signed main(){ios::sync_with_stdio(0);cin.tie(0);
cin>>n>>m>>p;
cout<<C(n,m,p)<<endl;
}这里大家可能会注意到这一行：
int ans=(T&1?pc-1:1);
这里是对上文中那个大柿子中的第三项的计算。
你发现
\large(\prod_{i=1,p\nmid i}^{p^k}i)^{\lfloor\frac{n}{p^k}\rfloor}
这东西中，除了 1 和 p^k-1 以外，每个数的逆元在这个乘积式中，且一个数和它的逆元可以两两配对。毕竟我们把所有带 p 的都扔掉了，所以剩下的每个数肯定都存在逆元。两两配对后乘起来得到一堆 1 ，忽略掉，剩下的是 1 和 p-1 ，乘起来是 p-1 即 -1 。
这个似乎可以优化算法常数（？
\
\
\
\
上点强度。
[SHOI2015] 超能粒子炮-改
多测，每次询问给定 n,k ，求 \sum\limits_{i=0}^kC_n^i \bmod 2333 。( 2333 是质数）。
T\leq 10^5,n,k\leq \textcolor{#FF0000}{\large 10^{18}}
很有脑子的一个题。
直接考虑与 \log 有关的做法。
考虑设 f(n,k) 为 \sum\limits_{i=0}^kC_n^i \bmod 2333 。
下文令 p=2333 。
我们直接套一个 \rm Lucas ：
\sum\limits_{i=0}^kC_n^i \bmod 2333
=\sum\limits_{i=1}^n\textcolor{#0000FF}{C_{n/p}^{i/p}} \times \textcolor{#FF0000}{C_{n \bmod p}^{i\bmod p}}
你发现右边的那一坨的值（标红的部分）随 i 的改变，以 p 为周期。并且对于 t:0\rightarrow \lfloor\frac{k}{p}\rfloor-1 ，每个 \textcolor{#0000FF}{C_{n/p}^t} （这里的 t 对应着原来的 i/p ）都对应着一个 \sum\limits_{i=0}^{p-1}C_{n\bmod p}^i 。
写的更具体些就是这样：
\large\sum\limits_{i=0}^{\lfloor\frac{k}{p}\rfloor-1}\textcolor{#0000ff}{C_{n/p}^i}\sum\limits_{j=i\times p}^{(i+1)\times p-1}\textcolor{#ff0000}{C_{n\bmod p}^{j \bmod p}}
剩下的，对于 i:\lfloor\frac{k}{p}\rfloor\times p\rightarrow k ，你发现他们它们对应的蓝色部分都是 C_{n/p}^{k/p} ，红色部分是 \sum\limits_{i=0}^{k\bmod p}C_{n\bmod p}^i 。
所以总的这个柿子就被我们拆成了
\sum\limits_{t=0}^{k/p-1}C_{n/p}^t\sum\limits_{j=0}^{p-1}C_{n\bmod p}^j + C_{n/p}^{k/p}\sum\limits_{i=0}^{k\bmod p}C_{n\bmod p}^i
=\sum\limits_{t=0}^{k/p-1}C_{n/p}^t\times f(n \bmod p,p-1)+C_{n/p}^{k/p}f(n\bmod p,k\bmod p)
=f(n \bmod p,p-1)\sum\limits_{t=0}^{k/p-1}C_{n/p}^t+C_{n/p}^{k/p}f(n\bmod p,k\bmod p)
=f(n \bmod p,p-1)\sum\limits_{t=0}^{k/p-1}C_{n/p}^t+C_{n/p}^{k/p}f(n\bmod p,k\bmod p)
=f(n \bmod p,p-1)\times f(n/p,k/p-1)+C_{n/p}^{k/p}f(n\bmod p,k\bmod p)
呐，这东西在 n,k\leq p 的情况可以 \mathcal{O}(p^2) 预处理出来，剩下的你发现这是一个递归状物，并且子问题的规模每次都是原来的 \large\frac{1}{p} 。这说明什么？这说明这个东西只会向下递归 \log_pn 次。那个组合数可以 \rm Lucas 做到 \log 。
总复杂度就是 \mathcal{O}(p^2+T\log^2n) ，\log 的底数是 2333 ，通过此题绰绰有余。
数论分块
也叫做整除分块。
用于快速求解形如 \sum\limits_{i=1}^nf(i)\lfloor{\frac{n}{i}}\rfloor 的和式，其中能够 \mathcal{O}(1) 计算 f 的前缀和。
当 n = 30 时，\lfloor{\frac{n}{i}}\rfloor 的值如下：
1\ \ 30
2\ \ 25
3\ \ 10
4\ \ 7
5\ \ 6
6\ \ 5
7\ \ 4
8 - 10\ \ 3
11 - 15\ \ 2
16 - 30\ \ 1
我们发现，后面 \lfloor{\frac{n}{i}}\rfloor 的值出现了块状！
经过瞪眼观察法，你发现：\lfloor{\frac{n}{i}}\rfloor 的取值至多只有 2\lfloor{\sqrt{n}}\rfloor 种！
简证：
对于 1 \leq i \leq \lfloor{\sqrt{n}}\rfloor ，\lfloor{\frac{n}{i}}\rfloor 有至多 \lfloor{\sqrt{n}}\rfloor 种取值。
对于 \lfloor{\sqrt{n}}\rfloor \leq i \leq n ，\lfloor{\frac{n}{i}}\rfloor \leq \lfloor{\sqrt{n}}\rfloor ，也有至多 \lfloor{\sqrt{n}}\rfloor 种取值！
得证！
所以，我们只要把值相同的 \lfloor{\frac{n}{i}}\rfloor 打包计算即可。
问题现在变成：使得 \lfloor{\frac{n}{i}}\rfloor = \lfloor{\frac{n}{j}}\rfloor 的最大的 j 。
\lfloor{\frac{n}{i}}\rfloor = \lfloor{\frac{n}{j}}\rfloor \leq \frac{n}{j}
j \leq \frac{n}{\lfloor{\frac{n}{i}}\rfloor}
所以最大的整数 j 是
\lfloor\frac{n}{\lfloor{\frac{n}{i}}\rfloor}\rfloor
所以我们可以写出如下代码：
int H(int n){
int ans=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
int t=f[r]-f[l-1];
ans+=t*(n/l);
}
return ans;
}
数论分块是莫反与其他数论推柿子题中的重要技巧。比如说，我们把上面的 GCD SUM 那个题变为多测，$T\leq 10^5,n\leq 2\times10^5$，这时我们所讲的第一个方法就寄了，但是使用第二个方法，使用数论分块将单次询问的时间复杂度优化到 $\sqrt n$ 就可以平安通过。
完结撒花~