专栏文章

数论学习笔记

TTH911
算法·理论参与者 11已保存评论 12

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
12 条
当前快照
1 份
快照标识符
@mio8zaiy
此快照首次捕获于
2025/12/02 15:19
3 个月前
此快照最后确认于
2025/12/02 15:19
3 个月前
查看原文
大概率更好的阅读体验
本文可能会不定期更新。若存在内容正确性上的问题，请联系修复。
更新/修正记录
  • 2025/9/132025/9/13 更新：修正 excrt 的参考代码。
  • 2025/11/282025/11/28 更新：修正杜教筛的推导过程。
本文中出现的代数式、数，如无特殊说明，均为整数。
对于 OI 中的数论题，答案不对请先尝试：
CPP
#define int __int128
typedef long long ll;
#define ll __int128
数论函数指的是定义域为正整数的函数，也可以视为一个序列。

基础部分

⊥\perp 在数论中表示“互质”，例如 2⊥32\perp3。
对于一个序列 fn=afn−1+bfn−2f_n=af_{n-1}+bf_{n-2}，若 gcd⁡(a,b)=1\gcd(a,b)=1，则有：
gcd⁡(fx,fy)=fgcd⁡(x,y)\gcd(f_x,f_y)=f_{\gcd(x,y)}

质数

Fermat 素性测试

由费马小定理，对于质数 pp，若 gcd⁡(a,p)=1\gcd(a,p)=1，满足 ap−1≡1(modp)a^{p-1}\equiv 1\pmod p，ap−2≡a−1(modp)a^{p-2}\equiv a^{-1}\pmod p。
但是 ap−1≡1(modp)a^{p-1}\equiv1\pmod p 并不能推出 pp 为质数。
因此可以随机选择来测试。

Miller–Rabin 素性测试

Miller-Rabin 素性测试，取质数集合 A={2,3,5,7,11,13,17,19,23,29,31,37}A=\lbrace{2,3,5,7,11,13,17,19,23,29,31,37}\rbrace，则可以通过随机化确定 long long 范围（[0,264)[0,2^{64})）内的任意整数 nn 是否为质数。
从 AA 中取一底数 aa，若：
  • n=an=a，nn 为质数。
  • n mod a=0∧n>an\bmod a=0\land n>a，nn 不为质数。
在都不成立的情况下，进行 Miller-Rabin 素性测试。
将 n−1n-1 转化：
n−1=d×2rn-1=d\times2^r
其中，d,r∈N∗d,r\in\N^*，也就是说，dd 是 n−1n-1 在二进制位上的一个前缀，满足前缀的后缀不为 00。

二次探测定理

x2≡1(modp)⇓x≡±1(modp)x^2\equiv1\pmod p\\ \Downarrow\\ x\equiv\pm1\pmod p
使用平方差公式可证。

Miller-Rabin 素性测试的实现

基于 aa 执行 rr 轮测试，通过二次探测定理判断。

参考代码

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;
}

最大公约数

裴蜀定理

对于 ∀a,b∈Z\forall a,b\in \Z：
  • ∃x,y∈Z\exist x,y\in \Z，使 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b)。
  • ∀x,y∈Z\forall x,y\in\Z，gcd⁡(a,b)∣(ax+by)\gcd(a,b)\mid(ax+by)。

逆定理

∀a,b∈Z\forall a,b\in \Z，若 d>0d>0 且 dd 为 a,ba,b 的公因数，且存在 x,yx,y，使得 ax+by=dax+by=d，则：
d=gcd⁡(a,b)d=\gcd(a,b)

乘法逆元

若存在 ax≡1(modp)ax\equiv 1\pmod p，则称 xx 为 aa 关于 pp 的逆元，记 x=a−1x=a^{-1}。
当 gcd⁡(a,p)=1\gcd(a,p)=1 的时候（pp 为质数）逆元一定存在。

费马小定理

对于质数 pp，若 gcd⁡(a,p)=1\gcd(a,p)=1，满足 ap−1≡1(modp)a^{p-1}\equiv 1\pmod p，ap−2≡a−1(modp)a^{p-2}\equiv a^{-1}\pmod p。
常用于分数取模，即：
ab≡abp−2(modp)\frac{a}{b}\equiv ab^{p-2}\pmod p
详见费马小定理。

欧拉定理

欧拉函数
令 n=p1c1p2c2p3c3⋯pkckn=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_k^{c_k}，满足 pip_i 为质数，ci>0c_i>0。
令 n=p1c1p2c2p3c3⋯pkckn=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_k^{c_k}，满足 pip_i 为质数，ci>0c_i>0。
则：
φ(n)=n(1−1p1)(1−1p2)(1−1p3)⋯(1−1pk) \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)
推导见容斥原理求解欧拉函数。（组合数学没搬，放的原 Blog 链接）
若 gcd⁡(a,p)=1\gcd(a,p)=1，则 aφ(p)≡1(modp)a^{\varphi(p)}\equiv1\pmod p。
显然，当 pp 为质数时，有 φ(p)=p−1\varphi(p)=p-1，即费马小定理。因此，费马小定理为欧拉定理的特殊形式。

扩展欧拉定理

