阶

定义整数 $a$ 如果存在一个正整数 $n$ 满足 $a^n \equiv 1 \pmod m$，那么最小的那个 $n$ 称作整数 $a$ 的阶，记作 $\operatorname{ord}_m(a)$ 或者 $\delta_m(a)$。

容易发现，根据欧拉定理若 $\gcd(a,m)=1$ 那么 $a^{\varphi(m)} \equiv 1 \pmod m$。

性质 1

$a^1,a^2,a^3,\dots,a^{\delta_m(a)}$ 在模 $m$ 意义下两两不同余。

证明

若设 $a^i \equiv a^j \pmod m$，那么 $a_{j-i} \equiv 1 \pmod m$，因为 $1 \le i < j \le \delta_m(a)$，所以 $j-i < \delta_m(a)$，与阶的定义相违背。

性质 2

若 $a^n \equiv 1 \pmod m$ 那么 $\delta_m(a) \mid n$。

证明

若设 $n=p\delta_m(a)+q$，则 $0 \le q < \delta_m(a)$，并且 $a^{p \delta_m(a)} \equiv 1 \pmod m$，所以 $a^q \equiv 1 \pmod m$，如果 $q \ge 1$，与阶的定义相违背，所以 $q=0$，即 $\delta_m(a) \mid n$。

推论

若 $a^p \equiv a^q \pmod m$，那么 $p \equiv q \pmod {\delta_m(a)}$。

性质 3

若 $m \in N^*$，$a,b \in Z$，$\gcd(a,m)=\gcd(b,m)=1$，则：

$$
\delta_m(a)\delta_m(b)=\delta_m(ab)
$$

的充要条件是 $\gcd(\delta_m(a),\delta_m(b))=1$。

证明

必要性

因为 $a^{\delta_m(a)} \equiv 1$，$b^{\delta_m(b)} \equiv 1$，所以 $(ab)^{[\delta_m(a),\delta_m(b)]} \equiv 1$。

那么根据性质 $2$ 可得 $\delta_m(ab) \mid [\delta_m(a),\delta_m(b)]$，所以 $\delta_m(a)\delta_m(b) \mid [\delta_m(a),\delta_m(b)]$，推出 $(\delta_m(a),\delta_m(b))=1$。

充分性

由 $(ab)^{\delta_m(ab)} \equiv 1$，那么 $(ab)^{\delta_m(ab)\delta_m(b)} \equiv a^{\delta_m(ab)\delta_m(b)}\equiv 1$。

所以 $\delta_m(a)\mid \delta_m(ab)\delta_m(b)$，因为 $\gcd(\delta_m(a),\delta_m(b))=1$，所以 $\delta_m(a) \mid \delta_m(ab)$。

对称地，那么 $\delta_m(b) \mid \delta_m(ab)$。

所以 $\delta_m(a)\delta_m(b) \mid \delta_m(ab)$。

另一方面有：$(ab)^{\delta_m(a)\delta_m(b)}\equiv (a^{\delta_m(a)})^{\delta_m(b)} (b^{\delta_m(b)})^{\delta_m(a)} \equiv 1$。

所以 $\delta_m(ab) \mid \delta_m(a)\delta_m(b)$。

综上所述，有 $\delta_m(a)\delta_m(b)=\delta_m(ab)$。

性质 4

对于合法的 $a,k,m$，有 $\delta_m(a^k)=\dfrac{\delta_m(a)}{\gcd(\delta_m(a),k)}$。

证明

首先有 $a^{k\delta_m(a^k)} \equiv (a^k)^{\delta_m(a^k)} \equiv 1$，所以 $\delta_m(a) \mid k\delta_m(a^k)$，推得 $\dfrac{\delta_m(a)}{\gcd(\delta_m(a),k)} \mid \delta_m(a^k)$。

还有 $a^{\delta_m(a)} \equiv 1$，可得 $(a^k)^{\frac{\delta_m(a)}{\gcd(\delta_m(a),k)}} \equiv (a^{\delta_m(a)})^{\frac{k}{\gcd(k,\delta_m(a))}} \equiv 1$。

所以 $\delta_m(a^k) \mid \dfrac{\delta_m(a)}{\gcd(\delta_m(a),k)}$，综上所述，$\delta_m(a^k) = \dfrac{\delta_m(a)}{\gcd(\delta_m(a),k)}$。

