题解 数学 数论 提高+/省选− 省选/NOI− CRT exCRT exgcd BSGS exBSGS Lucas exLucas 未完

数论学习笔记

同余 CRT exgcd BSGS Lucas

Posted by TH911 on February 5, 2025

参考资料 参考资料

本文中出现的代数式、数，如无特殊说明，均为整数。

对于 OI 中的数论题，答案不对请先尝试：

1
2
3

#define int __int128
typedef long long ll;
#define ll __int128

数论函数指的是定义域为正整数的函数，也可以视为一个序列。

基础部分

$\perp$ 在数论中表示“互质”，例如 $2\perp3$。

质数

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$ 轮测试，通过二次探测定理判断。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

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\]

详见费马小定理。

欧拉定理

欧拉函数

记 $\varphi(n)$ 为 $1\sim n$ 中与 $n$ 互质的数的个数。

令 $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) $$

推导见容斥原理求解欧拉函数。

若 $\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}\)。

参考代码

1
2
3
4
5
6
7
8
9
10

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] 同余方程

求出 $x$ 后，为了使 $x$ 为正整数，只需要让 $x\leftarrow(x\bmod p+p)\bmod p$ 即可。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

//#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=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$。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

//#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\]

参考代码

1
2
3
4
5
6
7
8
9
10
11

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，防止溢出。

尽管题目大多会保证 $\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）/ 曹冲养猪

题目中的 $a_i,b_i$ 即分别为 $p_i,a_i$。

直接套 CRT 即可，但是需要注意的是，计算 $a_iq_ic_i$ 时需要 __int128，会爆 long long。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

//#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 解决的是当模数 $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)$ 是为了将两个方程合并。

这样，两两顺次合并，最终合并为一个方程即可解。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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';
    }
}
例题参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

//#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);
//	cin>>n;
	for(int i=1;i<=n;i++){
		IO::scanf(p[i]);
		IO::scanf(a[i]);
//		cin>>p[i]>>a[i];
	}
	//合并(i-1,i) 
	for(int i=2;i<=n;i++){
//		cerr<<"------------------\nmerge("<<i-1<<","<<i<<")\n";
		ll w=(a[i]-a[i-1])/gcd(p[i-1],-p[i]);
//		cerr<<"w="<<a[i]-a[i-1]<<" / "<<gcd(p[i-1],-p[i])<<endl;
		ll r,s;
		exgcd(p[i-1],r,-p[i],s);
		r*=w,s*=w;
//		cerr<<"r="<<r<<" s="<<s<<endl;
//		cerr<<"::A="<<p[i-1]*r+a[i-1]<<"\n::B="<<p[i]*s+a[i]<<endl;
		a[i]=(p[i-1]*r+a[i-1])%lcm(p[i-1],p[i]);
		p[i]=lcm(p[i-1],p[i]);
//		cerr<<"x≡"<<a[i]<<"(mod "<<p[i]<<")\n";
		if(i==n){
			if(a[n]<0){
				a[n]+=p[n];
			}
			IO::printf(a[n]);
//			cout<<a[n]<<'\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 的方法两两合并解决，具体见此。

离散对数

所谓离散对数，即模意义下取对数。

形如：给定模数 $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)$。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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

例题链接

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

//#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$ 最小的两个。

证明

由扩展欧拉定理，$a^x$ 是在 $k\leq\varphi(p)$ 值后以 $\varphi(p)$ 为周期循环的。

循环周期内的 $a^{qB}$ 显然至多存储 $1$ 个，而非循环部分 $a^{qB}$ 也至多存储 $1$ 个。

即：保留最小的 $2$ 个。

时间复杂度：仍然为 $\mathcal O\left(\sqrt{\varphi(p)}\right)$。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

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

例题链接

参考代码 
1
2
3
4
5	
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

//#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)$。

推荐使用乘法逆元。

例题链接

很简单，注意是 $\dbinom{n+m}{n}$ 即可。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

//#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!)$：

1
2
3
4
5
6
7
8

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)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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

例题链接

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

//#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
*/

[国家集训队] 礼物

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

//#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)$ 种。

证明
  • 当 $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] 余数求和

给定 $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)$。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

//#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] 模积和

给定 $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$，只需要使用数论分块优化即可。

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

//#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）

    $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) $$ </p>

    时间复杂度 $\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}\]
证明

令 $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

详见此处。

莫比乌斯反演的本质其实就是高维差分。