对于任意的 a,p,b∈N∗a,p,b\in \N^*，若满足 gcd⁡(a,p)>1,b≥φ(p)\gcd(a,p)>1,b\geq \varphi(p)，有：
ab≡ab mod φ(p)+φ(p)(modp)a^b\equiv a^{b\bmod \varphi(p)+\varphi(p)}\pmod p
当 gcd⁡(a,p)=1\gcd(a,p)=1 时即欧拉定理。
可以使用初等证明或群论证明。
常用于降幂。

扩展欧几里得算法 exgcd

注意不是“类欧几里得算法”。
用于求关于 x,yx,y 的不定方程 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b) 的一组整数特解。
不妨令 a>ba>b。
考虑到 a=⌊ab⌋b+a mod ba=\left\lfloor\dfrac{a}{b}\right\rfloor b+a\bmod b，代入可得：
(⌊ab⌋b+a mod b)x+by=gcd⁡(a,b)\left(\left\lfloor\frac{a}{b}\right\rfloor b+a\bmod b\right)x+by=\gcd(a,b)
因此：
b(⌊ab⌋x+y)+(a mod b)x=gcd⁡(b,a mod b)b\left(\left\lfloor\frac ab\right\rfloor x+y\right)+(a\bmod b)x=\gcd(b,a\bmod b)
令 a′=b,x′=⌊ab⌋x+y,b′=a mod b,y′=xa'=b,x'=\left\lfloor\dfrac ab\right\rfloor x+y,b'=a\bmod b,y'=x，有：
a′x′+b′y′=gcd⁡(a′,b′)a'x'+b'y'=\gcd(a',b')
显然，可以递归求解。
那么直到 b=0b=0 时，可以直接解得一组整数特解 {x=1y=0\begin{cases}x=1\\y=0\end{cases}。
设最终递归出来的解为 {x=x0y=y0\begin{cases}x=x_0\\y=y_0\end{cases}。
令 Δb=∣bgcd⁡(a,b)∣，Δa=∣agcd⁡(a,b)∣\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert，\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert，则对于 ∀k∈N∗\forall k\in \N^*，方程存在一组解为 {x=x0+kΔby=y0−kΔa\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 求乘法逆元

如果需要讨论 aa 模 pp 意义下的乘法逆元 a−1a^{-1}，显然有 gcd⁡(a,p)=1\gcd(a,p)=1。
那么可列方程 ax+py=1ax+py=1，显然 ax≡1(modp)ax\equiv1\pmod p。
则通过扩展欧几里得算法求出 xx，即求出了 aa 的逆元。
例题链接：[NOIP 2012] 同余方程
求出 xx 后，为了使 xx 为正整数，只需要让 x←(x mod p+p) mod px\leftarrow(x\bmod p+p)\bmod p 即可。
参考代码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;
}

例题：正整数最大/小解

例题链接
给定不定方程 ax+by=cax+by=c，a,b,c∈N∗a,b,c\in \N^*，判断其是否有解。
若有正整数解，分别求 x,yx,y 的最大/小值。
否则，求 x,yx,y 得最小正整数解。
对于无解，即 gcd⁡(a,b)∤c\gcd(a,b)\nmid c，很好判断。
先通过 exgcd 求出一组特解 (x0,y0)(x_0,y_0)，则 {x=x0+kΔby=y0−kΔa\begin{cases}x=x_0+k\Delta b\\y=y_0-k\Delta a\end{cases} 也是原方程的解。Δb\Delta b 为 bb 的单次变化量，最小为 Δb=∣bgcd⁡(a,b)∣\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert，Δa=∣agcd⁡(a,b)∣\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert。
先求 xx 的最小正整数解 xmin⁡x_{\min}，即 x0+kΔbx_0+k\Delta b 最小，显然 xmin⁡=x0 mod Δbx_{\min}=x_0\bmod \Delta b。
但是为了保证 xmin⁡x_{\min} 为正整数，因此当 x0 mod Δb≤0x_0\bmod \Delta b\leq0 时，xmin⁡=x0 mod Δb+Δbx_{\min}=x_0\bmod \Delta b+\Delta b。
a,b>0a,b>0，因此 xx 最小时 yy 最大，因此可以确定 ymax⁡=c−axmin⁡by_{\max}=\dfrac{c-ax_{\min}}{b}。
同理，可求出 xmax⁡,ymin⁡x_{\max},y_{\min}。
当 xmin⁡x_{\min} 为正整数，此时若 ymax⁡>0y_{\max}>0，则存在正整数解 {x=xmin⁡y=ymax⁡\begin{cases}x=x_{\min}\\y=y_{\max}\end{cases}，因此可以判断正整数解（判断 xmax⁡>0x_{\max}>0 也行）。
正整数解的个数即：
xmax⁡−xmin⁡Δb+1\frac{x_{\max}-x_{\min}}{\Delta b}+1
因为从 xmin⁡x_{\min} 开始，xmin⁡+kΔbx_{\min}+k\Delta b 对应一组 x,yx,y。
参考代码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+by=\gcd(a,b)
而同余方程形如：
ax≡c(modb)ax\equiv c\pmod b
可以将同余方程转化为：
ax+by=cax+by=c
进而使用 exgcd 求解。

线性求逆元

有一种方法，可以在 O(log⁡V+n)\mathcal O(\log V+n) 的时间内求出任意 nn 个数 a1,a2,⋯ ,ana_1,a_2,\cdots,a_n 关于 pp 的逆元。
记：
prei≡∏j=1iaj(modp)\textit{pre}_i\equiv\prod_{j=1}^ia_j\pmod p
那么就可以 O(log⁡V)\mathcal O\left(\log V\right) 计算：
preinvn≡1pren(modp)\textit{preinv}_{n}\equiv\frac{1}{\textit{pre}_n}\pmod p
n−1∼1n-1\sim 1 递推，可以推出：
preinvi≡1prei≡1prei+1⋅ai+1(modp)\textit{preinv}_i\equiv\frac{1}{\textit{pre}_i}\equiv\frac{1}{\textit{pre}_{i+1} }\cdot a_{i+1}\pmod p
随后再次递推，可以得到 aia_i 关于 pp 的逆元 invi\textit{inv}_i：
invi≡preinvi⋅prei−1(modp)\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 中国剩余定理

CRT 中的运算溢出
在使用 CRT 或 exCRT 时，一律建议使用 __int128，防止溢出。
尽管题目大多会保证 lcm⁡(p1,p2,⋯ ,pn)≤1018\operatorname{lcm}(p_1,p_2,\cdots,p_n)\leq10^{18}，但是中间计算时会溢出。
CRT（中国剩余定理）被用于求解线性同余方程组：
{x≡a1(modp1)x≡a2(modp2)⋯x≡an(modpn)\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}
其中，p1,p2,p3,⋯ ,pnp_1,p_2,p_3,\cdots,p_n 两两互质。
令 L=∏i=1npiL=\prod\limits_{i=1}^np_i，则对于 k∈Nk\in\N 有 k⋅L≡0(modpi)k\cdot L\equiv0\pmod{p_i}。
记 qi=Lpiq_i=\dfrac{L}{p_i}，设 ci≡qi−1(modpi)c_i\equiv q_i^{-1}\pmod{p_i}，即 qiq_i 模 pip_i 意义下的逆元。
则，方程组的最小整数解为：
x=(∑i=1naiqici) mod Lx=\left(\sum_{i=1}^na_iq_ic_i\right)\bmod L
同时，对于 ∀k∈N∗\forall k\in\N^*，x+kLx+kL 均为原方程组的解。
CRT 其实就是一个构造式的做法，易证 ∀i≠j,aiqici≡0(modpj)\forall i\neq j,a_iq_ic_i\equiv0\pmod{p_j}。

例题：中国剩余定理（CRT）/ 曹冲养猪

题目中的 ai,bia_i,b_i 即分别为 pi,aip_i,a_i。
直接套 CRT 即可，但是需要注意的是，计算 aiqicia_iq_ic_i 时需要 __int128，会爆 long long。
参考代码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

例题链接
exCRT 实际上与 CRT 关系不大，如同 exLucas 与 Lucas 一般。
exCRT 解决的是当模数 p1,p2,p3,⋯ ,pnp_1,p_2,p_3,\cdots,p_n 不一定两两互质的时候（这时有可能无解）。
先考虑两个同余方程：
{x≡a1(modp1)x≡a2(modp2)\begin{cases} x\equiv a_1\pmod{p_1}\\ x\equiv a_2\pmod{p_2}\\ \end{cases}
设 r,sr,s，将其转化为两个不定方程：
{x=a1+r⋅p1x=a2+s⋅p2\begin{cases} x=a_1+r\cdot p_1\\ x=a_2+s\cdot p_2\\ \end{cases}
因此可得 r⋅p1−s⋅p2=a2−a1r\cdot p_1-s\cdot p_2=a_2-a_1。
由裴蜀定理，若 (a2−a1)∤gcd⁡(p1,−p2)(a_2-a_1)\nmid\gcd(p_1,-p_2)，则该不定方程无解，则原同余方程组无解。
否则通过 exgcd 求出 r,sr,s，则：
x=a1+r⋅p1=a2+s⋅p2x=a_1+r\cdot p_1=a_2+s\cdot p_2
既然 x=a1+r⋅p1x=a_1+r\cdot p_1，因此有：
x≡a1+r⋅p1(modlcm⁡(p1,p2))x\equiv a_1+r\cdot p_1\pmod{\operatorname{lcm}(p_1,p_2)}
这样取 lcm⁡(p1,p2)\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';
    }
}
参考代码
并没有判断无解。
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

即形如：
{f1x≡a1(modp1)f2x≡a2(modp2)f3x≡a3(modp3)⋯fnx≡an(modpn)\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}
多了一个系数 fif_i。
同样可以用类似 exCRT 的方法两两合并解决，具体见此。

离散对数

所谓离散对数，即模意义下取对数。
形如：给定模数 pp，及整数 a,ba,b，求整数 xx 使得 ax≡b(modp)a^x\equiv b\pmod p。
注意：离散对数可能不存在。

BSGS 算法

在 OI 中，BSGS 算法（Baby-Step Giant-Step，大步小步算法北上广深算法）常用来求解离散对数。
BSGS 算法要求 a⊥pa\perp p，求 xx 使得 ax≡b(modp)a^x\equiv b\pmod p。
若有解，则存在 x≤φ(p)x\leq\varphi(p) 的解；因为欧拉定理说明，a⊥pa\perp p 时，aφ(p)≡1(modp)a^{\varphi(p)}\equiv1\pmod p。
若枚举 xx 是 O(φ(p))\mathcal O(\varphi(p)) 的，当 pp 为质数的时候不能接受，因此可以考虑分块。
令 B=⌈φ(p)⌉B=\left\lceil\sqrt{\varphi(p)}\right\rceil，x=qB+rx=qB+r，且 0≤q,r≤B0\leq q,r\leq B。
则有：
aqB+r≡b(modp)a^{qB+r}\equiv b\pmod p
a⊥pa\perp p，则 aa 的乘法逆元 a−1a^{-1} 一定存在，有：
aqB≡b(a−1)r(modp)a^{qB}\equiv b\left(a^{-1}\right)^r\pmod p
枚举 rr，将其与其对应的 b(a−1)rb\left(a^{-1}\right)^r 一同存储在数据结构（一般是哈希表）中，随后枚举 qq，在数据结构中找 aqBa^{qB} 对应的 rr，找到了 rr 便找到了一个解 x=qB+rx=qB+r。同时，因为 a⊥pa\perp p，因此 a−1,a−2,a−3,⋯ ,a−Ba^{-1},a^{-2},a^{-3},\cdots,a^{-B} 模 pp 意义下互不相同。
时间复杂度：O(φ(p))\mathcal O\left(\sqrt{\varphi(p)}\right)，在 pp 为质数时最劣，O(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";
例题链接
参考代码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⊥pa\perp p 不成立的情况。
由扩展欧拉定理可知，若 xx 有解，则 x≤2φ(p)x\leq2\varphi(p) 范围内也有解。
首先特判掉 x=0x=0 的情况，即 b=1b=1 或 p=1p=1。
考虑分块，令 B=⌈2φ(p)⌉,x=qB−r,0≤q,r≤BB=\left\lceil\sqrt{2\varphi(p)}\right\rceil,x=qB-r,0\leq q,r\leq B，则有：
aqB−r≡b(modp)⇓aqB≡b⋅ar(modp)a^{qB-r}\equiv b\pmod p\\ \Downarrow\\ a^{qB}\equiv b\cdot a^r\pmod p
因此可以预处理 aqB mod pa^{qB}\bmod p，随后枚举 rr；但是注意到前者是后者的充分条件而不是充要条件，因此找出来的 x=qB+rx=qB+r 只是“可能的解”，需要检验。
对于多个不同的 qq 产生的模 pp 意义下相同的 aqBa^{qB}，可以只存储 qq 最小的两个。
证明
由扩展欧拉定理，axa^x 是在 k≤φ(p)k\leq\varphi(p) 值后以 φ(p)\varphi(p) 为周期循环的。
循环周期内的 aqBa^{qB} 显然至多存储 11 个，而非循环部分 aqBa^{qB} 也至多存储 11 个。
即：保留最小的 22 个。
时间复杂度：仍然为 O(φ(p))\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";
}
例题链接
参考代码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 定理的内容很简单：
(nm)≡(⌊np⌋⌊mp⌋)(n mod pm mod p)(modp)\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=∑i=0p(pi)aibp−i(a+b)^p=\sum_{i=0}^p\binom pia^ib^{p-i}
考虑 (pi)\dbinom pi 在模 pp 意义下的取值。
因为：
(pi)=p!i!(p−i)!\dbinom pi=\dfrac{p!}{i!(p-i)!}
那么如果化简之后，分子 p!p! 中的 pp 项没有被约分掉，则有 (pi)≡p⋅(p−1)!i!(p−i)!≡0(modp)\dbinom pi\equiv p\cdot\dfrac{(p-1)!}{i!(p-i)!}\equiv0\pmod p。
因为 pp 为质数，所以 pp 项能被约分掉当且仅当 i!i! 中含有 pp 项或 (p−i)!(p-i)! 中含有 pp 项。
即：(pi)≢0(modp)\dbinom pi\not\equiv0\pmod p 当且仅当 i≡0(modp)i\equiv0\pmod p 或 p−i≡0(modp)p-i\equiv 0\pmod p。
在二项式定理中，ii 满足 0≤i≤p0\leq i\leq p，所以 (pi)≢0\dbinom pi\not\equiv0 当且仅当 i=0i=0 或 i=pi=p。
在这两种情况中，都可以计算得到 (pi)=1\dbinom pi=1，即 (pi)≡[i=0∨i=p]\dbinom pi\equiv[i=0\lor i=p]。
重新带回二项式定理，可得：
(a+b)p≡(p0)a0bp+(pp)apb0≡ap+bp(modp)\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 mod p(1+x)^n\bmod p。
(1+x)n≡(1+x)p⌊np⌋+n mod p≡(1+x)p⌊np⌋(1+x)n mod p≡(1+xp)⌊np⌋(1+x)n mod p(modp)\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(1+x)^n 的 xmx^m 项系数为 (nm)\dbinom nm。
而想要从 (1+xp)⌊np⌋(1+x)n mod p(1+x^p)^{\lfloor\frac np\rfloor}(1+x)^{n\bmod p} 中得到 xmx^m，即从 (1+xp)⌊np⌋(1+x^p)^{\lfloor\frac np\rfloor} 中选取 ⌊mp⌋\lfloor\frac mp\rfloor 个 xpx^p，再从 (1+x)n mod p(1+x)^{n\bmod p} 中选取 m mod pm\bmod p 个 xx。
即：
(nm)≡(⌊np⌋⌊mp⌋)(n mod pm mod p)(modp)\dbinom nm\equiv\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p} \pmod p
你可能有一个问题：这看起来明明应当是一个等式，但为什么是同余呢？
即，Lucas 定理应当表述为：
(nm)=(⌊np⌋⌊mp⌋)(n mod pm mod p)\dbinom nm=\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}
但是，你要知道，只有在模 pp 意义下才有 (1+x)p⌊np⌋=(1+xp)⌊np⌋(1+x)^{p\lfloor\frac np\rfloor}=(1+x^p)^{\lfloor\frac np\rfloor}。
因此，只有在模 pp 意义下才有：
(nm)=(⌊np⌋⌊mp⌋)(n mod pm mod p)\dbinom nm=\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}
即：
(nm)≡(⌊np⌋⌊mp⌋)(n mod pm mod p)(modp)\dbinom nm\equiv\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\dbinom{n\bmod p}{m\bmod p}\pmod p

应用

当 nn 比较大而无法使用其他方法（例如预处理 1∼n1\sim n 的阶乘再利用乘法逆元）直接求解组合数时，可以使用 Lucas 定理。
Lucas 定理只需要递归使用即可，即递归计算 (⌊np⌋⌊mp⌋)\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}，递归终点即 (00)\dbinom{0}{0}。
其时间复杂度为：O(f(p)log⁡m)\mathcal O\left(f(p)\log m\right)。
其中，f(p)f(p) 表示单次计算 (n mod pm mod p)\dbinom{n\bmod p}{m\bmod p} 的复杂度，因写法而异。
可以使用乘法逆元，则时间复杂度为 O(log⁡plog⁡m)\mathcal O(\log p\log m)。
也可以 O(plog⁡p)\mathcal O(p\log p) 递推，时间复杂度为 O(plog⁡plog⁡m)\mathcal O(p\log p\log m)。
推荐使用乘法逆元。
例题链接
很简单，注意是 (n+mn)\dbinom{n+m}{n} 即可。
参考代码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 算法可以求 (nm) mod P\dbinom{n}{m}\bmod P，PP 不一定为质数。（但是，exLucas 实际上和 Lucas 没有多大关系……）
由唯一分解定理，可以对 PP 进行质因数分解：
P=p1c1p2c2p3c3⋯pncnP=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_n^{c^n}

