线性筛法及扩展

发表于 | 更新于 | 分类于 | 评论数： | 阅读次数：
本文字数： 6k | 阅读时长 ≈ 0:06

以下代码中的定义：

  • mindiv[i]: i的最小质因子
  • phi[i]: 欧拉函数i的值
  • mindivq[i]: i的最小质因子的个数
  • d[i]: i的约数个数
  • sumd[i]: i的约数和
  • miu[i]: 莫比乌斯函数i的值
  • inv[i]: i在mod n意义下的乘法逆元

标准筛法

欧拉筛法，可以保证每个数只被自己最小的质因子筛去，时间复杂度O(n)
两种等价实现
易于理解版：

1
2
3
4
5
6
7
8
9
10
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i])
			mindiv[i]=pris[++cnt]=i;
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0)break;
		}
	}
}

又短又快版（避免了取模运算）：（待更新）

1
2
3
4
5
6
7
8
9
10
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i])
			mindiv[i]=pris[++cnt]=i;
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			if(i%pris[j]==0)break;
			mindiv[k]=pris[j];
		}
	}
}

扩展：求积性函数

积性函数定义：对于一个数论函数 f(x)f(x)，满足：若 gcd(a,b)=1gcd(a,b)=1 ，则 f(ab)=f(a)f(b)f(ab)=f(a)f(b) 的函数 ff 。
由于这个性质，我们可以在筛法的同时求出积性函数的值
先看几个简单的积性函数：

id函数

定义 id(i)=iid(i)=i ，显然满足积性。
这就不用筛了，直接得到值

e函数

又称单位函数，定义

e(n)={1,if n=10,otherwisee(n)=\begin{cases} 1, & \text{if $ n=1 $}\\ 0, & \text{otherwise} \end{cases}

phi函数

欧拉函数，很有用，可以用来求逆元。
定义

ϕ(i)=∑j=1ie(gcd(j,i))\phi(i)=\sum_{j=1}^ie(gcd(j,i))

即 [1,n][1,n] 中，与 nn 互质的数的数量
易得对于质数，phi(p)=p−1phi(p)=p-1，
对于合数有以下公式：

n=∏i=1rpiqin=\prod_{i=1}^r p_i^{q_i}

ϕ(n)=n∗(1−1/p1)∗(1−1/p2)∗...\phi(n)=n*(1-1/p_1)*(1-1/p_2)*...

将n带入可得

phi(n)=∏i=1r(pi−1)∗piqi−1phi(n)=\prod_{i=1}^r(p_i-1)*p_i^{q_i-1}

这就很容易发现它的积性性质
用筛法怎么求呢？
稍微改一下就行了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i]){
			mindiv[i]=pris[++cnt]=i;
			phi[i]=i-1;
		}
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0){
				phi[k]=phi[i]*pris[j];
				break;
			}
			phi[k]=phi[i]*(pris[j]-1);
		}
	}
}

为什么呢？

  1. 对于 i%pris[j]!=0 的项，它的最小质因子一定是 pris[j] ，则由积性函数性质可得 phi[k]=phi[i]*phi[pris[j]] ，其中phi[pris[j]]=pris[j]-1
  2. 对于 i%pris[j]==0 的项，则说明 pris[j] 已经在 i 中已经出现了，而且是最小的，观察 ϕ\phi 函数的计算式，可知此时应乘 pris[j]

最小质因数的指数

这个就很简单了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i]){
			mindiv[i]=pris[++cnt]=i;
			mindivq[i]=1;
		}
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0){
				mindivq[k]=mindivq[i]+1;
				break;
			}
			mindivq[k]=1;
		}
	}
}

约数个数

首先约数个数如何求？

分解质因数，得到

n=∏i=1rpiqin=\prod_{i=1}^r p_i^{q_i}

则其 约数个数

=∏i=1r(qi+1)=\prod_{i=1}^r(qi+1)

这个结论由乘法原理易得，每个质因子有 qi+1qi+1 种选法。
积性就不证了，比较显然。

怎么求呢？

观察筛法的过程， i 是质数时或 i%pris[j]!=0 时都比较容易，主要是 i%pris[j]==0 时需要考虑好。

下面来分析一下：

关键在于如何由 d(i) 转移到 d(i*pris[j])
首先筛法保证了 pris[j] 一定是最小质因子，那么由于 i%pris[j]==0 ，则意味着最小质因子的指数会 +1+1，这一点在上面的 mindivq 求解过程中也体现了。

+1+1 会导致什么？

很自然会想到