阶的寻找

如果 $\gcd(a,m) \not=1$ 的话，说明 $a^k \not\equiv 1 \pmod m$，可以用 $a \equiv b \pmod m$ 与 $\dfrac{a}{\gcd(a,m)} \equiv \dfrac{b}{\gcd(a,m)} \pmod {\dfrac{m}{\gcd(a,m)}}$ 等价证明。

所以当 $\gcd(a,m)=1$ 时 $a$ 才有阶，此时 $a^{\varphi(m)} \equiv 1$，根据欧拉定理及其推论可以得 $\delta_m(a) \mid \varphi(m)$。

于是我们对于 $\varphi(m)$ 分解质因数，然后依次除以每个质因数，再用快速幂判断是不是为 $1$ 就可以了，如果为 $1$ 就除一下，否则继续试除即可。

寻找阶的代码（没有用试除优化）：

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

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define ll long long
using namespace std;
ll n,i,j,p,q,l[200005],tot;
inline ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
int main(){
	scanf("%lld",&n),q=n,p=n;
	for(i=2;i*i<=q;i++){
		if(q%i==0){
			p -= p/i;
			while(q%i==0) q/=i;
		}
	}
	if(q>1) p-=p/q;
	for(i=1;i<=p;i++) if(p%i==0) l[++tot]=i;
	for(i=1;i<n;i++){
		if(__gcd(i,n)!=1){
			printf("0 ");
			continue;
		}
		for(j=1;j<=tot;j++){
			if(qmi(i,l[j],n)==1){
				printf("%lld ",l[j]);
				break;
			}
		}
		if(j>tot) printf("0 ");
	}
	return 0;
}

原根

定义

对于合法的 $g,m$，如果其满足 $\delta_m(g) = \varphi(m)$，那么则称 $g$ 是模 $m$ 的原根。

容易发现，当 $m$ 为质数时对于任意 $0<i<m$ 的 $g_i \bmod m$ 互不相同。

判定

若 $g$ 是模 $m$ 的原根，首先满足 $\gcd(g,m)=1$（欧拉定理），然后对于 $m$ 的每个质因数 $t$ 都要满足 $g^{\frac{m}{t}} \not\equiv 1$。

证明显然，此处不展开叙述。

个数

如果一个数 $m$ 有原根，那么它的个数是 $\varphi(\varphi(m))$。

证明

设 $g$ 是 $m$ 的任意一个原根，则有：

$$
\delta_m(g^k)=\dfrac{\delta_m(g)}{\gcd(\delta_m(g),k)}
$$

所以当 $\gcd(\delta_m(g),k)=1$ 时，$\delta_m(g^k)=\varphi(m)$，因为 $\delta_m(g)=\varphi(m)$，所以个数是 $\varphi(\varphi(m))$。

原根存在定理

一个数 $m$ 存在原根当且仅当它等于 $2$ 或 $4$ 或 $p^k$ 或 $2p^k$。

证明略。

最小原根范围

可证 $g$ 大概在 $m^{\frac{1}{4}}$ 左右，所以暴力找原根的时间复杂度是可以接受的。

寻找原根的代码：

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

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define ll long long
#define N 1000005
using namespace std;
ll n,i,j,k,p,q,l[N],tot,ans[N],id,gcd[N],pri[N],ttt,vis[N];
inline ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
int main(){
	scanf("%lld",&n),q=n,p=n;
	for(i=2;i<=n;i++){
		if(!vis[i]){
			pri[++ttt]=i;
			if(n%i==0) gcd[i]=i;
			else gcd[i]=1;
		}
		for(j=1;j<=ttt;j++){
			if(pri[j]*i>n) break;
			vis[i*pri[j]]=1;
			gcd[i*pri[j]]=gcd[i]*gcd[pri[j]];
			if(i%pri[j]==0) break;
		}
	}
	for(i=2;i*i<=q;i++){
		if(q%i==0){
			p -= p/i;
			while(q%i==0) q/=i;
		}
	}
	if(q>1) p-=p/q;
	q=p;
	for(i=2;i*i<=q;i++){
		if(q%i==0){
			l[++tot] = p/i;
			while(q%i==0) q/=i;
		}
	}
	if(q>1) l[++tot] = p/q;
	for(i=1;i<n;i++){
		if(gcd[i]!=1) continue;
		for(j=1;j<=tot;j++) if(qmi(i,l[j],n)==1) break;
		if(j>tot){
			for(j=i,k=1;k<=p;j=j*i%n,k++) if(__gcd(k,p)==1) ans[++id]=j;
			break;
		}
	}
	sort(ans+1,ans+id+1);
	printf("%lld\n",id);
	for(i=1;i<=id;i++) printf("%lld ",ans[i]);
	return 0;
}