CRT 求解

我们如果能够求出 a1,a2,⋯ ,ana_1,a_2,\cdots,a_n 使得：
{(nm)≡a1(modp1c1)(nm)≡a2(modp1c2)⋯(nm)≡an(modpncn)\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 求解出 (nm) mod P\dbinom nm\bmod P。因为在 CRT 中，恰好也有 P=p1c1p2c2⋯pncnP=p_1^{c_1}p_2^{c_2}\cdots p_n^{c_n}。

求解模质数幂下余数

展开定义式：
(nm)=n!m!(n−m)!\dbinom nm=\frac{n!}{m!(n-m)!}
因此：
(nm)≡n!m!(n−m)!(modpici)\dbinom nm\equiv\frac{n!}{m!(n-m)!}\pmod{p_i^{c_i}}
因为 pici,m!,(n−m)!p_i^{c_i},m!,(n-m)! 不一定互质，因此乘法逆元不一定存在。
考虑先提出分子和分母中所有的 pip_i 次幂，随后便可以使用逆元求解。
记 xx 的质因数分解中质数 pp 的幂次为 νp(x)\nu_p(x)，剩余积为 (x)p(x)_p，即：
x=pνp(x)⋅(x)px=p^{\nu_p(x)}\cdot(x)_p
则有：
n!m!(n−m)!≡(n!)pi(m!)pi((n−m)!)pi⋅piνpi(n!)−νpi(m!)−νpi((n−m)!)(modpici)\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,px,p，如何高效地求出 νp(x!) mod pici,(x!)p mod pici\nu_{p}(x!)\bmod p_i^{c_i},(x!)_p\bmod p_i^{c_i}。
考虑：
x!=1×2×⋯×p×⋯×2p×⋯×⌊xp⌋p×⋯×xx!=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,⋯ ,⌊xp⌋pp,2p,3p,\cdots,\left\lfloor\dfrac xp \right\rfloor p 中 pp 的个数为 ⌊xp⌋+νp(⌊xp⌋!)\left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)。其中 νp(⌊xp⌋!)\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right) 是 1,2,3,⋯ ,⌊xp⌋1,2,3,\cdots,\left\lfloor\dfrac xp\right\rfloor 中的 pp 的个数。
pp 为质数，所以在 1,2,3,⋯ ,p−1,p+1,⋯1,2,3,\cdots,p-1,p+1,\cdots 中不会含有 pp。
将递推式展开：
νp(x!)=⌊xp⌋+νp(⌊xp⌋!)=⌊xp⌋+⌊xp2⌋+νp(⌊xp2⌋!)=⌊xp⌋+⌊xp2⌋+⌊xp3⌋+⋯\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}
因此可以 O(log⁡px)\mathcal O(\log_p x) 计算 ν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 mod pc(x!)_p\bmod p^c；显然不能利用定义式 (x!)p=x!νp(x!)(x!)_p=\dfrac{x!}{\nu_p(x!)}，而需要其他方法（因为无法得知 x!x!）。
不难进行递推：
(x!)p≡∏i=1n(i)p≡(∏1≤i≤x,i⊥p(i)p)(∏i=1⌊xp⌋(p⋅i)p)≡(∏1≤i≤x,i⊥p(i)p)(⌊xp⌋!)p≡(∏1≤i≤pc,i⊥pi)⌊xpc⌋(∏1≤i≤x mod pc,i⊥pi)(⌊xp⌋!)p\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,pc,2pcm=2,4,p^c,2p^c 时 k≡−1k\equiv-1，否则为 11）：
∏1≤k≤m,k⊥mk≡±1(modm)\prod_{1\leq k\leq m,k\perp m}k\equiv\pm1\pmod m
因此：
(x!)p≡(±1)⌊xpc⌋(∏1≤i≤x mod pc,i⊥pi)(⌊xp⌋!)p\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(x!)_p 时，pcp^c 是固定的，因此可以预处理：
fj≡∏1≤i≤j,i⊥pi(modp)cf_j\equiv\prod_{1\leq i\leq j,i\perp p}i\pmod p^c
因此，得到最终推导式：
(x!)p≡(±1)⌊xpc⌋fx mod pc(⌊xp⌋!)p\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}
可以递归或迭代处理，时间复杂度 O(pc+log⁡nx)\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;
}
例题链接
参考代码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
*/
[国家集训队] 礼物
参考代码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;
}