d(k)=d(i)∗(mindivq[k]=mindivq[i]+1)+1mindivq[i]+1d(k)=d(i)*\frac{(mindivq[k]=mindivq[i]+1)+1}{mindivq[i]+1}

就是这样。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i]){
			mindiv[i]=pris[++cnt]=i;
			mindivq[i]=1;
			d[i]=2;
		}
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0){
				mindivq[k]=mindivq[i]+1;
				d[k]=d[i]/(mindivq[i]+1)*(mindivq[k]+1);
				break;
			}
			mindivq[k]=1;
			d[k]=d[i]*d[pris[j]];
		}
	}
}

约数和

约数和又如何求？
分解质因数，得到

n=∏i=1rpiqin=\prod_{i=1}^r p_i^{q_i}

则其 约数和

=∏i=1r∑j=0qipij=\prod_{i=1}^r\sum_{j=0}^{q_i} p_i^j

看起来也是比较显然，我们直接讨论如何利用筛法来求：

依然是利用和上面的类似思路：
只讨论 i%pris[j]==0 的情况：
需要除以 ∑n=0mindivq[i](pris[j])n\sum_{n=0}^{mindivq[i]}(pris[j])^n
再乘以 ∑n=0mindivq[k](pris[j])n\sum_{n=0}^{mindivq[k]}(pris[j])^n
其中 mindivq[k]=mindiv[i]+1mindivq[k]=mindiv[i]+1
这样开两个辅助数组记录

t1[i]=∑n=0mindivq[i]pris[j]nt1[i]=\sum_{n=0}^{mindivq[i]}pris[j]^n

t2[i]=(mindiv[i])mindivq[i]t2[i]=(mindiv[i])^{mindivq[i]}

就可以做了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i]){
			mindiv[i]=pris[++cnt]=i;
			sumd[i]=1+i;
			t1[i]=1+i;
			t2[i]=i;
		}
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0){
				t2[k]=t2[i]*pris[j];
				t1[k]=t1[i]+t2[k];
				sumd[k]=sumd[i]/t1[i]*t1[k];
				break;
			}
			sumd[k]=sumd[i]*sumd[pris[j]];
			t1[k]=1+pris[j];
			t2[k]=pris[j];
		}
	}
}

miu函数

莫比乌斯函数，定义

μ(n)={1,if n=10,if  is square free numbern(−1)质因数个数otherwise\mu(n)=\begin{cases} 1, & \text{if $ n=1 $}\\ 0, & \text{if $n$ is square free number}\\ (-1)^{\text{质因数个数}}& \text{otherwise} \end{cases}

用筛法求也很简单：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void solve(int n){
	for(int i=2;i<=n;i++){
		if(!mindiv[i]){
			mindiv[i]=pris[++cnt]=i;
			miu[i]=-1;
		}
		for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
			mindiv[k]=pris[j];
			if(i%pris[j]==0){
				miu[k]=0;
				break;
			}
			miu[k]=-miu[i];
		}
	}
}

乘法逆元

由求逆元的欧拉定理方法易知乘法逆元也是积性的，而且有一个很优美的性质：它是完全积性的。

所以求它变得很简单：

inv[p]=pphi(n)−1modninv[p]=p^{phi(n)-1}\mod n

inv[k]=inv[i]∗inv[pris[j]]if gcd(k,n)=1inv[k]=inv[i]*inv[pris[j]] \quad\text{if } gcd(k,n)=1

然而乘法逆元有一个递推的方法，更加简单：

首先逆元满足以下等式（定义）：

inv[P mod i]∗(P mod i)=1inv[P\text{ mod }i]*(P\text{ mod }i)=1

将P % i展开：

inv[P mod i]∗(P−⌊P/i⌋∗i)=1inv[P\text{ mod }i]*(P-\lfloor P/i \rfloor *i)=1

进一步展开：

inv[P mod i]∗i∗(−⌊P/i⌋)=1inv[P\text{ mod }i]*i*(-\lfloor P/i \rfloor)=1

会发现一个神奇的等式：

i∗(−⌊P/i⌋)∗(inv[P mod i])=1i*(-\lfloor P/i \rfloor)*(inv[P\text{ mod }i])=1

此时 (−⌊P/i⌋)∗(inv[P mod i])(-\lfloor P/i \rfloor)*(inv[P\text{ mod }i]) 就是 ii 的逆元！

总结

线性欧拉筛法简洁，积性函数性质优美，充分结合才能发挥更大的效果。同时又要注意到不同函数的性质，用筛法求积性函数的本质是根据积性函数的定义和运算式，利用筛法的特点，构造合适的转移方法。

0%