BSGS 算法

定义

已知 $a,b,p$，求最小的 $t$ 满足下列公式：

$$
a^t \equiv b \pmod{p}
$$

普通的 BSGS 算法能够在 $\gcd(a,p)$ 互质时以 $O(\sqrt{p})$ 的时间复杂度内得出答案。

解决

核心思想是根号分治。

设 $r = \lceil{\sqrt{t}}\rceil$，则 $t=rk+q(0 \leq q < r)$。

写到式子上就是：

$$
a^{rk+q} \equiv b \pmod{p}
$$

因为 $a$ 有模 $p$ 意义下的逆元，我们可以将式子转化为下面的形式：

$$
a^{rk} \equiv b \times a^{-q} \pmod{p}
$$

注意到左边有不超过根号种取值，右边也有不超过根号种取值，于是预处理两边中的任意一边即可，剩下的那边暴力枚举就好。

时间复杂度看你存储一边的时候使用的数据结构，如果是 $\text{map}$，时间复杂度多一个 $\log$；如果是 $\text{unordered_map}$，时间复杂度多一个常数。（模运算具有较大的随机性）

除此之外，我们也可以写成 $a^t=a^{kr-q}$ 的形式，代码就是采用的这种结构，但是边界情况一定要判断好。

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

#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll p,a,b,t,i,j,cnt,ans=-1;
map<ll,ll> op;
ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t = a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
int main(){
    cin>>p>>a>>b;
    t = ceil(sqrt(p*1.0));
    for(j=0;j<t;j++) op[b*qmi(a,j,p)%p] = j;
    a = qmi(a,t,p);
	for(i=0;i<=t;i++){
		cnt = qmi(a,i,p);
		if(op.find(cnt)!=op.end()) ans=op[cnt];
		else ans=-1;
		if(ans>=0&&i*t-ans>=0){
			cout<<i*t-ans<<endl;
			return 0;
		}
	}
	cout<<"no solution";
}

广义 BSGS

如果定义一种运算 $p \times q = k$，并且对于 $a$ 存在某个结构 $a^{-1}$ 使得 $a \times a^{-1}=e$，$e$ 是这种运算的单位元，那么都可以采用 BSGS 求解。

比如矩阵就可以采用 BSGS 求解，即把上述运算全部改成矩阵的乘法、快速幂即可，时间复杂度为 $O(\sqrt{p}) \times k$，$k$ 是执行一次我们定义的运算的时间复杂度。

Millar-Rabin

判定素数经常被认为是一个 NP 问题，但是实际上有确定性多项式算法解决这个问题，但是它在算法竞赛中并不常用，而 Millar-Rabin 是一个非确定性的算法来解决素数判定问题。

时间复杂度为执行判定的次数乘上 $\log n$。

理论

有费马小定理当 $p$ 是质数的时候 $a^{p-1} \equiv 1 \pmod p(1 \le a < p)$，但是反过来是否如此呢？

事实上反过来会有一种数满足它仍然成立但是不是质数（$a$ 不整除它），这种数称作 Carmichael 数，例如 $561=3 \times 11 \times 17$ 就是 Carmichael 数。

Carmichael 数一定有超过三个不同的质因子，并且不包含平方因子。

于是根据这点，我们得到了一个最简单的判定方法，称作 Fermat 素性测试，代码如下：

1
2
3
4
5
6
7
8
9
10

bool millerRabin(int n) {
  if (n < 3) return n == 2;
  // test_time 为测试次数,建议设为不小于 8
  // 的整数以保证正确率,但也不宜过大,否则会影响效率
  for (int i = 1; i <= test_time; ++i) {
    int a = rand() % (n - 2) + 2;
    if (quickPow(a, n - 1, n) != 1) return 0;
  }
  return 1;
}

但是我们这个测试无法确定 Carmichael 数的素性，于是我们引进了二次探测定理。