数论分块

数论分块（整除分块）用于快速计算：
∑i=1nf(i)g(⌊ni⌋)\sum_{i=1}^nf(i)g\left(\left\lfloor\dfrac ni\right\rfloor\right)
数论分块的原理即本质不同的 ⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 只有至多 O(n)\mathcal O\left(\sqrt n\right) 种。
证明
  • 当 i≤ni\leq\sqrt n 时，ii 只有 O(n)\mathcal O\left(\sqrt n\right) 种取值。
  • 当 i>ni>\sqrt n 时，ni≤n\dfrac ni\leq\sqrt n，⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 只有 O(n)\mathcal O\left(\sqrt n\right) 种取值。
故，本质不同的 ⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 只有至多 O(n)\mathcal O\left(\sqrt n\right) 种取值。
因此 ⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 相同值肯定是连续的，那么就可以找出这 O(n)\mathcal O\left(\sqrt n\right) 区间 [l1,r1],[l2,r2],⋯ ,[lk,rk][l_1,r_1],[l_2,r_2],\cdots,[l_k,r_k]（ri+1=li+1r_i+1=l_{i+1}），求出：
g(⌊nl⌋)∑j=lirif(j)g\left(\left\lfloor\dfrac nl\right\rfloor\right)\sum_{j=l_i}^{r_i}f(j)
将 ii 个答案相加即可。
若分块后可以 O(1)\mathcal O(1) 求出，则数论分块是 O(n)\mathcal O\left(\sqrt n\right) 的。