线性筛

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

欧拉函数

若 $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$。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

//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 ni\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$，这样即使复杂度劣一些，使用时间也小一些。

例题链接

参考代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

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

群论

群论是近代数学的基础之一。

群即对称，对称即群。

“群”是一种代数结构，群 $(G,\circ)$ 由一个集合 $G$ 和一个二元运算符 $\circ$ 组成。

群 $(G,\circ)$ 需要满足以下性质：

  • 封闭性。即对于 $\forall a,b\in G$，有 $a\circ b\in G$。

  • 结合律。即对于二元运算 $\circ$ 满足 $a\circ(b\circ c)=(a\circ b)\circ c$。

  • 存在单位元。即存在 $e\in G$ 满足对于 $\forall a\in G$ 有 $e\circ a=a\circ e=a$。

    例如在 $\circ$ 为乘法的情况下，单位元 $e$ 为 $1$。

  • 存在逆元。即对于 $\forall a\in G$，存在 $x\in G$ 满足 $a\circ x=x\circ a=e$。

    记 $x$ 为 $a^{-1}$。

同时，记：

\[\underbrace{a\circ a\circ a\circ \cdots\circ a}_{x\text{ 个 }a}=a^x\]

有时为了强调单位元，也会将群 $(G,\circ)$ 记作群 $(G,\circ,e)$。

例如这两个群：

  • 群 $(\Z,+,0)$。

    显然满足上文所述性质。

    无限群与整数集加法群同构。

  • 群 \((\lbrace0,1,2,\cdots,p-1\rbrace,+_{\bmod p})\)，其中 $+_{\bmod p}$ 表示模 $p$ 意义下的加法。

    显然也满足上文所述性质。

    有限群与模意义下的整数集加法群同构。

需要注意的是，群不一定需要满足交换律。

群的性质

已知 $(G,\circ,e)$ 是一个群。

单位元唯一

假设存在两个不同的单位元 $e_1,e_2\in G$。

那么有 $e_1=e_1\circ e_2=e_2$，则 $e_1=e_2$。

故，单位元唯一。

消去律

  • 左消去律：$a\circ b=a\circ c\Rightarrow b=c$。
  • 右消去律：$a\circ b=c\circ b\Rightarrow a=c$。

逆元唯一

假设 $a$ 存在两个不同的逆元 $a_1,a_2$。

则有 $a\circ a_1=a\circ a_2=e$。

由消去律有 $a_1=a_2$。

故，逆元唯一。

逆元的拆解

对于 $\forall a,b\in G$，有 $(a\circ b)^{-1}=b^{-1}\circ a^{-1}$。

由群的性质，可以得到：

\[(a\circ b)^{-1}\circ (a\circ b)=e\\ b^{-1}\circ a^{-1}\circ (a\circ b)=b^{-1}\circ(a^{-1}\circ a)\circ b=b^{-1}\circ e\circ b=e\]

因此，有 $(a\circ b)^{-1}=b^{-1}\circ a^{-1}$。

群的有关概念

已知 $(G,\circ,e)$ 是一个群。

有限群

当 $G$ 为有限集时。

有限群与模意义下的整数集加法群同构。

群的阶

群 $(G,\circ,e)$ 的阶即 $\vert G\vert$。

无限群

当 $G$ 为无限集时。

无限群与整数集加法群同构。

阿贝尔群

若二元运算 $\circ$ 满足交换律，即满足 $a\circ b=b\circ a$，则称群 $(G,\circ,e)$ 为阿贝尔群、可交换群或加法群。

子群

设集合 $H\in G$ 且 $\vert H\vert>0$。

如果满足 $(H,\circ,e)$ 为一个群，则称群 $(H,\circ,e)$ 是群 $(G,\circ,e)$ 的子群，记作 $H\leq G$。

平凡子群

对于任意群 $(G,\circ,e)$，总有平凡子群 $(G,\circ,e)$ 和 $(\lbrace e\rbrace,\circ,e)$。

CATALOG
    FEATURED TAGS
    树状数组 题解 提高+/省选− 省选/NOI− DP 普及+/提高 基础算法 数学 组合数学 普及/提高− 贪心 字符串 线段树 未完 图论 平衡树 数论 树形 DP 搜索 Tarjan Trie 二分 位运算 二分答案 状态压缩 FHQ Treap 前端 无聊时候的一点思考 最短路 树链剖分 NOI/NOI+/CTSC
    FRIENDS