若 $p$ 是质数，则满足 $x^2 \equiv 1 \pmod p$ 的 $x$ 只可能等于 $1$ 或 $p-1$。

我们可以分解一下 $(x-1)(x+1) \equiv 0 \pmod p$，因为 $p$ 是质数，所以 $x-1$ 或者 $x+1$ 等于 $p$。

那么也有数可以通过二次探测定理的逆命题，但是这些数一定是 $p^k$，与 Carmichael 没有交集，故我们可以根据这两个定理设计素数判定的方案。

处理

将 $a^{n-1} \equiv 1 \pmod n$ 中的指数 $n-1$ 分解为 $n-1=u \times 2^t$，在每轮测试中对随机出来的 $a$ 先求出 $v = a^{u} \bmod n$，之后对这个值执行最多 $t$ 次平方操作，若发现非平凡平方根时即可判断出其不是素数，否则再使用 Fermat 素性测试判断。

有一个问题，就是如果 $v$ 在进入循环之前就等于 $1$ 了，那就直接跳过即可，如果 $t$ 次平方操作之后结果仍然不为 $1$（数字不为 $n-1$），那就返回不是素数即可。

代码如下：

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

inline bool check(ll x){
	if(x==2) return 1;
	if(x<2||x%2==0) return 0;
	ll u = n-1,now = 1,temp = 0;
	while(u%2==0) now*=2,u/=2,temp++;
	for(ll i=0;i<12;i++){
		ll cnt = (rnd()%1000000000+1)*(rnd()%1000000000+1)%(x-2)+2;
		cnt = qmi(cnt,u,x);
		if(cnt==1) continue;
		for(ll j=1;j<=temp;j++){
			if(cnt==x-1) break;
			cnt = (__int128)cnt*cnt%x;
		}
		if(cnt!=x-1) return 0;
	}
	return 1;
}

选择的数字看判断素数的范围，如果是 int 类型，那么取判定的底数 $a=2,7,61$ 即可；否则取 $a$ 为前 $12$ 个质数即可。

Pollard-Rho

给定一个正整数 $N \in \mathbf{N}_{+}$，试快速找到它的一个非平凡因数。

首先，有一个悖论称作“生日悖论”，如果一年有 $365$ 天，那么一个班如果有至少 $23$ 个人，那么有两个人生日相同的概率将达到 $50%$，如果有至少 $56$ 个人，概率近乎 $100%$。

我们还可以通过最大公因数来找出某个数的非平凡因子，即如果 $\gcd(k,n)>1$，那么 $\gcd(k,n)$ 一定是大于 $1$ 的 $n$ 的因子。

我们考虑一个启发式伪随机 $f(x)=(x^2+c)\bmod n$ 来生成一个序列 $x_i$：随机取一个 $x_1$，令 $x_2=f(x_1),x_3=f(x_2),\cdots,x_i=f(x_{i-1})$。其中 $c$ 是一个随机选取的常数。

因为这个序列生成的数列是一个 rho 形态的，就是一个环上面掉了一个链，所以成为 Pollard-Rho。

这个函数生成序列还有一个性质，如果 $x \equiv y \pmod p$，那么 $f(x) \equiv f(y) \pmod p$。

证明显然。

那么根据生日悖论，我们可以在这个环上任取两个数使得这两个数的差的 $n$ 的最大公因数大于 $1$，我们就可以找出来一个因数，这样的时间复杂度就是 $O(n^{\frac{1}{4}})$。

我们就可以设置两个指针，一个走慢一点，一个走快一点，大概是一倍的速度差，然后判断这两个指针的差值即可。

同时因为 $\gcd(a,n)>1$，那么 $\gcd(ac,n)>1$，我们可以通过倍增计算，即每隔 $2^k-1$ 个数，计算两个指针的差值，如果两个指针相等了，那么说明它们处在同一个位置，重新执行即可。

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