寻找区间

对于一个 ll 所对应的 ⌊nl⌋\left\lfloor\dfrac nl\right\rfloor，要寻找一个最大的 rr 使得 ⌊nl⌋=⌊nr⌋\left\lfloor\dfrac nl\right\rfloor=\left\lfloor\dfrac nr\right\rfloor，从而找到区间 [l,r][l,r]。
则有：
r=⌊n⌊nl⌋⌋r=\left\lfloor\dfrac{n}{\left\lfloor\dfrac nl\right\rfloor}\right\rfloor
而对于向上取整的情形，即对于 ll 对应的 ⌈nl⌉\left\lceil\dfrac nl\right\rceil，要寻找一个最大的 rr 使得 ⌈nr⌉\left\lceil\dfrac nr\right\rceil，从而找到区间 [l,r][l,r]。
则有：
r=⌊n−1⌊n−1l⌋⌋r=\left\lfloor\dfrac{n-1}{\left\lfloor\dfrac {n-1}l\right\rfloor}\right\rfloor

单一参数

[CQOI2007] 余数求和
给定 n,k∈N∗n,k\in\N^*，满足 1≤n,k≤1091\leq n,k\leq10^9，求：
G(n,k)=∑i=1nk mod iG(n,k)=\sum_{i=1}^nk\bmod i
显然，k mod i=k−i⌊ki⌋k\bmod i=k-i\left\lfloor\dfrac ki\right\rfloor。
因此有：
G(n,k)=∑i=1n(k−i⌊ki⌋)=nk−∑i=1ni⌊ki⌋\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}
考虑快速求出 ∑i=1ni⌊ki⌋\sum\limits_{i=1}^ni\left\lfloor\dfrac ki\right\rfloor，可以进行数论分块。
对于数论分块中一组确定的 l,rl,r，有：
∑i=lri⌊ki⌋=⌊kl⌋∑i=lri=⌊ki⌋(l+r)(r−l+1)2\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}
这可以 O(1)\mathcal O(1) 计算，因此总计算时间复杂度为 O(n)\mathcal O\left(\sqrt n\right)。
参考代码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;
}

