数论题单与总结(By Vector_net)
题单介绍
在 OI 界，数学一直是非常重要的一个板块，也是最让人头疼的板块之一。但实际上，只要你搞懂了数论的原理，数论就会非常简单。
本题单将会实时更新。
本题单将会介绍：(注：本题单中所有要使用龟速乘处均使用```__int128```代替。)
- 快速幂
- 筛法
- 矩阵快速幂
- Exgcd
- 乘法逆元
- Lucas定理
- ExCRT，~~因为CRT没用~~
- 扩展欧拉定理
- 原根
- BSGS算法
- exLucas定理
- exBSGS
- 二次剩余
- 数论分块
- 莫比乌斯反演
## Part-1：快速幂
我们考虑这样一个问题，求 $a^b\bmod p$，其中$1 \leq a,b \leq 10^9$，暴力的想法自然是一个个去乘，但庞大的数据范围使暴力不可能通过。我们考虑优化这个想法。
首先我们的一个想法是分治，由初中数学可知
$a^b = \begin{cases}
(a^\frac{b}{2})^2 & x \equiv 0 \pmod{2} \\
(a^\frac{b-1}{2})^2 \cdot a & x \equiv 1 \pmod{2} \\
\end{cases}$
按照这个直接分治递归即可，时间复杂度 $O(\log n)$。
但分治有时候会TLE(~~但凡数据大一点再有个多组就炸了~~)，我们可以考虑不使用递归，就是迭代版了。我们对 $b$ 进行二进制拆分。从低位(从右往左数)向高位计算，答案可以写成 $a^b=\prod_{i=1}^n a^{2^{k_i}}$ 按照这个式子递推计算即可，时间复杂度 $O(\log n)$。(更快)
### 代码如下：
分治版：
```cpp
inline int qpow(int a,int b,int p){
if(!b)return 1;
int z=qpow(a,b>>1,p);z=z*z%p;
if(b&1)return z*a%p;else return z;
}
```
迭代版：
```cpp
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
```
## Part-2 筛法
1.Eratosthenes筛法
基本思想：质数的倍数不是质数，遇到质数时枚举筛一下即可。
代码：
```cpp
inline int make_p(int n){
vis[0]=vis[1]=1;
for(int i=2;i<=n;i+=(i&1)+1)
if(!vis[i])
for(int j=i*i;j<=n;j+=i)vis[j]=1;
}
```
时间复杂度 $O(n \log \log n)$。
2.线性筛(欧拉筛)
Eratosthenes筛法即使在优化后也会重复标记合数，根本原因是算法找不到唯一产生数的方式，于是线性筛应运而生。具体而来，我们每次只让当前数乘上一个质因子，并让她是这个合数的最小质因子，具体见代码。(这个东西在分解质因数(默认大家都会 $O(\sqrt n)$ 的分解)中很有用，可以 $O(n)$ 预处理 $O(\log n)$分解)
代码：
```cpp
inline void make_p(int n){
vis[0]=vis[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i])prime[++cnt]=i,vis[i]=i;
for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0)break;
}
}
}
```
## Part-3 矩阵快速幂
矩阵快速幂其实就是矩阵乘法+快速幂，快速幂在此不加赘述，讲一讲矩阵乘法。
两个大小分别为 $m \times n$ 和 $n \times p$ 的矩阵 $A, B$ 相乘的结果为一个大小为 $m \times p$ 的矩阵。将结果矩阵记作 $C$，则
$$ c_{i j} = \sum_{k = 1}^{n} a_{i k} b_{k j} \text{,\qquad($1 \le i \le m$, $1 \le j \le p$).} $$
而如果 $A$ 的列数与 $B$ 的行数不相等，则无法进行乘法。
可以验证，矩阵乘法满足结合律，即 $(A B) C = A (B C)$。
完全可以打一个重载运算符。
矩阵乘法代码：
```cpp
struct Matrix{
int a[maxn][maxn],n,m;
inline Matrix(){memset(a,0,sizeof a);n=m=0;}
inline Matrix operator *(const Matrix B){
Matrix C;C.n=n,C.m=B.m;
for(int i=1;i<=n;i++)
for(int j=1;j<=B.m;j++)
for(int k=1;k<=m;k++)C.a[i][j]=(C.a[i][j]+a[i][k]*B.a[k][j]%mod)%mod;
return C;
}
inline void debug(){
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++)printf("%lld ",a[i][j]);
printf("\n");
}
}
}MatA,Mat_1;//Mat_1是单位矩阵，即对角线全为1，其他全为0的矩阵，性质是Mat_1*Mat_A=Mat_A
```
矩阵快速幂代码：
```cpp
inline Matrix qpow(Matrix a,int b){
Matrix res=Mat_1;
for(;b;a=a*a,b>>=1)
if(b&1)res=res*a;
return res;
}
```
矩阵快速幂最大的用处是优化递推与DP。(把递推式抽象成一个初始矩阵和一个转移矩阵，然后根据结合律写快速幂即可，斐波那契数列最典型了)时间复杂度 $O(m^3 \log n)$。
## Part-4 Exgcd
首先我们知道最大公约数(gcd)的 $O(\log V)$ 求法，即辗转相除。扩展欧几里得算法就是通过其中的 $gcd(a,b)=gcd(b,a\bmod b)$ 写出的。
Exgcd 算法主要用于求不定方程 $ax+by=c$ 的解，由裴蜀定理知，此时一定有 $gcd(a,b) \mid c$，于是我们只需要解方程 $ax+by=gcd(a,b)$，这个可以扩欧。
具体证明：
$\because ax+by=gcd(a,b)$
$\therefore ax+by=gcd(b,a \bmod b)$
令 $bx^,+(a \bmod b)y^,=gcd(b,a \bmod b)$
则有 $x=y^,$,$y=x- \lfloor \frac{a}{b} \rfloor y^,$
注意到这个问题可以由一个更小的子问题得来，于是可以递归求解，时间复杂度 $O(\log V)$。由此可以得到一组解 $x_0$ 与 $y_0$，令 $x=\frac{x_0c}{gcd(a,b)}$，$y=\frac{y_0c}{gcd(a,b)}$，$d_x=\frac{b}{gcd(a,b)}$，$d_y=\frac{a}{gcd(a,b)}$。则特解为 $(x,y)$，通解为 $(x+td_x,y-td_y) (t \in \mathbb Z)$
对于同余方程，我们可以把她写成 $ax+by=c$ 的形式，可以使用Exgcd求解。
代码：
```cpp
inline void exgcd(int a,int b,int &x,int &y){
if(!b)return x=1,y=0,void();
else exgcd(b,a%b,y,x),y-=a/b*x;
}
```
## Part-5 乘法逆元
对于模 $p$ 域下整数的加法、减法、乘法我们都有相应的解决方案，唯独除法有些难，乘法逆元便用于解决这个问题。
定义：模 $p$ 域下数 $a$ 的逆元记作 $a^{-1}$，定义为同余方程 $ax \equiv 1 \pmod p$ 的最小正整数解。当 $gcd(p,a) \ne 1$ 时不存在逆元。
求法：1.扩展欧几里得：这个直接解就可以了，而且非常快，基本可以看作小常数，时间复杂度 $O(\log p)$。(但估计复杂度时还是估计上为好，~~万一毒瘤给你造个斐波那契数列~~)
2.快速幂+欧拉函数 $\varphi$：($\varphi(p)$ 表示 $1$ 到 $p-1$ 内与 $p$ 互质的数的数量，计算公式是 $\varphi(p)=n\prod_{i=1}^m \frac{p_i-1}{p_i}$)首先，由欧拉定理可知有 $a^{\varphi(p)} \equiv 1 \pmod{p}$，于是有 $a^{-1} \equiv a^{\varphi(p)-1} \pmod{p}$，快速幂求解即可。当然可以先 $O(n)$ 预处理欧拉函数(欧拉函数是积性函数，可以使用欧拉筛 $O(n)$ 预处理)使时间复杂度变为 $O(\log p)$。
代码：
$O(n)$ 预处理欧拉函数代码：
```cpp
inline void make_Phi(int n){
vis[0]=vis[1]=1,phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i])prime[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
else phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
}
```
Exgcd：
```cpp
inline void exgcd(int a,int b,int &x,int &y){
if(!b)return x=1,y=0,void();
else exgcd(b,a%b,y,x),y-=a/b*x;
}
inline int inv(int a,int p){
int x,y;exgcd(a,p,x,y);
return (x%p+p)%p;
}
```
快速幂：
```cpp
inline int phi(int n){
int ret=n;
for(int i=2;i*i<=n;i+=(i&1)+1){
if(n%i==0){
ret=ret/i*(i-1);
while(n%i==0)n/=i;
}
}
if(n>1)ret=ret/n*(n-1);
return ret;
}
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
inline int inv(int a,int p){
int Phi=phi(p);
return qpow(a,Phi-1,p);
}
```
当然，对于 $n$ 个数，我们不需要傻傻的对每个数分别求逆元，存在一种可以做到 $O(n+\log p)$ 的优秀方法。我们可以先预处理出数列在模 $p$ 意义下的前缀积 $L[i]$ 与后缀积 $R[i]$，还有数列中所有数的积的逆元 $Q$。那么有 $a_i^{-1} \equiv L[i-1] \times R[i+1] \times Q \pmod{p}$，这样我们就得到了一个非常优秀的做法。
## Part-6 Lucas 定理
Lucas 定理其实就是一个公式：
$C_n^m \equiv C_{\lfloor \frac{n}{p} \rfloor}^{\lfloor \frac{m}{p} \rfloor} \times C_{n \bmod p}^{m \bmod p} \pmod{p}$
证明详见[OI-wiki](https://oi-wiki.org/math/number-theory/lucas/)
特别的，当 $p=2$ 时，有 $C_n^m \bmod 2 = [n \land m = m] = [n \lor m = n]$
代码：
```cpp
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
inline int C(int n,int m,int p){//fac是阶乘数组
if(n<0||m<0||n<m)return 0;
return fac[n]*qpow(fac[m],p-2,p)%p*qpow(fac[n-m],p-2,p)%p;;
}
inline int Lucas(int n,int m,int p){
if(!m)return 1;
return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}
```
## Part-7 ExCRT
~~为什么不讲 CRT？因为没用。~~
考虑一个形如：
$$\begin{cases}x\equiv b_1\pmod{a_1}\\x\equiv b_2\pmod{a_2}\\\dots\\x\equiv b_n\pmod{a_n}\end{cases}$$
的同余方程组，其中对于 $a_i$，$b_i$并没有什么限制条件。
我们考虑合并同余方程，假设前 $i-1$ 个方程的特解为 $r_1$，通解为 $r_1+tM (t \in \mathbb Z)$，易知 $M=LCM_{j=1}^{j=i-1} a_j$。则对于当前方程 $x \equiv b_i \pmod{a_i}$ 我们需要找一个 $t \in \mathbb Z$ 使得 $r_1+tM \equiv b_i \pmod{a_i}$，这个同余方程可以转化成 $ax+by=c$，用 Exgcd 解决，如果遇到 $gcd(a,b)\nmid c$ 的情况则无解。
代码：[P4777 【模板】扩展中国剩余定理（EXCRT）](https://www.luogu.com.cn/problem/P4777)
```cpp
#include <bits/stdc++.h>
#define int __int128
using namespace std;
const int maxn=1e5+5;
int n,m[maxn],a[maxn],ans,M;
inline int read(){
int ret=0,f=1;
char c=getchar();
while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
while(isdigit(c)){ret=ret*10+c-'0';c=getchar();}
return ret*f;
}
inline void write(int x){
if(x<0)x=-x,putchar('-');
if(x>9)write(x/10);
putchar(x%10+'0');
}
inline int exgcd(int a,int b,int &x,int &y){
if(!b){x=1,y=0;return a;}
int gcd=exgcd(b,a%b,x,y);
int t=x;x=y,y=t-(a/b)*y;
return gcd;
}
signed main(){
n=read();
for(int i=1;i<=n;i++)m[i]=read(),a[i]=read();
ans=a[1],M=m[1];
for(int i=2;i<=n;i++){
int mod=m[i],re=a[i];
int c=((re-ans)%mod+mod)%mod;
int x,y,gcd=exgcd(M,mod,x,y);
if(c%gcd){printf("-1\n");return 0;}
int t=c/gcd*x%mod;ans+=t*M;
M=mod/gcd*M;ans=(ans%M+M)%M;
}
write(ans);
return 0;
}
```
## Part-8 扩展欧拉定理
我们重拾一下快速幂的例题，如果那一题 $1 \leq a,b \leq 10^{2 \times 10^7}$ 怎么办？这时普通的快速幂难以解决问题，于是我们有请扩展欧拉定理登场。
~~其实也是一个公式，但我扯了一大堆废话。~~
$a^b \equiv \begin{cases}
a^b &b<\varphi(p) \\
a^{b \bmod \varphi(p)+\varphi(p)} &b \ge \varphi(p)
\end{cases}$
这个的证明并不难：
1.第一项显然。
2.令 $t=\varphi(p)$，$b=k\varphi(p)+r$
则 $a^b=a^{k\varphi(p)+r}$
由欧拉定理 $a^{\varphi(p)} \equiv 1 \pmod{p}$ 知 $a^{k\varphi(p)+r}=({a^{\varphi(p)}})^k \cdot a^r \equiv a^{r+\varphi(p)} \pmod{p}$，得证。
代码：[P5091 【模板】扩展欧拉定理](https://www.luogu.com.cn/problem/P5091)
```cpp
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int INF=1e9+7;
int a,m,b;bool flg;
inline int read(const int mod){
int ret=0,f=1;
char c=getchar();
while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
while(isdigit(c)){ret=ret*10+c-'0';if(mod!=INF&&ret>=mod)ret%=mod,flg=1;c=getchar();}
return ret*f;
}
inline int phi(int x){
int ret=x;
for(int i=2;i*i<=x;i++){
if(x%i==0){
ret=ret*(i-1)/i;
while(x%i==0)x/=i;
}
}
if(x>1)ret=ret*(x-1)/x;
return ret;
}
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
signed main(){
a=read(INF),m=read(INF),b=read(phi(m));
if(flg)b+=phi(m);
printf("%lld\n",quick_pow(a,b,m));
return 0;
}
```
## Part-9 原根
注：Part-9 的讨论均在模 $p$ 意义下进行，若想了解更多有关原根的知识及证明，请移步[codecode的题解](https://www.luogu.com.cn/article/py24mb57)。
首先我们引入“阶”的概念，记作 $\delta_p(a)$，她表示 $a^x \equiv 1 \pmod{p}$ 的最小正整数解。
阶的性质：
1.若 $a^n \equiv 1 \pmod{p}$，则有$\delta_p(a) \mid n$。
2.若 $gcd(a,p)=1$，则 $\delta_p(a^k)= \frac{\delta_p(a)}{gcd(\delta_p(a),k)}$
3.若 $gcd(a,p)=1$，$gcd(b,p)=1$，则 $\delta_p(ab)=\delta_p(a)\delta_p(b)$ 的充要条件是 $gcd(\delta_p(a),\delta_p(b))=1$
我们可以在 $O(\log^2 p)$ 的复杂度下求出一个数的阶，我们可以使用试除法。首先，我们可以判断 $\delta_p(a) \mid k$ 这样一个命题的真伪，即使用快速幂判断是否有 $a^k \equiv 1 \pmod{p}$。我们可以枚举一个数的质因子，如果出去之后条件满足，则除掉，最后剩下的就是阶。
关于阶，她还有一个性质。$a^x\equiv y \pmod{p}$ 有解的充要条件是 $\delta_p(y) \mid \delta_p(x)$ 。
代码：
```cpp
inline int Get_delta(int a){
int Tp=phi_p;
for(int i=1;i<=tot;i++)
while(Tp%key[i]==0&&qpow(a,Tp/key[i],p)==1)Tp/=key[i];//key[i]表示phi(p)的质因子
return Tp;
}
```
接下来我们讲一讲原根的定义，对于一个数 $g$，若满足 $gcd(g,p)=1$ 且 $\delta_p(g)=\varphi(p)$ 即为原根。
回到正题，接下来来看关于怎样的 $p$，其原根存在。
结论：当且仅当，$p=1,2,4,m^k,2 \times m^k$，其中 $p$ 为奇素数，$k$ 为正整数时有原根。
对于原根的求解，我们可以先考虑如何判断一个数是否是原根，我们令判断数为 $a=\prod_{i=1}^{n} p_i^{k_i}$，则当且仅当所有 $a^{\frac{\varphi(p)}{p_i}} \not\equiv 1 \pmod{p}$ 时 $a$ 为原根。
随后，我们可以求出最小原根，王元教授于1959年证明了她是 $O(n^{0.25})$ 级别的，复杂度为 $O(n^{0.25} \log^2p)$，设其为 $g$
然后，其他的原根都可以由 $g^k (k \in \mathbb Z)$ 得到，由阶的性质2，当 $gcd(\varphi(p),i)=1$ 时 $g^i$ 为原根。时间复杂度 $O(\varphi(p) \log \varphi(p))$
由上述可得，一个数 $p$ 如果有原根，则原根个数为 $\varphi(\varphi(p))$
代码：[【模板】原根](https://www.luogu.com.cn/problem/P6091)
```cpp
#include <bits/stdc++.h>
#define int long long
using namespace std;
int T,p,d,Phi,key[105],tot;
inline int read(){
int ret=0,f=1;
char c=getchar();
while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
while(isdigit(c)){ret=ret*10+c-'0';c=getchar();}
return ret*f;
}
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
inline int phi(int x){
int ret=x;tot=0;
for(int i=2;i*i<=x;i++){
if(x%i==0){
ret=ret*(i-1)/i,key[++tot]=i;
while(x%i==0)x/=i;
}
}
if(x>1)ret=ret*(x-1)/x,key[++tot]=x;
return ret;
}
inline bool check(int x){
if(qpow(x,Phi,p)!=1)return 0;
for(int i=1;i<=tot;i++)
if(qpow(x,Phi/key[i],p)==1)return 0;
return 1;
}
signed main(){
T=read();
while(T--){
p=read(),d=read();
Phi=phi(p);int c=phi(Phi),g=p;
for(int i=1;i<p;i++)
if(check(i)){g=i;break;}
if(g==p){printf("0\n\n");continue;}
vector<int> ans;ans.push_back(-1);
for(int i=1;i<=Phi;i++)
if(__gcd(i,Phi)==1)ans.push_back(qpow(g,i,p));
sort(ans.begin(),ans.end());
printf("%lld\n",c);
for(int i=1;i<=c/d;i++)printf("%lld ",ans[i*d]);
putchar('\n');
}
return 0;
}
```
## Part-10 BSGS 算法
模 $p$ 意义下的离散对数，记作 $\log_a b$，定义为 $a^x \equiv b \pmod{p}$ 的解。对于这样的方程且 $gcd(a,p)=1$ 的情况，我们可以使用大步小步算法(Baby Step Giant Step,简称 BSGS)求解。
本质上，BSGS 是一个优雅的暴力(分块)。我们令 $t=\lfloor \sqrt p \rfloor$，则有 $x=kt-r$，
则只需求 $a^{kt-r} \equiv b \pmod{p}$ 即 $(a^t)^k \equiv ba^r\pmod{p}$ 的解。我们将每个 $r \in [0,t]$ 对应的 $ba^r \bmod p$ 存入一个 gp_hash_table 中，然后在枚举 $k$ 进行查询即可。还有一些特判见代码吧，大家应该能理解。
离散对数也具有和对数一样的性质，即 $\log a+\log b=\log ab$。
代码：
```cpp
inline int BSGS(int a,int b,int p){
__gnu_pbds::gp_hash_table<int,int> Hash;
int t=sqrt(p)+1,tmp=b;
for(int i=0;i<=t;i++)Hash[tmp]=i,tmp=tmp*a%p;
int now=1;a=qpow(a,t,p);
if(!a)return !b?1:-1;
for(int i=1;i<=t;i++){
now=now*a%p;
if(Hash[now])return i*t-Hash[now];
}
return -1;
}
```
## Part-11 exLucas定理
当 $p$ 不为质数时，如何求 $C_n^m \bmod p$ 呢？显然，此时 Lucas 定理不再适用，于是我们有了 exLucas。
首先对 $p$ 分解质因数，有 $p=\prod_{i=1}^n p_i^{k_i}$，我们考虑求出每个 $C_n^m \bmod p_i^{k_i}$ 的值，再使用 ExCRT 合并。
现在，我们要求 $C_n^m \bmod p^k$ 的值，其中 $p$ 为素数。我们将 $C_n^m$ 展开：$C_n^m=\frac{n!}{m!(n-m)!}$。
显然，我们不能直接求解逆元，因为可能不存在，我们可以将原式化为这样：$\frac{n!}{m!(n-m)!} \bmod p^k=\frac{\frac{n!}{p^{k_1}}}{\frac{m!}{p^{k_2}} \cdot \frac{(n-m)!}{p^{k_3}}} \times p^{k_1-k_2-k_3} \bmod p^k$，其中 $k_1$，$k_2$，$k_3$ 分别表示 $n!$，$m!$，$(n-m)!$ 含有质因子 $p$ 的个数。现在，分母项中的数均与模数互质，可以直接求逆元。
我们再来考虑一下如何求 $\frac{n!}{p^a} \bmod p^k$ 的值。
举个例子，$n=22,p=3,k=2$，把这个写出来： $22!=1 \times 2 \times 3 … \times 22$，将其中所有 $p$ 的倍数提取出来，则有 $22!=3^7 \times (1 \times 2 … \times 7) \times (1 \times 2 \times 4 \times 5 \times 7 \times 8 \times 10 \times 11 \times 13 \times 14 \times 16 \times 17 \times 19 \times 20 \times 22)$，容易发现第一部分是 $p^{\lfloor \frac{n}{p} \rfloor}$，第二部分是 ${\lfloor \frac{n}{p} \rfloor}!$，第三部分是 $1$ 到 $n$ 中与 $n$ 互质的数的乘积。关于第三部分，有 $\prod_{i=1,gcd(i,p)=1}^{p^k} i \equiv \prod_{i=1,gcd(i,p)=1}^{p^k} (i+tp^k) \pmod{p^k}$，$\prod_{i=1,gcd(i,p)=1}^{p^k} i$ 一共循环了 $\lfloor \frac{n}{p^k} \rfloor$ 次，可以直接暴力求一个循环节在加上快速幂，最后在乘上零头的贡献。而显然第一部分与第二部分的 $p$ 因子都被除掉，无需考虑，于是可以直接递归 ${\lfloor \frac{n}{p} \rfloor}!$ 与当前答案相乘。
在将结论回代，就写好了。最后在补上一个 ExCRT 即可。(好像 CRT 也可以(大雾))
代码：[P4720 【模板】扩展卢卡斯定理/exLucas](https://www.luogu.com.cn/problem/P4720)
```cpp
#include <bits/stdc++.h>
#define int __int128
using namespace std;
const int maxn=1005;
int N,M,P,tot,m[maxn],a[maxn];
inline int read(){
int ret=0,f=1;
char c=getchar();
while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
while(isdigit(c)){ret=ret*10+c-'0';c=getchar();}
return ret*f;
}
inline void write(int x){
if(x<0)x=-x,putchar('-');
if(x>9)write(x/10);
putchar(x%10+'0');
}
inline int qpow(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
inline int exgcd(int a,int b,int &x,int &y){
if(!b){x=1,y=0;return a;}
int gcd=exgcd(b,a%b,y,x);
y-=a/b*x;return gcd;
}
inline int inv(int a,int p){
int x,y;exgcd(a,p,x,y);
return (x%p+p)%p;
}
inline int fac(int n,int p,int pk){
if(!n)return 1;
int ans=1;
for(int i=1;i<pk;i++)
if(i%p)ans=ans*i%pk;
ans=qpow(ans,n/pk,pk);
for(int i=1;i<=n%pk;i++)
if(i%p)ans=ans*i%pk;
return ans*fac(n/p,p,pk)%pk;
}
inline int C(int n,int m,int p,int pk){
if(n<m)return 0;
int fn=fac(n,p,pk),fm=fac(m,p,pk),fnm=fac(n-m,p,pk),cnt=0;
for(int i=n;i;i/=p)cnt+=i/p;
for(int i=m;i;i/=p)cnt-=i/p;
for(int i=n-m;i;i/=p)cnt-=i/p;
return fn*inv(fm,pk)%pk*inv(fnm,pk)%pk*qpow(p,cnt,pk)%pk;
}
inline int exCRT(){
int ans=a[1],M=m[1];
for(int i=2;i<=tot;i++){
int mod=m[i],re=a[i];
int c=((re-ans)%mod+mod)%mod;
int x,y,gcd=exgcd(M,mod,x,y);
int t=c/gcd*x%mod;ans+=t*M;
M=mod/gcd*M;ans=(ans%M+M)%M;
}
return ans;
}
inline int exLucas(int N,int M,int P){
tot=0;
for(int i=2;i*i<=P;i+=(i&1)+1){
int tmp=1;
while(P%i==0)P/=i,tmp*=i;
if(tmp!=1)a[++tot]=C(N,M,i,tmp),m[tot]=tmp;
}
if(P>1)a[++tot]=C(N,M,P,P),m[tot]=P;
return exCRT();
}
signed main(){
N=read(),M=read(),P=read();
write(exLucas(N,M,P)%P);
return 0;
}
```
## Part-12 exBSGS
当 $gcd(a,p) \ne 1$ 时，BSGS 算法将不能正确求解离散对数，于是，exBSGS 应要求而生。
不妨设 $g=gcd(a,p)$ 那么一定有 $\frac{a}{g} \cdot a^{x-1} \equiv \frac{b}{g} \pmod{\frac{p}{g}}$
当然，可能存在 $g \nmid b$，此时方程无解。
但即使能够整除，我们也无法保证 $gcd(a,\frac{p}{g})=1$，怎么办？
我们可以对上述操作进行重复迭代，直到 $gcd(a,\frac{p}{\prod_{i=1}^k g_i})=1$ 为止。即：
$\frac{a}{\prod_{i=1}^k g_i} \cdot a^{x-k} \equiv \frac{b}{\prod_{i=1}^k g_i} \pmod{\frac{p}{\prod_{i=1}^k g_i}}$
这个时候，两边同时乘以 ${(\frac{a}{\prod_{i=1}^k g_i})}^{-1}$，得：
$a^{x-k} \equiv {(\frac{a}{\prod_{i=1}^k g_i})}^{-1} \cdot \frac{b}{\prod_{i=1}^k g_i} \pmod{\frac{p}{\prod_{i=1}^k g_i}}$
这样一个方程，就可以使用传统的 BSGS 算法解决了！
代码：
```cpp
inline int exBSGS(int a,int b,int p){
a%=p,b%=p;if(b==1)return 0;if(!a)return !b?p!=1:-1;
int x,y,gcd=exgcd(a,p,x,y),k=0,tmp=1;
while(gcd!=1){
if(b%gcd)return -1;
k++,b/=gcd,p/=gcd,tmp=tmp*(a/gcd)%p;
if(tmp==b)return k;
gcd=exgcd(a,p,x,y);
}
int ans=BSGS(a,b*inv(tmp,p)%p,p);
if(ans==-1)return -1;else return ans+k;
}
```
## Part-13 二次剩余
二次剩余俗称模意义开根，定义为 $x^2 \equiv a \pmod{p}$ 的解，其中 $p$ 为奇素数。
首先我们先介绍一下 $Legendre$ 符号。定义如下：
$(\frac{a}{p}) = \begin{cases}
0 & p \mid a\\
1 & \exists x,x^2 \equiv a \pmod{p}\\
-1 & otherwise
\end{cases}$
接下来是欧拉准则，用于判断是否存在二次剩余。
当 $p$ 为奇素数时，由费马小定理有 $a^{p-1} \equiv 1 \pmod{p}$，即 $({a^2})^{\frac{p-1}{2}} \equiv 1 \pmod{p}$，可得 $(a^{p-1}-1)(a^{p-1}+1) \equiv 0 \pmod{p}$，所以，在模 $p$ 意义下，$a^{p-1} \bmod p$ 的取值只有两种，即 $1$ 和 $p-1$。
欧拉准则给出公式：$(\frac{a}{p})=a^{\frac{p-1}{2}}$。
我们来证明一下她。
首先，我们先证“$a$ 是模 $p$ 意义下的二次剩余”是“$a^{\frac{p-1}{2}} \equiv 1 \pmod{p}$”的充要条件。
- 充分性：
当 $a \equiv x^2 \pmod{p}$ 成立时，有 $a^{\frac{p-1}{2}} \equiv x^{p-1} \equiv 1 \pmod{p}$，证毕。
- 必要性：
设 $g$ 是 $p$ 的一个原根(由于 $p$ 是奇素数一定存在)，设 $a=g^k$，则有 $a^{\frac{p-1}{2}} \equiv 1 \equiv g^{\frac{k}{2}(p-1)} \pmod{p}$，由于 $g^{p-1} \equiv 1 \pmod{p}$，所以此时 $k$ 一定为偶数，令 $x=g^{\frac{k}{2}} \bmod p$ 则有 $x^2 \equiv a \pmod{p}$，证毕。
既然“$a$ 是模 $p$ 意义下的二次剩余”是“$a^{\frac{p-1}{2}} \equiv 1 \pmod{p}$”的充要条件，那么另一种情况自然对应了非二次剩余，剩下的一个也是正确的，于是，我们证明了欧拉准则。
对于解的数量，我们有在模 $p$ 意义下共有 $\frac{p-1}{2}$ 个二次剩余，且若 $(\frac{a}{p})=1$ 则两个解互为相反数(相加为 $p$)。这个非常好证。
至此，我们已经有了一种求二次剩余的方法，即先求出最小原根 $g$，在使用 BSGS 解 $g^k \equiv a \pmod{p}$ 的解，答案就是 $g^{\frac{k}{2}}$，时间复杂度 $O(\sqrt p)$。
但是，这种方法并不够优秀，我们有一种方法，可以 $O(\log p)$ 求解二次剩余，他就是 Cipolla 算法。
流程如下：
首先，我们先找出一个数 $n$ 满足 $n-a^2$ 是非二次剩余，令 $i^2 \equiv n-a^2 \pmod{p}$，这个 $i$ 就类似我们学过的复数单位，可以把一个数表示为 $a+bi$ 的形式。
这个找 $n$ 个过程可以直接随机化+快速幂判断，由于非二次剩余的比例接近 $50\%$，所以这样做期望 $2$ 次就可以找到 $n$。
那么，我们有 $(n+i)^{p+1} \equiv a \pmod{p}$，考虑证明。
我们先证两个引理：
引理1：$i^p \equiv -i \pmod{p}$
证明：$i^p \equiv (i^2)^{\frac{p-1}{2}}i \equiv (n-a^2)^{\frac{p-1}{2}}i \equiv -i \pmod{p}$。
引理2：$(A+B)^p \equiv A^p+B^p \pmod{p}$
证明：二项式定理展开后，中间项模 $p$ 全部为 $0$。
那么，就有 $(n+i)^{p+1} \equiv (n+i)^p(n+i) \equiv (n^p+i^p)(n+i) \equiv (n-i)(n+i) \equiv n^2-i^2 \equiv a^2 \pmod{p}$。
最后还有一个问题，$(n+i)^{p+1}$ 的虚部一定为 $0$ 吗？幸运的是，的确如此，可以使用反证法证明，在此不加赘述。
代码：[【模板】二次剩余](https://www.luogu.com.cn/problem/P5491)
```cpp
#include <bits/stdc++.h>
#define int long long
using namespace std;
int T,mod,n,p,ImI;
struct Complex{
int real,imag;
bool operator ==(const Complex &B){return real==B.real&&imag==B.imag;}
Complex operator *(const Complex &B){
Complex C;
C.real=(real*B.real%mod+ImI*imag%mod*B.imag%mod+mod)%mod;
C.imag=(real*B.imag%mod+imag*B.real%mod+mod)%mod;
return C;
}
};
inline int read(){
int ret=0,f=1;
char c=getchar();
while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
while(isdigit(c)){ret=ret*10+c-'0';c=getchar();}
return ret*f;
}
inline int qpow_r(int a,int b,int p){
int ret=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)ret=ret*a%p;
return ret;
}
inline Complex qpow_i(Complex a,int b){
Complex ret={1,0};
for(;b;b>>=1,a=a*a)
if(b&1)ret=ret*a;
return ret;
}
inline bool check(int n,int p){return qpow_r(n,p-1>>1,p)==p-1;}
inline void Cipolla(int n,int p){
if(!n)return printf("0\n"),void();n%=p;
if(qpow_r(n,p-1>>1,p)==p-1)return printf("Hola!\n"),void();
mod=p;int a=rand()%p;
while(1){
a=rand()%p,ImI=((a*a%p-n)%p+p)%p;
if(check(ImI,p))break;
}
Complex ans=qpow_i({a,1},p+1>>1);
int ans1=(ans.real%p+p)%p,ans2=p-ans1;
if(ans1==ans2)return printf("%lld\n",ans1),void();
if(ans1>ans2)swap(ans1,ans2);
printf("%lld %lld\n",ans1,ans2);
}
signed main(){
srand(time(0));
T=read();
while(T--){
n=read(),p=read();
Cipolla(n,p);
}
return 0;
}
```
## Part-14 数论分块
数论分块，又称整除分块，核心思想在于利用整除相同值具有一段连续区间的特性优化暴力。
我们来考虑数论分块的基础问题：
给定 $n$，求 $\sum_{i=1}^n \lfloor \frac{n}{i} \rfloor$，满足 $1 \leq n \leq 10^9$
我们不妨令 $t=\sqrt n$，则当 $i\leq t$ 时，满足 $\lfloor \frac{n}{i} \rfloor \leq \sqrt n$，即只有 $\sqrt n$ 个取值，当 $i>t$ 时同理。于是 $\lfloor \frac{n}{i} \rfloor \leq \sqrt n$ 只有 $O(\sqrt n)$ 个取值。
我们考虑对取值分块，我们枚举满足 $\lfloor \frac{n}{i} \rfloor$ 的连续段区间 $[l,r]$，当我们已经知道 $l$ 时，我们可以快速算出 $r=\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor} \rfloor$，统计答案时只需加上贡献值以及贡献出现次数即可。
代码：
```cpp
for(int l=1,r=1;l<=n;l=r+1)r=n/(n/l),ans+=(r-l+1)*(n/l);
```
## Part-15 莫比乌斯反演
在学习莫比乌斯反演之前，我们先来学习一下狄利克雷卷积。
- 定义：对于两个数论函数 $f(x)$ 与 $g(x)$，她们的狄利克雷卷积得到的结果 $h(x)$ 定义为 $h(x)=\sum_{d \mid x} f(d)g(\frac{n}{d})=\sum_{ab=x} f(a)g(b)$。上式可以简记为 $h=f*g$。
- 性质：
交换律：$f*g=g*f$
结合律：$f*g*h=f*(g*h)$
分配律：$(f+g)*h=f*h+g*h$
接下来，我们定义几个新函数：
- $\varepsilon$ 函数：$\varepsilon(n)=[n=1]$(又被称为单位元，性质是 $\varepsilon * f=f$)
- $id^k$ 函数：$id^k(n)=n^k$
- $\sigma_k$ 函数：$\sigma_k(n)=\sum_{d\mid n} d^k=(id^k * 1)(n)$
以及积性函数的定义：
- 若数论函数 $f(x)$ 满足 $\forall gcd(a,b)=1,a \in I,b \in I,f(ab)=f(a)f(b)$，则 $f(x)$ 被称为积性函数。其中，若 $f(x)$ 满足 $\forall a \in I,b \in I,f(ab)=f(a)f(b)$，则 $f(x)$ 被称为完全积性函数。
还有狄利克雷卷积的重要结论：
- 两个积性函数的狄利克雷卷积仍是积性函数。
接下来，我们就可以开始讲莫比乌斯反演了。
- 定义：
$\mu$ 为莫比乌斯函数，定义为
$\mu(n) = \begin{cases}
1 & n=1 \\
0 & \exists x,x^2 \mid n\\
(-1)^k & n=\prod_{i=1}^k p_i
\end{cases}$
- 性质：
$\mu(x)$ 为积性函数。
$\mu*1=\varepsilon$
反演结论：$[n=1]=\sum_{d \mid n} \mu(d)$
根据她的性质及积性函数，可以使用欧拉筛 $O(n)$ 计算莫比乌斯函数的值，代码：
```cpp
inline void make_Mu(int n){
mu[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i])prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
}
```
当然，这样的好东西也少不了欧拉函数的参与，对于 $\varphi$ 函数，我们有 $\varphi*1=id$ 和 $\sum_{d \mid n} \varphi(d)=n$，这些都有大用。
接下来，我们一起推导一下莫比乌斯反演的例子：
- 形式 $1$:
$$\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^m [gcd(i,j)=1]&=\sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{x \mid j} \mu(x)\\
&=\sum_{x=1}^n \mu(x)\sum_{i=1}^n\sum_{j=1}^m [x \mid i][x \mid j]\\
&=\sum_{x=1}^n \mu(x)\lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor
\end{aligned}$$
- 形式 $2$：
$$\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^m gcd(i,j) &=\sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{x \mid j} \varphi(x)\\
&=\sum_{x=1}^n \varphi(x)\sum_{i=1}^n\sum_{j=1}^m [x \mid i][x \mid j]\\
&=\sum_{x=1}^n \varphi(x)\lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor
\end{aligned}$$
- 形式 $3$：
$$\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^m f(gcd(i,j))&=\sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{x \mid j} (f*\mu)(x)\\
&=\sum_{x=1}^n (f*\mu)(x)\sum_{i=1}^n\sum_{j=1}^m [x \mid i][x \mid j]\\
&=\sum_{x=1}^n (f*\mu)(x)\lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor\\
&=\sum_{x=1}^n \lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor \sum_{d \mid x}f(x)\mu(\frac{x}{d})
\end{aligned}$$
以上都可以使用线性筛 $O(n)$ 预处理，数论分块 $O(\sqrt n)$ 查询。