inline ll pollard_rho(ll x){
	ll s=0,t=0,c=(rnd()%1000000000+1)*(rnd()%1000000000+1)%(x-1)+1;
	for(ll i=1,val=1;;i*=2,s=t,val=1){
		for(ll j=1;j<=i;j++){
			t=((__int128)t*t+c)%x;
			val=(__int128)val*(ll)(abs(s-t))%x;
			if(!val) return x;
			if(j%127==0){
				ll temp = __gcd(val,x);
				if(temp>1) return temp;
			}
		}
		ll temp = __gcd(val,x);
		if(temp>1) return temp;
	}
}
inline void found(ll x){
	if(check(x)){
		d[++tot] = x;
		return ;
	}
	ll temp = x;
	while(temp>=x) temp=pollard_rho(x);
	found(x/temp),found(temp);
}

总的时间复杂度是小常 $\log$ 乘上 $n^{\frac{1}{4}}$。

MIN_25 筛

基础算法

给定 $n$，求 $\pi(n)$，$n \le 10^{11}$。

我们考虑埃氏筛法，那么每一轮枚举 $1 \sim n$ 的质数，但是事实上只需要枚举 $1 \le \sqrt{n}$ 的质数就可以了，因为一个数是合数的话其最小质因子 $\le \sqrt{n}$。

所以对于每一轮枚举完之后剩下没有被标记的数有 $g(n,k)$ 个，$n$ 是值域，$k$ 是枚举完的轮数，设 $p_k$ 是从小到大第 $k$ 个质数，则有下面的公式：

$$
g(n,k)=g(n,k-1)-g(\lfloor \dfrac{n}{p_k} \rfloor,k-1)+k-1
$$

考虑意义，首先我们考虑 $k$ 这一轮筛掉了多少数字，设它们分别为 $a_1,a_2,a_3,\dots,a_m$，那么它们的最小质因子等于 $p_k$，则让它们除以 $p_k$ 之后得到 $\dfrac{a_1}{p_k},\dfrac{a_2}{p_k},\dfrac{a_3}{p_k},\dots,\dfrac{a_m}{p_k}$，这个时候它们的最小质因子大于等于 $p_k$，则这些数的个数是 $g(\lfloor \dfrac{n}{p_k} \rfloor,k-1)-(k-1)$，为什么要减 $k-1$ 呢，是因为这一坨同时算上了质数 $2,3,5,7,\dots$，为了去除多算的，就要减去 $k$ 之前的质数个数 $k-1$ 个。

代码如下：

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

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define ll long long
#define N 1000005
using namespace std;
ll n,K,i,j,r,pri[N],vis[N],tot,f1[N],f2[N],used[N],used_inv[N],utot;
inline ll found(ll x){
	if(x<=K) return f1[x];
	else return f2[n/x];
}
int main(){
	ios::sync_with_stdio(false);
	cin>>n;
	K=sqrt(n);
	for(i=2;i<=K;i++){
		if(!vis[i]) pri[++tot] = i;
		for(j=1;j<=tot;j++){
			if(i*pri[j]>K) break;
			vis[i*pri[j]] = 1;
			if(i%pri[j]==0) break;
		}
	}
	for(i=1;i<=n;){
		r = n/(n/i);
		if(n/i>1) used[++utot] = n/i;
		i = r+1;
	}
	reverse(used+1,used+utot+1);
	for(j=1;j<=utot;j++){
		used_inv[j] = n/used[j];       //过多用除法会导致常数爆炸
		if(used[j]<=K) f1[used[j]]=used[j]-1;
		else f2[used_inv[j]]=used[j]-1;
	}
	for(i=1,r=1;i<=tot;i++){
		for(j=utot;j>=1;j--){
			if(used[j]<pri[i]*pri[i]) break;
			if(used[j]<=K) f1[used[j]]=f1[used[j]]-found(used[j]/pri[i])+i-1;
			else f2[used_inv[j]]=f2[used_inv[j]]-found(used[j]/pri[i])+i-1;
		}
	}
	cout<<f2[1]<<endl;
	return 0;
}
/*
Input:
10000000000

Output:
455052511
f(pk,n)=f(pk-1,n)-f(pk-1,n/pk)+k-1
*/

时间复杂度证明此处就不说了，但是总的时间复杂度大概是 $O(\dfrac{n^{\frac{3}{4}}}{\ln n})$。

特别的，如果是求 $1 \sim n$ 中质数的和，把公式改为下面即可（其实就是动态规划）：

$$
f(n,k)=f(n,k-1)-p_k \times f(\lfloor \dfrac{n}{p_k} \rfloor,k-1)+\sum_{i=1}^{k-1} p_i \times p_k
$$