多参数

此时寻找 rr，会取 min⁡\min，即：
r=min⁡(⌊n⌊nl⌋⌋,⌊m⌊ml⌋⌋,⋯ )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] 模积和
给定 1≤n,m≤1091\leq n,m\leq10^9，求：
(∑i=1n∑j=1m[i≠j](n mod i)(m mod j)) mod 19940417\left(\sum_{i=1}^n\sum_{j=1}^m[i\neq j](n\bmod i)(m\bmod j)\right)\bmod 19940417
不妨令 n≤mn\leq m，否则交换 n,mn,m。
模运算不好处理，转换一下：
∑i=1n∑j=1m[i≠j](n−i⌊ni⌋)(m−j⌊mj⌋)\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≠j][i\neq j] 不好处理，因此可以先全部求出来，再减去相等的部分：
∑i=1n∑j=1m(n−i⌊ni⌋)(m−j⌊mj⌋)−∑i=1n(n−i⌊ni⌋)(m−i⌊mi⌋)\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)
不难发现 jj 与 ii 的取值无关，因此转换为：
∑i=1n(n−i⌊ni⌋)∑j=1m(m−j⌊mj⌋)−∑i=1n(n−i⌊ni⌋)(m−i⌊mi⌋)\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)
于是可以 O(max⁡(n,m))\mathcal O\left(\max(n,m)\right) 计算，但是考虑到 1≤n,m≤1091\leq n,m\leq10^9，只需要使用数论分块优化即可。
参考代码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(n) 满足 f(1)=1f(1)=1，且对于 ∀x,y∈N∗∧x⊥y\forall x,y\in\N^*\land x\perp y，有 f(xy)=f(x)⋅f(y)f(xy)=f(x)\cdot f(y)，称 f(n)f(n) 为积性函数。
若数论函数 f(n)f(n) 满足 f(1)=1f(1)=1，且对于 ∀x,y∈N∗\forall x,y\in\N^*，有 f(xy)=f(x)⋅f(y)f(xy)=f(x)\cdot f(y)，称 f(n)f(n) 为完全积性函数。

