数论学习笔记
TH911
·
·
算法·理论
大概率更好的阅读体验
本文可能会不定期更新。若存在内容正确性上的问题，请联系修复。
:::info[更新/修正记录]
本文中出现的代数式、数，如无特殊说明，均为整数。
对于 OI 中的数论题，答案不对请先尝试：
#define int __int128
typedef long long ll;
#define ll __int128
数论函数指的是定义域为正整数的函数，也可以视为一个序列。
基础部分
对于一个序列 $f_n=af_{n-1}+bf_{n-2}$，若 $\gcd(a,b)=1$，则有：
$$
\gcd(f_x,f_y)=f_{\gcd(x,y)}
$$
# 质数
## Fermat 素性测试
由费马小定理，对于质数 $p$，若 $\gcd(a,p)=1$，满足 $a^{p-1}\equiv 1\pmod p$，$a^{p-2}\equiv a^{-1}\pmod p$。
但是 $a^{p-1}\equiv1\pmod p$ 并不能推出 $p$ 为质数。
因此可以随机选择来测试。
## Miller–Rabin 素性测试
Miller-Rabin 素性测试，取质数集合 $A=\lbrace{2,3,5,7,11,13,17,19,23,29,31,37}\rbrace$，则可以通过随机化确定 `long long` 范围（$[0,2^{64})$）内的任意整数 $n$ 是否为质数。
从 $A$ 中取一底数 $a$，若：
* $n=a$，$n$ 为质数。
* $n\bmod a=0\land n>a$，$n$ 不为质数。
在都不成立的情况下，进行 Miller-Rabin 素性测试。
将 $n-1$ 转化：
$$
n-1=d\times2^r
$$
其中，$d,r\in\N^*$，也就是说，$d$ 是 $n-1$ 在二进制位上的一个前缀，满足前缀的后缀不为 $0$。
### 二次探测定理
$$
x^2\equiv1\pmod p\\
\Downarrow\\
x\equiv\pm1\pmod p
$$
使用平方差公式可证。
### Miller-Rabin 素性测试的实现
基于 $a$ 执行 $r$ 轮测试，通过二次探测定理判断。
### 参考代码
```cpp
constexpr const ll A[12]={2,3,5,7,11,13,17,19,23,29,31,37};
bool MillerRabin(ll n){
if(n<=1){
return false;
}
for(ll a:A){
if(n==a){
return true;
}
if(n%a==0){
return false;
}
bool no=true;
ll d=n-1,r=0;
while(!(d&1)){
r++;
d>>=1;
}
ll pl=qpow(a,d,n);
if(pl==1){
no=false;
}
for(int i=0;no&&i<r;i++){
if(pl==n-1){
no=false;
}
pl=(__int128)pl*pl%n;
}
if(no){
return false;
}
}
return true;
}
```
# 最大公约数
## 裴蜀定理
对于 $\forall a,b\in \Z$：
* $\exist x,y\in \Z$，使 $ax+by=\gcd(a,b)$。
* $\forall x,y\in\Z$，$\gcd(a,b)\mid(ax+by)$。
### 逆定理
$\forall a,b\in \Z$，若 $d>0$ 且 $d$ 为 $a,b$ 的公因数，且存在 $x,y$，使得 $ax+by=d$，则：
$$
d=\gcd(a,b)
$$
# 乘法逆元
若存在 $ax\equiv 1\pmod p$，则称 $x$ 为 $a$ 关于 $p$ 的逆元，记 $x=a^{-1}$。
当 $\gcd(a,p)=1$ 的时候（$p$ 为质数）逆元**一定存在**。
## 费马小定理
对于质数 $p$，若 $\gcd(a,p)=1$，满足 $a^{p-1}\equiv 1\pmod p$，$a^{p-2}\equiv a^{-1}\pmod p$。
常用于分数取模，即：
$$
\frac{a}{b}\equiv ab^{p-2}\pmod p
$$
详见[费马小定理](https://www.cnblogs.com/TH911/p/-/Fermat-Little-Therorem)。
## 欧拉定理
:::info[欧拉函数]
令 $n=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_k^{c_k}$，满足 $p_i$ 为质数，$c_i>0$。
令 $n=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_k^{c_k}$，满足 $p_i$ 为质数，$c_i>0$。
则：
$$
\varphi(n)=n\left(1-\frac 1{p_1}\right)\left(1-\frac 1{p_2}\right)\left(1-\frac 1{p_3}\right)\cdots\left(1-\frac 1{p_k}\right)
$$
推导见[容斥原理求解欧拉函数](https://www.cnblogs.com/TH911/p/-/Combinatorics#求解欧拉函数-varphin)。（组合数学没搬，放的原 Blog 链接）
:::
若 $\gcd(a,p)=1$，则 $a^{\varphi(p)}\equiv1\pmod p$。
显然，当 $p$ 为质数时，有 $\varphi(p)=p-1$，即费马小定理。因此，**费马小定理为欧拉定理的特殊形式**。
## 扩展欧拉定理
对于任意的 $a,p,b\in \N^*$，若满足 $\gcd(a,p)>1,b\geq \varphi(p)$，有：
$$
a^b\equiv a^{b\bmod \varphi(p)+\varphi(p)}\pmod p
$$
当 $\gcd(a,p)=1$ 时即欧拉定理。
可以使用初等证明或群论证明。
常用于降幂。
## 扩展欧几里得算法 exgcd
**注意不是“类欧几里得算法”。**
用于求关于 $x,y$ 的不定方程 $ax+by=\gcd(a,b)$ 的**一组整数特解**。
不妨令 $a>b$。
考虑到 $a=\left\lfloor\dfrac{a}{b}\right\rfloor b+a\bmod b$，代入可得：
$$
\left(\left\lfloor\frac{a}{b}\right\rfloor b+a\bmod b\right)x+by=\gcd(a,b)
$$
因此：
$$
b\left(\left\lfloor\frac ab\right\rfloor x+y\right)+(a\bmod b)x=\gcd(b,a\bmod b)
$$
令 $a'=b,x'=\left\lfloor\dfrac ab\right\rfloor x+y,b'=a\bmod b,y'=x$，有：
$$
a'x'+b'y'=\gcd(a',b')
$$
显然，可以**递归求解**。
那么直到 $b=0$ 时，可以直接解得一组整数特解 $\begin{cases}x=1\\y=0\end{cases}$。
设最终递归出来的解为 $\begin{cases}x=x_0\\y=y_0\end{cases}$。
令 $\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert，\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert$，则对于 $\forall k\in \N^*$，方程存在一组解为 $\begin{cases}x=x_0+k\Delta b\\y=y_0-k\Delta a\end{cases}$。
### 参考代码
```cpp
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
```
### exgcd 求乘法逆元
如果需要讨论 $a$ 模 $p$ 意义下的乘法逆元 $a^{-1}$，显然有 $\gcd(a,p)=1$。
那么可列方程 $ax+py=1$，显然 $ax\equiv1\pmod p$。
则通过扩展欧几里得算法求出 $x$，即求出了 $a$ 的逆元。
> [例题链接：[NOIP 2012] 同余方程](https://www.luogu.com.cn/problem/P1082)
求出 $x$ 后，为了使 $x$ 为正整数，只需要让 $x\leftarrow(x\bmod p+p)\bmod p$ 即可。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int a,b;
cin>>a>>b;
int x0,y0;
exgcd(a,x0,b,y0);
x0%=b;
if(x0<0){
x0+=b;
}
cout<<x0<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
### 例题：正整数最大/小解
> [例题链接](https://www.luogu.com.cn/problem/P5656)
>
> 给定不定方程 $ax+by=c$，$a,b,c\in \N^*$，判断其是否有解。
>
> 若有正整数解，分别求 $x,y$ 的最大/小值。
>
> 否则，求 $x,y$ 得最小正整数解。
对于无解，即 $\gcd(a,b)\nmid c$，很好判断。
先通过 exgcd 求出一组特解 $(x_0,y_0)$，则 $\begin{cases}x=x_0+k\Delta b\\y=y_0-k\Delta a\end{cases}$ 也是原方程的解。$\Delta b$ 为 $b$ 的单次变化量，最小为 $\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert$，$\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert$。
先求 $x$ 的最小正整数解 $x_{\min}$，即 $x_0+k\Delta b$ 最小，显然 $x_{\min}=x_0\bmod \Delta b$。
但是为了保证 $x_{\min}$ 为正整数，因此当 $x_0\bmod \Delta b\leq0$ 时，$x_{\min}=x_0\bmod \Delta b+\Delta b$。
$a,b>0$，因此 $x$ 最小时 $y$ 最大，因此可以确定 $y_{\max}=\dfrac{c-ax_{\min}}{b}$。
同理，可求出 $x_{\max},y_{\min}$。
当 $x_{\min}$ 为**正整数**，此时若 $y_{\max}>0$，则存在正整数解 $\begin{cases}x=x_{\min}\\y=y_{\max}\end{cases}$，因此可以判断正整数解（判断 $x_{\max}>0$ 也行）。
正整数解的个数即：
$$
\frac{x_{\max}-x_{\min}}{\Delta b}+1
$$
因为从 $x_{\min}$ 开始，$x_{\min}+k\Delta b$ 对应一组 $x,y$。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
#define int ll
typedef long long ll;
int gcd(int a,int b){
while(b){
int tmp=a;
a=b;
b=tmp%b;
}
return a;
}
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int T;
cin>>T;
while(T--){
int a,b,c;
cin>>a>>b>>c;
int gcdAB=gcd(a,b);
if(c%gcdAB){
cout<<"-1\n";
}else{
int w=c/gcdAB;
int x0,y0;
exgcd(a,x0,b,y0);
x0*=w,y0*=w;
// cerr<<"x0="<<x0<<" y0="<<y0<<endl;
int xMin,xMax,yMin,yMax;
int deltaA=a/gcdAB;
int deltaB=b/gcdAB;
xMin=x0%deltaB;
if(xMin<=0){
xMin+=deltaB;
}
yMax=(c-a*xMin)/b;
yMin=y0%deltaA;
if(yMin<=0){
yMin+=deltaA;
}
xMax=(c-b*yMin)/a;
// cerr<<"A="<<a<<" B="<<b<<endl;
// cerr<<"ΔA="<<deltaA<<" ΔB="<<deltaB<<endl;
// cerr<<"xMax="<<xMax<<" xMin="<<xMin<<" yMax="<<yMax<<" yMin="<<yMin<<endl;
if(yMax>0){
cout<<(xMax-xMin)/deltaB+1<<' '<<xMin<<' '<<yMin<<' '<<xMax<<' '<<yMax<<'\n';
}else{
cout<<xMin<<' '<<yMin<<'\n';
}
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
/*
1
3 18 6
2 1
*/
```
:::
### exgcd 与同余方程
exgcd 求解的是不定方程：
$$
ax+by=\gcd(a,b)
$$
而同余方程形如：
$$
ax\equiv c\pmod b
$$
可以将同余方程转化为：
$$
ax+by=c
$$
进而使用 exgcd 求解。
## 线性求逆元
有一种方法，可以在 $\mathcal O(\log V+n)$ 的时间内求出任意 $n$ 个数 $a_1,a_2,\cdots,a_n$ 关于 $p$ 的逆元。
记：
$$
\textit{pre}_i\equiv\prod_{j=1}^ia_j\pmod p
$$
那么就可以 $\mathcal O\left(\log V\right)$ 计算：
$$
\textit{preinv}_{n}\equiv\frac{1}{\textit{pre}_n}\pmod p
$$
$n-1\sim 1$ 递推，可以推出：
$$
\textit{preinv}_i\equiv\frac{1}{\textit{pre}_i}\equiv\frac{1}{\textit{pre}_{i+1} }\cdot a_{i+1}\pmod p
$$
随后再次递推，可以得到 $a_i$ 关于 $p$ 的逆元 $\textit{inv}_i$：
$$
\textit{inv}_i\equiv\textit{preinv}_i\cdot\textit{pre}_{i-1}\pmod p
$$
### 参考代码
```cpp
pre[0]=1;
for(int i=1;i<=n;i++){
pre[i]=1ll*pre[i-1]*a[i]%P;
}
inv[n]=qpow(pre[n],P-2);
for(int i=n-1;i>=1;i--){
inv[i]=1ll*inv[i+1]*(a[i+1])%P;
}
for(int i=1;i<=n;i++){
inv[i]=1ll*inv[i]*pre[i-1]%P;
}
```
# CRT 中国剩余定理
:::info[CRT 中的运算溢出]
**在使用 CRT 或 exCRT 时，一律建议使用 `__int128`，防止溢出**。
尽管题目大多会保证 $\operatorname{lcm}(p_1,p_2,\cdots,p_n)\leq10^{18}$，但是中间计算时会溢出。
:::
CRT（中国剩余定理）被用于求解线性同余方程组：
$$
\begin{cases}
x\equiv a_1\pmod{p_1}\\
x\equiv a_2\pmod{p_2}\\
\cdots\\
x\equiv a_n\pmod{p_n}
\end{cases}
$$
其中，$p_1,p_2,p_3,\cdots,p_n$ **两两互质**。
***
令 $L=\prod\limits_{i=1}^np_i$，则对于 $k\in\N$ 有 $k\cdot L\equiv0\pmod{p_i}$。
记 $q_i=\dfrac{L}{p_i}$，设 $c_i\equiv q_i^{-1}\pmod{p_i}$，即 $q_i$ 模 $p_i$ 意义下的逆元。
则，方程组的**最小**整数解为：
$$
x=\left(\sum_{i=1}^na_iq_ic_i\right)\bmod L
$$
同时，对于 $\forall k\in\N^*$，$x+kL$ 均为原方程组的解。
***
CRT 其实就是一个**构造式**的做法，易证 $\forall i\neq j,a_iq_ic_i\equiv0\pmod{p_j}$。
## 例题：[中国剩余定理（CRT）/ 曹冲养猪](https://www.luogu.com.cn/problem/P1495)
题目中的 $a_i,b_i$ 即分别为 $p_i,a_i$。
直接套 CRT 即可，但是需要注意的是，计算 $a_iq_ic_i$ 时需要 `__int128`，会爆 `long long`。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
constexpr const int N=10;
int n,a[N+1],p[N+1];
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
ll L=1;
for(int i=1;i<=n;i++){
cin>>p[i]>>a[i];
L*=p[i];
}
ll x=0;
for(int i=1;i<=n;i++){
ll q=L/p[i],c=inverse(q,p[i]);
x=(x+(__int128)1*a[i]%L*q%L*c)%L;
}
if(x<0){
x+=L;
}
cout<<x<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
## 扩展 CRT/exCRT
> [例题链接](https://www.luogu.com.cn/problem/P4777)
exCRT 实际上与 CRT 关系不大，如同 exLucas 与 Lucas 一般。
exCRT 解决的是当模数 $p_1,p_2,p_3,\cdots,p_n$ **不一定两两互质**的时候（这时有可能无解）。
先考虑两个**同余方程**：
$$
\begin{cases}
x\equiv a_1\pmod{p_1}\\
x\equiv a_2\pmod{p_2}\\
\end{cases}
$$
设 $r,s$，将其转化为两个**不定方程**：
$$
\begin{cases}
x=a_1+r\cdot p_1\\
x=a_2+s\cdot p_2\\
\end{cases}
$$
因此可得 $r\cdot p_1-s\cdot p_2=a_2-a_1$。
由裴蜀定理，若 $(a_2-a_1)\nmid\gcd(p_1,-p_2)$，则该不定方程**无解**，则**原同余方程组无解**。
否则通过 exgcd 求出 $r,s$，则：
$$
x=a_1+r\cdot p_1=a_2+s\cdot p_2
$$
既然 $x=a_1+r\cdot p_1$，因此有：
$$
x\equiv a_1+r\cdot p_1\pmod{\operatorname{lcm}(p_1,p_2)}
$$
这样取 $\operatorname{lcm}(p_1,p_2)$ 是为了将两个方程合并。
这样，两两顺次合并，**最终合并为一个方程**即可解。
### 参考代码
```cpp
for(int i=2;i<=n;i++){
ll w=(a[i]-a[i-1])/gcd(p[i-1],-p[i]);
ll r,s;
exgcd(p[i-1],r,-p[i],s);
r*=w,s*=w;
a[i]=(p[i-1]*r+a[i-1])%lcm(p[i-1],p[i]);
p[i]=lcm(p[i-1],p[i]);
if(i==n){
if(a[n]<0){
a[n]+=p[n];
}
cout<<a[n]<<'\n';
}
}
```
:::success[参考代码]
并没有判断无解。
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
namespace IO{
inline char getchar(){
static int p1,p2;
static char buf[1<<20];
if(p1==p2)p2=fread(buf,1,1<<20,stdin),p1=0;
return (p1==p2?EOF:buf[p1++]);
}
template<typename T>
inline void scanf(T &x){
x=0;
register int f=1;
register char ch=IO::getchar();
for(;ch<'0'||'9'<ch;ch=IO::getchar());
for(;'0'<=ch&&ch<='9';ch=IO::getchar())x=(x<<3)+(x<<1)+ch-'0';
x*=f;
}
static int p;
static char pbuf[1<<20];
inline void putchar(char ch){
pbuf[p++]=ch;
if(p==(1<<20)){
fwrite(pbuf,1,1<<20,stdout);
p=0;
}
}
template<typename T>
inline void printf(T x){
static char s[101];
int top=0;
do{
s[++top]=x%10+'0';
x/=10;
}while(x);
while(top)IO::putchar(s[top--]);
}
inline void end(){
fwrite(pbuf,1,p,stdout);
}
struct Tool{
~Tool(){
end();
}
}tool;
}
typedef __int128 lll;
#define int lll
#define ll lll
//typedef long long ll;
constexpr const int N=1e5;
int n;
ll a[N+1],p[N+1];
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll gcd(ll a,ll b){
while(b){
ll tmp=a;
a=b;
b=tmp%b;
}
return a;
}
ll lcm(ll a,ll b){
return a/gcd(a,b)*b;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
/*ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);*/
IO::scanf(n);
for(int i=1;i<=n;i++){
IO::scanf(p[i]);
IO::scanf(a[i]);
}
//合并(i-1,i)
for(int i=2;i<=n;i++){
ll w=(a[i]-a[i-1])/gcd(p[i-1],-p[i]);
ll r,s;
exgcd(p[i-1],r,-p[i],s);
r*=w,s*=w;
a[i]=(p[i-1]*r+a[i-1])%lcm(p[i-1],p[i]);
p[i]=lcm(p[i-1],p[i]);
if(i==n){
if(a[n]<0){
a[n]+=p[n];
}
IO::printf(a[n]);
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
### exexCRT
即形如：
$$
\begin{cases}
f_1x&\equiv a_1\pmod{p_1}\\
f_2x&\equiv a_2\pmod{p_2}\\
f_3x&\equiv a_3\pmod{p_3}\\
&\cdots\\
f_nx&\equiv a_n\pmod{p_n}\\
\end{cases}
$$
多了一个**系数** $f_i$。
同样可以用类似 exCRT 的方法两两合并解决，[具体见此](https://www.cnblogs.com/TH911/p/-/P4774)。
# 离散对数
所谓离散对数，即模意义下取对数。
形如：给定模数 $p$，及**整数** $a,b$，求**整数** $x$ 使得 $a^x\equiv b\pmod p$。
**注意：离散对数可能不存在**。
## BSGS 算法
在 OI 中，BSGS 算法（Baby-Step Giant-Step，大步小步算法~~北上广深算法~~）常用来求解离散对数。
BSGS 算法要求 $a\perp p$，求 $x$ 使得 $a^x\equiv b\pmod p$。
若有解，则存在 $x\leq\varphi(p)$ 的解；因为欧拉定理说明，$a\perp p$ 时，$a^{\varphi(p)}\equiv1\pmod p$。
若枚举 $x$ 是 $\mathcal O(\varphi(p))$ 的，当 $p$ 为质数的时候**不能接受**，因此可以考虑分块。
令 $B=\left\lceil\sqrt{\varphi(p)}\right\rceil$，$x=qB+r$，且 $0\leq q,r\leq B$。
则有：
$$
a^{qB+r}\equiv b\pmod p
$$
$a\perp p$，则 $a$ 的乘法逆元 $a^{-1}$ 一定存在，有：
$$
a^{qB}\equiv b\left(a^{-1}\right)^r\pmod p
$$
枚举 $r$，将其与其对应的 $b\left(a^{-1}\right)^r$ 一同存储在数据结构（一般是哈希表）中，随后枚举 $q$，在数据结构中找 $a^{qB}$ 对应的 $r$，找到了 $r$ 便找到了一个解 $x=qB+r$。同时，因为 $a\perp p$，因此 $a^{-1},a^{-2},a^{-3},\cdots,a^{-B}$ 模 $p$ 意义下**互不相同**。
时间复杂度：$\mathcal O\left(\sqrt{\varphi(p)}\right)$，在 $p$ 为质数时最劣，$\mathcal O\left(\sqrt p\right)$。
### 参考代码
```cpp
B=ceil(sqrt(euler(P)));
c=qpow(a,B);
for(int r=0,powR=1,R=qpow(a,P-2);r<B;r++){
m[1ll*b*powR%P]=r;
powR=1ll*powR*R%P;
}
for(int q=0,powC=1;q<=B;q++){
if(m.count(powC)){
cout<<q*B+m[powC]<<endl;
return 0;
}
powC=1ll*powC*c%P;
}
cout<<"no solution\n";
```
> [例题链接](https://www.luogu.com.cn/problem/P3846)
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
using namespace std;
int a,b,P,B,c;
unordered_map<int,int>m;
int euler(int n){
int ans=n;
for(int i=2;1ll*i*i<=n;i++){
if(n%i==0){
ans=ans/i*(i-1);
while(n%i==0){
n/=i;
}
}
}
if(n>1){
ans=ans/n*(n-1);
}
return ans;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>P>>a>>b;
B=ceil(sqrt(euler(P)));
c=qpow(a,B);
for(int r=0,powR=1,R=qpow(a,P-2);r<B;r++){
m[1ll*b*powR%P]=r;
powR=1ll*powR*R%P;
}
for(int q=0,powC=1;q<=B;q++){
if(m.count(powC)){
cout<<q*B+m[powC]<<endl;
return 0;
}
powC=1ll*powC*c%P;
}
cout<<"no solution\n";
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
## exBSGS 扩展 BSGS 算法
exBSGS 可解决 $a\perp p$ 不成立的情况。
由扩展欧拉定理可知，若 $x$ 有解，则 $x\leq2\varphi(p)$ 范围内也有解。
首先特判掉 $x=0$ 的情况，即 $b=1$ 或 $p=1$。
考虑分块，令 $B=\left\lceil\sqrt{2\varphi(p)}\right\rceil,x=qB-r,0\leq q,r\leq B$，则有：
$$
a^{qB-r}\equiv b\pmod p\\
\Downarrow\\
a^{qB}\equiv b\cdot a^r\pmod p
$$
因此可以预处理 $a^{qB}\bmod p$，随后枚举 $r$；但是注意到前者是后者的**充分条件而不是充要条件**，因此找出来的 $x=qB+r$ 只是“**可能的解**”，**需要检验**。
对于多个不同的 $q$ 产生的模 $p$ 意义下相同的 $a^{qB}$，可以只存储 **$q$ 最小的**两个。
:::info[证明]
由扩展欧拉定理，$a^x$ 是在 $k\leq\varphi(p)$ 值后以 $\varphi(p)$ 为周期循环的。
循环周期内的 $a^{qB}$ 显然至多存储 $1$ 个，而非循环部分 $a^{qB}$ 也至多存储 $1$ 个。
即：保留最小的 $2$ 个。
:::
时间复杂度：仍然为 $\mathcal O\left(\sqrt{\varphi(p)}\right)$。
### 参考代码
```cpp
a%=P;b%=P;
if(b==1||P==1){
cout<<"0\n";
continue;
}
B=ceil(2*sqrt(euler(P)));
c=qpow(a,B);
for(int q=0,powC=1;q<B;q++){
auto &pl=m[powC];
if(pl.size()<2){
pl.push_back(q);
}
powC=1ll*c*powC%P;
}
int x=2147483647;
for(int r=0,powA=1;r<B;r++){
auto &pl=m[1ll*b*powA%P];
for(int q:pl){
int x0=(1ll*q*B-r)%P;
if(x0<0){
continue;
}
if(qpow(a,x0)==b){
x=min(x,x0);
}
}
powA=1ll*powA*a%P;
}
if(x<2147483647){
cout<<x<<'\n';
}else{
cout<<"No Solution\n";
}
```
> [例题链接](https://www.luogu.com.cn/problem/P4195)
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
using namespace std;
int a,b,P,B,c;
//pb_ds 哈希表常数较 unordered_map 更小
__gnu_pbds::gp_hash_table<int,vector<int> >m;
//unordered_map<int,vector<int> >m;
//map<int,vector<int> >m;
int euler(int n){
int ans=n;
for(int i=2;1ll*i*i<=n;i++){
if(n%i==0){
ans=ans/i*(i-1);
while(n%i==0){
n/=i;
}
}
}
if(n>1){
ans=ans/n*(n-1);
}
return ans;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
while(true){
m.clear();
cin>>a>>P>>b;
if(a==P&&P==b&&b==0){
break;
}
a%=P;b%=P;
if(b==1||P==1){
cout<<"0\n";
continue;
}
B=ceil(2*sqrt(euler(P)));
c=qpow(a,B);
for(int q=0,powC=1;q<B;q++){
auto &pl=m[powC];
if(pl.size()<2){
pl.push_back(q);
}
powC=1ll*c*powC%P;
}
int x=2147483647;
for(int r=0,powA=1;r<B;r++){
auto &pl=m[1ll*b*powA%P];
for(int q:pl){
int x0=(1ll*q*B-r)%P;
if(x0<0){
continue;
}
if(qpow(a,x0)==b){
x=min(x,x0);
}
}
powA=1ll*powA*a%P;
}
if(x<2147483647){
cout<<x<<'\n';
}else{
cout<<"No Solution\n";
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
# Lucas 定理
## Lucas 定理
Lucas 定理用于求解大组合数取**质数模**。（对于模数不为质数的情况，请参考 exLucas。
Lucas 定理的内容很简单：
$$
\binom{n}{m}\equiv\binom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\binom{n\bmod p}{m\bmod p}\pmod p
$$
考虑如何证明。
### 证明
由二项式定理：
$$
(a+b)^p=\sum_{i=0}^p\binom pia^ib^{p-i}
$$
考虑 $\dbinom pi$ 在模 $p$ 意义下的取值。
因为：
$$
\dbinom pi=\dfrac{p!}{i!(p-i)!}
$$
那么如果化简之后，分子 $p!$ 中的 $p$ 项**没有**被约分掉，则有 $\dbinom pi\equiv p\cdot\dfrac{(p-1)!}{i!(p-i)!}\equiv0\pmod p$。
因为 $p$ 为质数，所以 $p$ 项能被约分掉当且仅当 $i!$ 中含有 $p$ 项或 $(p-i)!$ 中含有 $p$ 项。
即：$\dbinom pi\not\equiv0\pmod p$ 当且仅当 $i\equiv0\pmod p$ 或 $p-i\equiv 0\pmod p$。
在二项式定理中，$i$ 满足 $0\leq i\leq p$，所以 $\dbinom pi\not\equiv0$ 当且仅当 $i=0$ 或 $i=p$。
在这两种情况中，都可以计算得到 $\dbinom pi=1$，即 $\dbinom pi\equiv[i=0\lor i=p]$。
重新带回二项式定理，可得：
$$
\begin{aligned}
(a+b)^p&\equiv \dbinom{p}{0}a^0b^p+\dbinom ppa^pb^0\\
&\equiv a^p+b^p
\end{aligned}
\pmod p
$$
***
考虑一个二项式 $(1+x)^n\bmod p$。
$$
\begin{aligned}
(1+x)^n&\equiv (1+x)^{p\lfloor\frac np\rfloor+n\bmod p}\\
&\equiv(1+x)^{p\lfloor\frac np\rfloor}(1+x)^{n\bmod p}\\
&\equiv(1+x^p)^{\lfloor\frac np\rfloor}(1+x)^{n\bmod p}\\
\end{aligned}
\pmod p
$$
由二项式定理，$(1+x)^n$ 的 $x^m$ 项系数为 $\dbinom nm$。
而想要从 $(1+x^p)^{\lfloor\frac np\rfloor}(1+x)^{n\bmod p}$ 中得到 $x^m$，即从 $(1+x^p)^{\lfloor\frac np\rfloor}$ 中选取 $\lfloor\frac mp\rfloor$ 个 $x^p$，再从 $(1+x)^{n\bmod p}$ 中选取 $m\bmod p$ 个 $x$。
即：
$$
\dbinom nm\equiv\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}
\pmod p
$$
***
你可能有一个问题：这看起来明明应当是一个**等式**，但为什么是**同余**呢？
即，Lucas 定理应当表述为：
$$
\dbinom nm=\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}
$$
但是，你要知道，只有在模 $p$ 意义下才有 $(1+x)^{p\lfloor\frac np\rfloor}=(1+x^p)^{\lfloor\frac np\rfloor}$。
因此，只有在模 $p$ 意义下才有：
$$
\dbinom nm=\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}
$$
即：
$$
\dbinom nm\equiv\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}\pmod p
$$
### 应用
当 $n$ 比较大而无法使用其他方法（例如预处理 $1\sim n$ 的阶乘再利用乘法逆元）直接求解组合数时，可以使用 Lucas 定理。
Lucas 定理只需要递归使用即可，即递归计算 $\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}$，递归终点即 $\dbinom{0}{0}$。
其时间复杂度为：$\mathcal O\left(f(p)\log m\right)$。
其中，$f(p)$ 表示单次计算 $\dbinom{n\bmod p}{m\bmod p}$ 的复杂度，因写法而异。
可以使用乘法逆元，则时间复杂度为 $\mathcal O(\log p\log m)$。
也可以 $\mathcal O(p\log p)$ 递推，时间复杂度为 $\mathcal O(p\log p\log m)$。
推荐使用乘法逆元。
> [例题链接](https://www.luogu.com.cn/problem/P3807)
很简单，注意是 $\dbinom{n+m}{n}$ 即可。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
const int N=1e5;
int ksm(int a,int n,int p){
if(n==0)return 1;
int t=ksm(a,n>>1,p);
t=1ll*t*t%p;
if(n&1)t=1ll*t*a%p;
return t;
}
int C(int n,int m,int p){
static int ans[N+1]={1};
for(int i=1;i<=m;i++){
ans[i]=1ll*ans[i-1]*(n-i+1)%p*ksm(i,p-2,p)%p;
}return ans[m];
}
int Lucas(int n,int m,int p){
if(m==0)return 1;
return (1ll*Lucas(n/p,m/p,p)*C(n%p,m%p,p))%p;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
int T,n,m,p;
scanf("%d",&T);
while(T--){
scanf("%d %d %d",&n,&m,&p);
printf("%d\n",Lucas(n+m,n,p));
}
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
## exLucas 算法（exLucas 定理）
exLucas 算法可以求 $\dbinom{n}{m}\bmod P$，$P$ **不一定为质数**。（但是，exLucas 实际上和 Lucas 没有多大关系……）
由唯一分解定理，可以对 $P$ 进行质因数分解：
$$
P=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_n^{c^n}
$$
### CRT 求解
我们如果能够求出 $a_1,a_2,\cdots,a_n$ 使得：
$$
\begin{cases}
\dbinom nm&\equiv a_1\pmod{p_1^{c_1}}\\
\dbinom nm&\equiv a_2\pmod{p_1^{c_2}}\\
&\cdots\\
\dbinom nm&\equiv a_n\pmod{p_n^{c_n}}\\
\end{cases}
$$
那么我们就可以通过 **CRT** 求解出 $\dbinom nm\bmod P$。因为在 CRT 中，恰好也有 $P=p_1^{c_1}p_2^{c_2}\cdots p_n^{c_n}$。
### 求解模质数幂下余数
展开定义式：
$$
\dbinom nm=\frac{n!}{m!(n-m)!}
$$
因此：
$$
\dbinom nm\equiv\frac{n!}{m!(n-m)!}\pmod{p_i^{c_i}}
$$
因为 $p_i^{c_i},m!,(n-m)!$ 不一定互质，因此**乘法逆元不一定存在**。
考虑先**提出**分子和分母中所有的 $p_i$ 次幂，随后便可以使用逆元求解。
记 $x$ 的质因数分解中**质数** $p$ 的幂次为 $\nu_p(x)$，剩余积为 $(x)_p$，即：
$$
x=p^{\nu_p(x)}\cdot(x)_p
$$
则有：
$$
\frac{n!}{m!(n-m)!}\equiv\frac{(n!)_{p_i}}{(m!)_{p_i}((n-m)!)_{p_i}}\cdot p_i^{\nu_{p_i}(n!)-\nu_{p_i}(m!)-\nu_{p_i}((n-m)!)}\pmod{p_i^{c_i}}
$$
那么，现在只需要考虑对于 $x,p$，如何**高效地**求出 $\nu_{p}(x!)\bmod p_i^{c_i},(x!)_p\bmod p_i^{c_i}$。
***
考虑：
$$
x!=1\times2\times\cdots\times p\times\cdots\times2p\times\cdots\times\left\lfloor\dfrac{x}{p}\right\rfloor p\times\cdots\times x
$$
容易发现，在 $p,2p,3p,\cdots,\left\lfloor\dfrac xp \right\rfloor p$ 中 $p$ 的个数为 $\left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)$。其中 $\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)$ 是 $1,2,3,\cdots,\left\lfloor\dfrac xp\right\rfloor$ 中的 $p$ 的个数。
$p$ 为**质数**，所以在 $1,2,3,\cdots,p-1,p+1,\cdots$ 中不会含有 $p$。
将递推式展开：
$$
\begin{aligned}
\nu_p\left(x!\right)&=\left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)\\
&=\left\lfloor\dfrac xp\right\rfloor+\left\lfloor\dfrac x{p^2}\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp^2 \right\rfloor!\right)\\
&=\left\lfloor\dfrac xp\right\rfloor+\left\lfloor\dfrac x{p^2}\right\rfloor+\left\lfloor\dfrac x{p^3}\right\rfloor+\cdots
\end{aligned}
$$
因此可以 $\mathcal O(\log_p x)$ 计算 $\nu_p(x!)$：
```cpp
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;//这里取不取模其实都无所谓,因为幂次不会太多,想取也可以.
}
```
***
现在考虑如何计算 $(x!)_p\bmod p^c$；显然不能利用定义式 $(x!)_p=\dfrac{x!}{\nu_p(x!)}$，而需要其他方法（因为无法得知 $x!$）。
不难进行递推：
$$
\begin{aligned}
(x!)_p&\equiv\prod_{i=1}^n(i)_p\\
&\equiv\left(\prod_{1\leq i\leq x,i\perp p}(i)_p\right)\left(\prod_{i=1}^{\left\lfloor\frac xp\right\rfloor}(p\cdot i)_p\right)\\
&\equiv\left(\prod_{1\leq i\leq x,i\perp p}(i)_p\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p\\
&\equiv\left(\prod_{1\leq i\leq p^c,i\perp p}i\right)^{\left\lfloor\frac{x}{p^c}\right\rfloor}\left(\prod_{1\leq i\leq x\bmod p^c,i\perp p}i\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p\\
\end{aligned}
$$
由 Wilson 定理的推论（$m=2,4,p^c,2p^c$ 时 $k\equiv-1$，否则为 $1$）：
$$
\prod_{1\leq k\leq m,k\perp m}k\equiv\pm1\pmod m
$$
因此：
$$
\begin{aligned}
(x!)_p&\equiv(\pm1)^{\left\lfloor\frac{x}{p^c}\right\rfloor}\left(\prod_{1\leq i\leq x\bmod p^c,i\perp p}i\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p
\end{aligned}
$$
可以发现，每次计算 $(x!)_p$ 时，$p^c$ 是固定的，因此可以预处理：
$$
f_j\equiv\prod_{1\leq i\leq j,i\perp p}i\pmod p^c
$$
因此，得到最终推导式：
$$
\begin{aligned}
(x!)_p&\equiv(\pm1)^{\left\lfloor\frac{x}{p^c}\right\rfloor}f_{x\bmod p^c}\left(\left\lfloor\frac xp\right\rfloor!\right)_p
\end{aligned}
$$
可以递归或迭代处理，时间复杂度 $\mathcal O(p^c+\log n x)$。
```cpp
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc)
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
```
> [例题链接](https://www.luogu.com.cn/problem/P4720)
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
#define ll __int128
#define int __int128
constexpr const int PMAX=1e6;
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base;
}
base=1ll*base*base;
n>>=1;
}
return ans;
}
int CRT(vector<int>a,vector<int>p){
ll L=1;
for(int &i:p){
L*=i;
}
ll x=0;
for(int i=0;i<a.size();i++){
ll q=L/p[i];
x=(x+1ll*a[i]*q%L*inverse(q,p[i])%L)%L;
}
if(x<0){
x+=L;
}
return x;
}
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;
}
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc);
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
void breakDown(int P,vector<int>&p,vector<int>&pc){
int pp=P;
for(int i=2;i*i<=pp;i++){
if(pp%i==0){
p.push_back(i);
pc.push_back(1);
while(pp%i==0){
pc.back()*=i;
pp/=i;
}
}
}
if(pp>1){
p.push_back(pp);
pc.push_back(pp);
}
}
long long exLucas(ll n,ll m,ll P){
if(n<m){
return 0;
}
vector<int>p,pc;
breakDown(P,p,pc);
vector<int>a(pc.size());
for(int i=0;i<p.size();i++){
int nWilson=Wilson(n,p[i],pc[i]);
int mWilson=Wilson(m,p[i],pc[i]);
int nmWilson=Wilson(n-m,p[i],pc[i]);
a[i]=nWilson*inverse(mWilson*nmWilson%pc[i],pc[i])*qpow(p[i],v(n,p[i])-v(m,p[i])-v(n-m,p[i]));
}
return CRT(a,pc);
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
long long n,m,P;
cin>>n>>m>>P;
cout<<exLucas(n,m,P)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
/*
1000000000000000000 500000000000000000 998243
0
*/
```
:::
> [[国家集训队] 礼物](https://www.luogu.com.cn/problem/P2183)
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
#define ll __int128
#define int __int128
constexpr const int M=5;
long long w[M+1];
constexpr const int PMAX=1e5;
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base;
}
base=1ll*base*base;
n>>=1;
}
return ans;
}
int CRT(vector<int>a,vector<int>p){
ll L=1;
for(int &i:p){
L*=i;
}
ll x=0;
for(int i=0;i<a.size();i++){
ll q=L/p[i];
x=(x+1ll*a[i]*q%L*inverse(q,p[i])%L)%L;
}
if(x<0){
x+=L;
}
return x;
}
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;
}
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc);
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
void breakDown(int P,vector<int>&p,vector<int>&pc){
int pp=P;
for(int i=2;i*i<=pp;i++){
if(pp%i==0){
p.push_back(i);
pc.push_back(1);
while(pp%i==0){
pc.back()*=i;
pp/=i;
}
}
}
if(pp>1){
p.push_back(pp);
pc.push_back(pp);
}
}
long long exLucas(ll n,ll m,ll P){
if(n<m){
return 0;
}
vector<int>p,pc;
breakDown(P,p,pc);
vector<int>a(pc.size());
for(int i=0;i<p.size();i++){
int nWilson=Wilson(n,p[i],pc[i]);
int mWilson=Wilson(m,p[i],pc[i]);
int nmWilson=Wilson(n-m,p[i],pc[i]);
a[i]=nWilson*inverse(mWilson*nmWilson%pc[i],pc[i])*qpow(p[i],v(n,p[i])-v(m,p[i])-v(n-m,p[i]));
}
return CRT(a,pc);
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
long long n,m,P;
cin>>P>>n>>m;
long long ans=1;
for(int i=1;i<=m;i++){
cin>>w[i];
ans=ans*exLucas(n,w[i],P)%P;
n-=w[i];
if(n<0){
cout<<"Impossible"<<endl;
return 0;
}
}
cout<<ans<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
# 数论分块
数论分块（整除分块）用于快速计算：
$$
\sum_{i=1}^nf(i)g\left(\left\lfloor\dfrac ni\right\rfloor\right)
$$
数论分块的原理即本质不同的 $\left\lfloor\dfrac ni\right\rfloor$ 只有至多 $\mathcal O\left(\sqrt n\right)$ 种。
:::info[证明]
* 当 $i\leq\sqrt n$ 时，$i$ 只有 $\mathcal O\left(\sqrt n\right)$ 种取值。
* 当 $i>\sqrt n$ 时，$\dfrac ni\leq\sqrt n$，$\left\lfloor\dfrac ni\right\rfloor$ 只有 $\mathcal O\left(\sqrt n\right)$ 种取值。
故，本质不同的 $\left\lfloor\dfrac ni\right\rfloor$ 只有至多 $\mathcal O\left(\sqrt n\right)$ 种取值。
:::
因此 $\left\lfloor\dfrac ni\right\rfloor$ 相同值肯定是连续的，那么就可以找出这 $\mathcal O\left(\sqrt n\right)$ 区间 $[l_1,r_1],[l_2,r_2],\cdots,[l_k,r_k]$（$r_i+1=l_{i+1}$），求出：
$$
g\left(\left\lfloor\dfrac nl\right\rfloor\right)\sum_{j=l_i}^{r_i}f(j)
$$
将 $i$ 个答案相加即可。
若分块后可以 $\mathcal O(1)$ 求出，则数论分块是 $\mathcal O\left(\sqrt n\right)$ 的。
## 寻找区间
对于一个 $l$ 所对应的 $\left\lfloor\dfrac nl\right\rfloor$，要寻找一个最大的 $r$ 使得 $\left\lfloor\dfrac nl\right\rfloor=\left\lfloor\dfrac nr\right\rfloor$，从而找到区间 $[l,r]$。
则有：
$$
r=\left\lfloor\dfrac{n}{\left\lfloor\dfrac nl\right\rfloor}\right\rfloor
$$
而对于向上取整的情形，即对于 $l$ 对应的 $\left\lceil\dfrac nl\right\rceil$，要寻找一个最大的 $r$ 使得 $\left\lceil\dfrac nr\right\rceil$，从而找到区间 $[l,r]$。
则有：
$$
r=\left\lfloor\dfrac{n-1}{\left\lfloor\dfrac {n-1}l\right\rfloor}\right\rfloor
$$
## 单一参数
> [[CQOI2007] 余数求和](https://www.luogu.com.cn/problem/P2261)
>
> 给定 $n,k\in\N^*$，满足 $1\leq n,k\leq10^9$，求：
>
> $$
> G(n,k)=\sum_{i=1}^nk\bmod i
> $$
显然，$k\bmod i=k-i\left\lfloor\dfrac ki\right\rfloor$。
因此有：
$$
\begin{aligned}
G(n,k)&=\sum_{i=1}^n\left(k-i\left\lfloor\dfrac ki\right\rfloor\right)\\
&=nk-\sum_{i=1}^ni\left\lfloor\dfrac ki\right\rfloor
\end{aligned}
$$
考虑快速求出 $\sum\limits_{i=1}^ni\left\lfloor\dfrac ki\right\rfloor$，可以进行**数论分块**。
对于数论分块中一组确定的 $l,r$，有：
$$
\sum_{i=l}^ri\left\lfloor\dfrac ki\right\rfloor=\left\lfloor\dfrac kl\right\rfloor\sum_{i=l}^ri=\left\lfloor\dfrac ki\right\rfloor\dfrac{(l+r)(r-l+1)}{2}
$$
这可以 $\mathcal O(1)$ 计算，因此总计算时间复杂度为 $\mathcal O\left(\sqrt n\right)$。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
#define int ll
typedef long long ll;
ll G(int n,int k){
ll ans=1ll*n*k;
for(int l=1,r;l<=n;l=r+1){
int t=k/l;
if(!t){
r=n;
}else{
r=min(k/t,n);
}
ans-=1ll*(r-l+1)*t*(l+r)>>1;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,k;
cin>>n>>k;
cout<<G(n,k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
## 多参数
此时寻找 $r$，会取 $\min$，即：
$$
r=\min\left(\left\lfloor\dfrac{n}{\left\lfloor\dfrac nl\right\rfloor}\right\rfloor,\left\lfloor\dfrac{m}{\left\lfloor\dfrac ml\right\rfloor}\right\rfloor,\cdots\right)
$$
> [[清华集训 2012] 模积和](https://www.luogu.com.cn/problem/P2260)
>
> 给定 $1\leq n,m\leq10^9$，求：
>
> $$
> \left(\sum_{i=1}^n\sum_{j=1}^m[i\neq j](n\bmod i)(m\bmod j)\right)\bmod 19940417
> $$
不妨令 $n\leq m$，否则交换 $n,m$。
模运算不好处理，转换一下：
$$
\sum_{i=1}^n\sum_{j=1}^m[i\neq j]\left(n-i\left\lfloor\dfrac ni\right\rfloor\right)\left(m-j\left\lfloor\dfrac mj\right\rfloor\right)
$$
$[i\neq j]$ 不好处理，因此可以先全部求出来，再减去相等的部分：
$$
\sum_{i=1}^n\sum_{j=1}^m\left(n-i\left\lfloor\dfrac ni\right\rfloor\right)\left(m-j\left\lfloor\dfrac mj\right\rfloor\right)-\sum_{i=1}^n\left(n-i\left\lfloor\dfrac ni\right\rfloor\right)\left(m-i\left\lfloor\dfrac mi\right\rfloor\right)
$$
不难发现 $j$ 与 $i$ 的取值**无关**，因此转换为：
$$
\sum_{i=1}^n\left(n-i\left\lfloor\dfrac ni\right\rfloor\right)\sum_{j=1}^m\left(m-j\left\lfloor\dfrac mj\right\rfloor\right)-\sum_{i=1}^n\left(n-i\left\lfloor\dfrac ni\right\rfloor\right)\left(m-i\left\lfloor\dfrac mi\right\rfloor\right)
$$
于是可以 $\mathcal O\left(\max(n,m)\right)$ 计算，但是考虑到 $1\leq n,m\leq10^9$，只需要使用数论分块优化即可。
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int P=19940417,inv2=9970209,inv6=3323403;
#define ll __int128
//typedef long long ll;
ll g(ll n,ll m){
ll ans=0;
for(ll l=1,r=n;l<=n;l=r+1,r=n){
ll t=m/l;
if(t){
r=min(n,m/t);
}
ans=(ans+t*(l+r)%P*(r-l+1)%P*inv2%P)%P;
}
return ans;
}
ll g2(ll n,ll m){
ll ans=0;
for(ll l=1,r=n;l<=n;l=r+1,r=n){
ll tn=n/l,tm=m/l;
r=min(n/tn,m/tm);
ans+=m*(l+r)%P*(r-l+1)*inv2*tn%P;
ans+=n*(l+r)%P*(r-l+1)*inv2*tm%P;
ans-=(r*(r+1)*(2*r+1)%P-(l-1)*l*(2*l-1)%P)*inv6*tn*tm%P;
ans%=P;
}
return ans;
}
int f(ll n,ll m){
ll ans=(n*n%P-g(n,n))*(m*m%P-g(m,m))%P;
ans=(ans-n*n%P*m%P)%P;
ans=(ans+g2(n,m))%P;
if(ans<0){
ans+=P;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,m;
cin>>n>>m;
if(n>m){
swap(n,m);
}
cout<<f(n,m)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::
# 积性函数
## 定义
若数论函数 $f(n)$ 满足 $f(1)=1$，且对于 $\forall x,y\in\N^*\land x\perp y$，有 $f(xy)=f(x)\cdot f(y)$，称 $f(n)$ 为**积性函数**。
若数论函数 $f(n)$ 满足 $f(1)=1$，且对于 $\forall x,y\in\N^*$，有 $f(xy)=f(x)\cdot f(y)$，称 $f(n)$ 为**完全积性函数**。
## 基础部分
若 $f(n),g(n)$ 均为积性函数，则 $f(x^p),f^p(x),f(x)g(x),(f*g)(x)$ 均为积性函数。其中 $f*g$ 为迪利克雷卷积：
$$
(f*g)(x)=\sum_{d\mid x}f(d)g\left(\dfrac xd\right)
$$
### 常见积性函数
* **单位函数**：$\varepsilon(n)=[n=1]$，完全积性函数。（`\varepsilon`）
* **恒等函数**：$\mathrm{id}_{k}(n)=n^k$，$\mathrm{id}_1(n)$ 通常记为 $\mathrm{id}(n)$，完全积性函数。
* 常数函数：$1(n)=1$，完全积性函数。
* 除数函数：$\sigma_k(n)=\sum\limits_{d\mid n}d^k$，$\sigma_0(n)$ 通常记作 $d(n)$ 或 $\tau(n)$，$\sigma_1(n)$ 通常记作 $\sigma(n)$。（`\sigma`）
:::info[$d(n)$ 的求法]
令 $n=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}$，其中 $p_1,p_2,\cdots,p_k$ 均为质数。
对于每一个指数 $c_i$ 的更改，都会产生新的因数。每一个 $i$，都有 $c_i+1$ 种选择。显然有：
$$
d(n)=\prod_{i=1}^k(c_i+1)
$$
时间复杂度 $\mathcal O\left(\sqrt n\right)$。
:::
* 欧拉函数：$\varphi(n)$。（`\varphi`）
* 莫比乌斯函数：
$$
\mu(n)=
\begin{cases}
1&n=1\\
0&\exist d>1,d^2\mid n\\
(-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}
\end{cases}
$$
$\mu(n)=0$ 即 $n$ 含有平方因子。（`\mu`）
同时，除数函数 $d(n)$ 即 $n$ 的因数个数，且 $d(ij)$ 满足：
$$
\begin{aligned}
d(ij)&=\sum_{x\mid i}\sum_{y\mid j}[x\perp y]
\end{aligned}
$$
:::info[证明]
令 $i=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k},j={p'_1}^{c'_1}{p'_2}^{c'_2}\cdots {p'_{k'}}^{c'_{k'}}$。如果将指数对位相加，则可以得到 $ij$ 的质因数分解形式。
将每一个质因子 $p_l$ 枚举 $c_l+c'_l+1$ 次，便可以得到全部的因数。因此枚举 $x\mid i,y\mid j$，但是需要保证 $x,y$ 不包含重复质因子，即 $x\perp y$。
:::
## 迪利克雷卷积
对于数论函数 $f,g$，记 $h=f*g$ 为 $f,g$ 通过迪利克雷卷积相乘的到的结果，即：
$$
h(n)=\sum_{d\mid n}f(d)g\left(\dfrac nd\right)
$$
### 常见迪利克雷卷积
* $\varphi*1=\mathrm{id}$。
* $\mu*1=\varepsilon$。
* $\mu*1=\varphi$。
## 莫比乌斯函数
莫比乌斯函数 $\mu(n)$ 满足：
$$
\sum_{d\mid n}\mu(d)=[n=1]
$$
即 $\mu*1=\varepsilon$。
因此，$\mu(n)$ 可用于处理**互质**相关信息：
$$
[i\perp j]=[\gcd(i,j)=1]=\sum_{d\mid\gcd(i,j)}\mu(d)=\sum_{d\mid i\land d\mid j}\mu(d)
$$
### 莫比乌斯反演/莫比乌斯变换
若 $g(n)=\sum\limits_{d\mid n}f(d)$，则有：
$$
f(n)=\sum_{d\mid n}\mu(d)g\left(\dfrac nd\right)
$$
因为：
$$
\begin{aligned}
\sum_{d\mid n}\mu(d)g\left(\dfrac nd\right)&=\sum_{d\mid n}\mu(d)\sum_{d'\mid\large\frac nd}f(d')\\
&=\sum_{d\mid n}\mu(d)\sum_{d'\mid\large\frac nd}f(d')\\
&=\sum_{d\mid n}f(d)\sum_{d'\mid\large \frac nd}\mu(d')\\
&=\sum_{d\mid n}f(d)\left[n=d\right]\\
&=f(n)
\end{aligned}
$$
若 $g(n)=\sum\limits_{n|d}f(d)$，则有：
$$
f(n)=\sum_{n\mid d}\mu\left(\dfrac dn\right)g(d)
$$
> [[国家集训队] Crash的数字表格 / JZPTAB](https://www.luogu.com.cn/problem/P1829)
>
> [详见此处](https://www.cnblogs.com/TH911/p/-/P1829)。
莫比乌斯反演的本质其实就是高维差分。
## 线性筛
线性筛可以用于求解欧拉函数及莫比乌斯函数。实际上，**线性筛几乎可以求解所有的积性函数**，求解方法视具体函数而定。
### 欧拉函数
若 $i$ 为质数，显然 $\varphi(i)=i-1$。
设质数 $\textit{prime}_j$。
若 $i\bmod\textit{prime}_j\neq0$，则 $i\perp\textit{prime}_j$，因为 $\varphi$ 是一个积性函数，于是有：
$$
\varphi\left(i\cdot\textit{prime}_j\right)=\varphi(i)\cdot\varphi\left(\textit{prime}_j\right)
$$
否则，当 $i\equiv0\pmod{\textit{prime}_j}$ 时，有：
$$
\varphi\left(i\cdot\textit{prime}_j\right)=\varphi(i)\cdot\textit{prime}_j
$$
因为与 $i\cdot\textit{prime}_j$ 互质，即与 $i$ 互质。$1,2,3,\cdots,i\cdot\textit{prime}_j$ 在模 $i$ 意义下以 $i$ 为周期循环了 $\textit{prime}_j$ 次。
### 莫比乌斯函数
若 $i$ 为质数，显然 $\mu(i)=-1$。
设质数 $\textit{prime}_j$。
若 $i\bmod\textit{prime}_j\neq0$ 时，则 $i\perp\textit{prime}_j$，因为 $\mu$ 是一个积性函数，于是有：
$$
\mu\left(i\cdot\textit{prime}_j\right)=\mu(i)\cdot\mu\left(\textit{prime}_j\right)=-\mu(i)
$$
否则，当 $i\equiv0\pmod{\textit{prime}_j}$ 时，有：
$$
\mu\left(i\cdot\textit{prime}_j\right)=0
$$
因为 $i$ 含有 $\textit{prime}_j$，则 $i\cdot\textit{prime}_j$ 至少含有两个 $\textit{prime}_j$，故 $\mu\left(i\cdot\textit{prime}_j\right)=0$。
### 参考代码
```cpp
//ans1 为欧拉函数，ans2 为莫比乌斯函数
void pre(){
static int vis[N+1],prime[N+1],size;
ans1[1]=ans2[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++size]=i;
vis[i]=i;
ans1[i]=i-1;
ans2[i]=-1;
}
for(int j=1;j<=size&&i*prime[j]<=N;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0){
//ans2 为 0,可以不写
ans1[i*prime[j]]=ans1[i]*prime[j];
break;
}
ans1[i*prime[j]]=ans1[i]*ans1[prime[j]];
ans2[i*prime[j]]=-ans2[i];
}
}
//前缀和,可不写
for(int i=1;i<=N;i++){
ans1[i]+=ans1[i-1];
ans2[i]+=ans2[i-1];
}
return;
}
```
## 杜教筛
杜教筛用于快速求解积性函数 $f$ 的前缀和 $S(n)=\sum\limits_{i=1}^nf(i)$。（实际上，只要能够构造恰当的函数，均可以使用杜教筛求解。）
构造一组 $h=f*g$，有：
$$
\begin{aligned}
\sum_{i=1}^nh(i)&=\sum_{i=1}^n\sum_{d\mid i}f\left(\dfrac id\right)g(d)\\
&=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}f\left(i\right)\\
&=\sum_{i=1}^ng(i)S\left(\left\lfloor\frac ni\right\rfloor\right)
\end{aligned}
$$
可得：
$$
\begin{aligned}
g(1)S(n)&=\sum_{i=1}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right)\\
&=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right)
\end{aligned}
$$
因此，只要 $h,g$ 适当，能够快速地求出来，就能够在较短时间内求出 $S(n)$。
### 数论分块
可以发现，$\sum\limits_{i=2}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right)$ 直接计算的复杂度较高，而 $\left\lfloor\dfrac ni\right\rfloor$ 就很适合用于数论分块来加速。
### 记忆化
每当我们求出了 $S(n)$ 之后，都应当将其记录下来，避免重复计算。
这样的杜教筛的时间复杂度时 $\mathcal O\left(n^{\frac34}\right)$ 的。
### 线性筛预处理
可以预处理较小的一部分的 $S(n)$，从而节约时间。
当预处理前 $\mathcal O\left(n^\frac23\right)$ 的 $S(n)$ 时，时间复杂度最小，为 $\mathcal O\left(n^\frac23\right)$。
#### 实际常数影响
在 OI 代码中，递归求解 $S(n)$ 会有常数的影响。
因此最好的方法应当是实际测试预处理部分大小 $m$ 的值从而确定时间最小的 $m$，这样即使复杂度劣一些，使用时间也小一些。
> [例题链接](https://www.luogu.com.cn/problem/P4213)
:::success[参考代码]
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace std;
typedef long long ll;
__gnu_pbds::gp_hash_table<int,ll>mem1,mem2;
constexpr const int N=5e6/*1664510*/;
ll ans1[N+1],ans2[N+1];
ll S1(int n){
if(n<=N){
return ans1[n];
}
if(mem1[n]){
return mem1[n];
}
ll ans=n*(n+1ll)>>1;
for(ll l=2,r=n;l<=n;l=r+1,r=n){
ll t=n/l;
r=n/t;
ans-=(r-l+1ll)*S1(t);
}
return mem1[n]=ans;
}
ll S2(int n){
if(n<=N){
return ans2[n];
}
if(mem2[n]){
return mem2[n];
}
ll ans=1;
for(ll l=2,r=n;l<=n;l=r+1,r=n){
ll t=n/l;
r=n/t;
ans-=(r-l+1ll)*S2(t);
}
return mem2[n]=ans;
}
void pre(){
static int vis[N+1],prime[N+1],size;
ans1[1]=ans2[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++size]=i;
vis[i]=i;
ans1[i]=i-1;
ans2[i]=-1;
}
for(int j=1;j<=size&&i*prime[j]<=N;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0){
ans1[i*prime[j]]=ans1[i]*prime[j];
break;
}
ans1[i*prime[j]]=ans1[i]*ans1[prime[j]];
ans2[i*prime[j]]=-ans2[i];
}
}
for(int i=1;i<=N;i++){
ans1[i]+=ans1[i-1];
ans2[i]+=ans2[i-1];
}
return;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
pre();
ll T;
cin>>T;
while(T--){
ll n;
cin>>n;
mem1.clear();
mem2.clear();
cout<<S1(n)<<' '<<S2(n)<<'\n';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
```
:::