基础部分

若 f(n),g(n)f(n),g(n) 均为积性函数，则 f(xp),fp(x),f(x)g(x),(f∗g)(x)f(x^p),f^p(x),f(x)g(x),(f*g)(x) 均为积性函数。其中 f∗gf*g 为迪利克雷卷积：
(f∗g)(x)=∑d∣xf(d)g(xd)(f*g)(x)=\sum_{d\mid x}f(d)g\left(\dfrac xd\right)

常见积性函数

  • 单位函数：ε(n)=[n=1]\varepsilon(n)=[n=1]，完全积性函数。（\varepsilon）
  • 恒等函数：idk(n)=nk\mathrm{id}_{k}(n)=n^k，id1(n)\mathrm{id}_1(n) 通常记为 id(n)\mathrm{id}(n)，完全积性函数。
  • 常数函数：1(n)=11(n)=1，完全积性函数。
  • 除数函数：σk(n)=∑d∣ndk\sigma_k(n)=\sum\limits_{d\mid n}d^k，σ0(n)\sigma_0(n) 通常记作 d(n)d(n) 或 τ(n)\tau(n)，σ1(n)\sigma_1(n) 通常记作 σ(n)\sigma(n)。（\sigma）
    d(n)d(n) 的求法
    令 n=p1c1p2c2⋯pkckn=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}，其中 p1,p2,⋯ ,pkp_1,p_2,\cdots,p_k 均为质数。
    对于每一个指数 cic_i 的更改，都会产生新的因数。每一个 ii，都有 ci+1c_i+1 种选择。显然有：
    d(n)=∏i=1k(ci+1)d(n)=\prod_{i=1}^k(c_i+1)
    时间复杂度 O(n)\mathcal O\left(\sqrt n\right)。
  • 欧拉函数：φ(n)\varphi(n)。（\varphi）
  • 莫比乌斯函数：
    μ(n)={1n=10∃d>1,d2∣n(−1)kk 为 n 的本质不同质因子个数\mu(n)= \begin{cases} 1&n=1\\ 0&\exist d>1,d^2\mid n\\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数} \end{cases}
    μ(n)=0\mu(n)=0 即 nn 含有平方因子。（\mu）
同时，除数函数 d(n)d(n) 即 nn 的因数个数，且 d(ij)d(ij) 满足：
d(ij)=∑x∣i∑y∣j[x⊥y]\begin{aligned} d(ij)&=\sum_{x\mid i}\sum_{y\mid j}[x\perp y] \end{aligned}
证明
令 i=p1c1p2c2⋯pkck,j=p1′c1′p2′c2′⋯pk′′ck′′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'}}。如果将指数对位相加，则可以得到 ijij 的质因数分解形式。
将每一个质因子 plp_l 枚举 cl+cl′+1c_l+c'_l+1 次，便可以得到全部的因数。因此枚举 x∣i,y∣jx\mid i,y\mid j，但是需要保证 x,yx,y 不包含重复质因子，即 x⊥yx\perp y。

迪利克雷卷积

对于数论函数 f,gf,g，记 h=f∗gh=f*g 为 f,gf,g 通过迪利克雷卷积相乘的到的结果，即：
h(n)=∑d∣nf(d)g(nd)h(n)=\sum_{d\mid n}f(d)g\left(\dfrac nd\right)

常见迪利克雷卷积

  • φ∗1=id\varphi*1=\mathrm{id}。
  • μ∗1=ε\mu*1=\varepsilon。
  • μ∗1=φ\mu*1=\varphi。

莫比乌斯函数

莫比乌斯函数 μ(n)\mu(n) 满足：
∑d∣nμ(d)=[n=1]\sum_{d\mid n}\mu(d)=[n=1]
即 μ∗1=ε\mu*1=\varepsilon。
因此，μ(n)\mu(n) 可用于处理互质相关信息：
[i⊥j]=[gcd⁡(i,j)=1]=∑d∣gcd⁡(i,j)μ(d)=∑d∣i∧d∣jμ(d)[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)=∑d∣nf(d)g(n)=\sum\limits_{d\mid n}f(d)，则有：
f(n)=∑d∣nμ(d)g(nd)f(n)=\sum_{d\mid n}\mu(d)g\left(\dfrac nd\right)
因为：
∑d∣nμ(d)g(nd)=∑d∣nμ(d)∑d′∣ndf(d′)=∑d∣nμ(d)∑d′∣ndf(d′)=∑d∣nf(d)∑d′∣ndμ(d′)=∑d∣nf(d)[n=d]=f(n)\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)=∑n∣df(d)g(n)=\sum\limits_{n|d}f(d)，则有：
f(n)=∑n∣dμ(dn)g(d)f(n)=\sum_{n\mid d}\mu\left(\dfrac dn\right)g(d)
[国家集训队] Crash的数字表格 / JZPTAB
详见此处。
莫比乌斯反演的本质其实就是高维差分。

线性筛

线性筛可以用于求解欧拉函数及莫比乌斯函数。实际上，线性筛几乎可以求解所有的积性函数，求解方法视具体函数而定。

欧拉函数

若 ii 为质数，显然 φ(i)=i−1\varphi(i)=i-1。
设质数 primej\textit{prime}_j。
若 i mod primej≠0i\bmod\textit{prime}_j\neq0，则 i⊥primeji\perp\textit{prime}_j，因为 φ\varphi 是一个积性函数，于是有：
φ(i⋅primej)=φ(i)⋅φ(primej)\varphi\left(i\cdot\textit{prime}_j\right)=\varphi(i)\cdot\varphi\left(\textit{prime}_j\right)
否则，当 i≡0(modprimej)i\equiv0\pmod{\textit{prime}_j} 时，有：
φ(i⋅primej)=φ(i)⋅primej\varphi\left(i\cdot\textit{prime}_j\right)=\varphi(i)\cdot\textit{prime}_j
因为与 i⋅primeji\cdot\textit{prime}_j 互质，即与 ii 互质。1,2,3,⋯ ,i⋅primej1,2,3,\cdots,i\cdot\textit{prime}_j 在模 ii 意义下以 ii 为周期循环了 primej\textit{prime}_j 次。

莫比乌斯函数

若 ii 为质数，显然 μ(i)=−1\mu(i)=-1。
设质数 primej\textit{prime}_j。
若 i mod primej≠0i\bmod\textit{prime}_j\neq0 时，则 i⊥primeji\perp\textit{prime}_j，因为 μ\mu 是一个积性函数，于是有：
μ(i⋅primej)=μ(i)⋅μ(primej)=−μ(i)\mu\left(i\cdot\textit{prime}_j\right)=\mu(i)\cdot\mu\left(\textit{prime}_j\right)=-\mu(i)
否则，当 i≡0(modprimej)i\equiv0\pmod{\textit{prime}_j} 时，有：
μ(i⋅primej)=0\mu\left(i\cdot\textit{prime}_j\right)=0
因为 ii 含有 primej\textit{prime}_j，则 i⋅primeji\cdot\textit{prime}_j 至少含有两个 primej\textit{prime}_j，故 μ(i⋅primej)=0\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;
}

杜教筛

杜教筛用于快速求解积性函数 ff 的前缀和 S(n)=∑i=1nf(i)S(n)=\sum\limits_{i=1}^nf(i)。（实际上，只要能够构造恰当的函数，均可以使用杜教筛求解。）
构造一组 h=f∗gh=f*g，有：
∑i=1nh(i)=∑i=1n∑d∣if(id)g(d)=∑d=1ng(d)∑i=1⌊nd⌋f(i)=∑i=1ng(i)S(⌊ni⌋)\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}
可得：
g(1)S(n)=∑i=1ng(i)S(⌊ni⌋)−∑i=2ng(i)S(⌊ni⌋)=∑i=1nh(i)−∑i=2ng(i)S(⌊ni⌋)\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,gh,g 适当，能够快速地求出来，就能够在较短时间内求出 S(n)S(n)。

数论分块

可以发现，∑i=2ng(i)S(⌊ni⌋)\sum\limits_{i=2}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right) 直接计算的复杂度较高，而 ⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 就很适合用于数论分块来加速。

记忆化

每当我们求出了 S(n)S(n) 之后，都应当将其记录下来，避免重复计算。
这样的杜教筛的时间复杂度时 O(n34)\mathcal O\left(n^{\frac34}\right) 的。

线性筛预处理

可以预处理较小的一部分的 S(n)S(n)，从而节约时间。
当预处理前 O(n23)\mathcal O\left(n^\frac23\right) 的 S(n)S(n) 时，时间复杂度最小，为 O(n23)\mathcal O\left(n^\frac23\right)。

实际常数影响

在 OI 代码中，递归求解 S(n)S(n) 会有常数的影响。
因此最好的方法应当是实际测试预处理部分大小 mm 的值从而确定时间最小的 mm，这样即使复杂度劣一些，使用时间也小一些。
例题链接
参考代码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;
}

评论

共 12 条评论，欢迎与作者交流。

正在加载评论...