筛选素数:埃氏筛选、欧拉筛选(线性筛选)、区间筛选

一、知识储备

基本的素数知识

  • 素数除了1和它本身以外没有别的因数
  • 任何一个大于1的自然数 N,如果N不为质数,那么N可以唯一分解成有限个质数的乘积
  • 一个数c=a*b,那么要么a=b,要么a和b一个大于根号c,一个小于根号c

二、朴素筛选

先讲一下最基本的朴素筛选
判断n是不是质数,最基本的方式i从2~n-1枚举,判断n是否整除i,若对与每一个i,n都不能整除,则说明n是素数。
很明显2以上的偶数都不是素数,又a*b=c,则总有一个因素小于根号c,故可以缩小枚举范围,具体见下方代码。

const int MAXN = 1e5;
bool primeArr[MAXN]; //表示primeArr[i]=True,表示i为素数	
//判断k是不是素数
bool Judge(int k)
{
	if (k == 2)return true;
	if (k < 2 || k%2==0)return false;
	for (int i = 3; i <= sqrt(k); i += 2)
	{
		if (k%i == 0)return false;
	}
	return true;
}
//朴素筛选
void PurifiedSieve(int n)
{
	for (int i = 0; i <= n; i++)
	{
		primeArr[i] = Judge(i);
	}
}

三、埃氏筛选

埃氏筛选原理:
给出要筛数值的范围n,找出sqrt(n)以内的素数 p_{1},p_{2},…,p_{k}。先用2去筛,即把2留下,把2的倍数剔除掉;
再用下一个素数,也就是3筛,把3留下,把3的倍数剔除掉,考虑4,4被2筛除了;
接下去用下一个素数5筛,把5留下,把5的倍数剔除掉;不断重复下去直至根号n
可以得到小于等于n的所有素数,它的复杂度是O(n log log n)
为什么是根号n呢?
以100为例,sqrt(102)=10,
依次筛选2、3、(4为2的倍数,已被筛选)5、(6同理)、7、(8)(9)(10)的倍数
接下来考虑11,11的2倍被2筛选过了,同理3倍被3筛选,4倍被2筛选···
到10倍的时候已经超过原理的102了
下方代码注释加深理解

const int MAXN = 1e5;
bool primeArr[MAXN];

void EraSieve(int n)
{
	for (int i = 0; i < n; i++)
		primeArr[i] = true; //初始化每个数都为素数
	primeArr[0] = primeArr[1] = false;	//初始化0和1不是素数
	for (int i = 2; i <= sqrt(n); i++)	//从i从2到根号n作为质数,依次剔除它的2、3、4……倍
		if (primeArr[i])	//若这个数是素数则剔除它的整数倍
			for (int j = i*2; j <= n; j += i)	//剔除整数倍,上限是n(这里会有重复剔除的数,欧拉筛可优化)
				primeArr[j] = false;
}

四、欧拉筛(线性筛)

在埃氏中也有提到,一个数在埃氏筛中可能会被多次筛选考虑,多余的考虑就增加了算法的复杂度,而欧拉筛选算法可优化至线性时间复杂度。
先看一下原理:

  • 设minP[i]为i最小质因子。
  • i从从2开始枚举筛选
  • 若minP[i]没有在之前得到它的最小质因子,就说明i是质数,记录下i为质数,即pr[i]=i,prime[++len]=i
  • 对于i,枚举每一个不超过minP[i]的质数prime[j],并将i的prime[j]倍标记为非质数,即pr[i∗prime[j]]=prime[j]记录下了pr[i∗prime[j]]的最小质因子
    这样可以发现,每一个数只会被它最小的质因子筛去(下方举例说明),所以保证了算法为O(n)的。

原理看不懂没关系,可以抛开原理,结合下方例子以及代码进行理解:

值得注意的是:
埃氏筛是在一个循环中依次剔除素数1,素数2,素数3…的2倍,3倍,4倍……
而欧拉筛是在一个循环中依次剔除2,3,4…n的素数1倍,素数2倍,素数3倍(枚举已确定的素数)

int minP[MAXN];//minP[i] 为i的最小素因子 
int prime[MAXN];//prime[i] 代表n以内从小到大,第i个质数 
int cnt;//素数计数 

void EulerSieve(int n)
{
	for (int i = 2; i <= n; i++)
	{
		//从2开始,若minP[i]不为0,则表示i为质数,eg:2、3(质数)、4(循环体内被2筛除,minP[4]=2)
		if (!minP[i])prime[cnt++] = minP[i] = i;
		for (int j = 0; j < cnt&&i*prime[j] <= n; j++)
		{
			minP[i*prime[j]] = prime[j]; //分别筛除 i*prime[j],即将当前已知质数乘以i倍筛去,记录下它存在最小质因子
			if (i%prime[j] == 0)break;	//等价于//if(prime[j]>minP[i])break;
			//这一步是保证每个合数仅被其最小因子筛选一次的关键
		}
	}
}

eg:12=2*3*3,即12可被其质因子,2(2*6=12)和3(3*4=12)筛去
当prime中存在质数2,3(此时筛除掉了4,6,9)时,i循环到为4,
2*i=8筛除掉了8,因为i*prime[j]即4%2=0,退出循环,
假设没有这一break语句,下一步就是3*i=12,筛除掉了12,然后内层循环结束
然后i=5,筛除(10,15,25),接下来i=6,筛除(12,18,30)
有没有发现i=6的时候12又被筛除了一次。当i=4时被3筛除,当i=6时被2筛除
如果加上if (i%prime[j] == 0)break;就能保证12只会被其最小质因子2筛除
该break语句也可等价于if(prime[j]>minP[i])break;
还是以4为例子,当4*3筛去12时,发现3>minP[4](即2),即3大于12的另一个质因子2,故break;
因此可保证每个数i只被它最小质因子筛选。

五、区间筛

给定整数 a 和 b,请问区间[a, b) 内有多少个素数?(a < b≤1012, b−a≤106

区间筛实际上是利用了埃氏筛,只不过把范围偏移和缩小了一下。
根据素数的储备知识,b以内的合数的最小质因数一定不超过sqrt(b),
如果有sqrt(b)以内的素数表的话,就可以把埃氏筛法运用在[a, b)上了。
也就是说,先分别做好[2, sqrt(b) )的表和[a, b)d的表,
然后从[2, sqrt(b) )的表中筛得素数的同时,也将其倍数从[a, b)的表中划去,
其中[sqrt(b),a)范围内的数不需要筛选
最后剩下的就是区间[a, b)内的素数了。(《挑战程序设计竞赛(第二版)》)

当[a,b)中b的值很大,而b-a不大时可以利用区间偏移减少内存开销
例如要求 [1e8,1e8+10) 区间内的素数,难道我们要开 1e8+8 大小的数组吗?
利用下标偏移,可以减少空间开销。即将 [1e8,1e8+10) 对应到[0,10),整体区间向左偏移 1e8
代码中也有详细的注释

typedef long long ll;
const int MAXN = 1e5;
bool isPrime[MAXN];	//根据isPrime[i]判断i+a是不是质数,isPrime[i-a]=true 表示i是素数(下标偏移了a)
bool isSqrtPrime[MAXN];	//根据isSqrtPrime[i]判断i是不是质数

//对区间[a,b)内的整数执行筛法 
void SegmentSieve(ll a, ll b)	//对区间[a,b)进行筛选 
{
	for (ll i = 0; i * i < b; i++) //初始化[2,sqrt(b))的全为质数
		isSqrtPrime[i] = true;
	for (ll i = 0; i < b - a; i++) //初始化偏移后的[a,b)全为质数
		isPrime[i] = true;

	for (ll i = 2; i * i < b; i++) //埃氏筛选
	{ 
		if (isSqrtPrime[i]) {
			for (ll j = 2 * i; j  <= sqrt(b); j += i)//筛选上限变为sqrt(b)
				isSqrtPrime[j] = false;
			//筛除[a,b)里的非素数
			//(a+i-1)/i 得到最接近(a/i)这个数,也即(a的i倍数),最低是2LL是表示最低2倍
			for (ll j = max(2LL, (a + i - 1) / i) * i; j < b; j += i)
				isPrime[j - a] = false;
		}
	}
}

六、完整的测试代码:

#include<iostream>
#include<algorithm>
typedef long long ll;
using namespace std;
const int MAXN = 1e5;
bool primeArr[MAXN];	//用于朴素筛和埃氏筛

//用于欧拉筛
int minP[MAXN];//minP[i] 为i的最小素因子 
int prime[MAXN];//prime[i] 代表n以内从小到大,第i个质数 
int cnt;//素数计数 

//用于区间筛
bool isPrime[MAXN]; //isPrime[i-a]=true 表示i是素数(下标偏移了a)
bool isSqrtPrime[MAXN];

//朴素筛选n以内的素数
//判断k是不是素数
bool Judge(int k)
{
	if (k == 2)return true;
	if (k < 2 || k%2==0)return false;
	for (int i = 3; i <= sqrt(k); i += 2)
	{
		if (k%i == 0)return false;
	}
	return true;
}
//朴素筛
void PurifiedSieve(int n)
{
	for (int i = 0; i <= n; i++)
	{
		primeArr[i] = Judge(i);
	}
}

//埃氏筛
void EraSieve(int n)
{
	for (int i = 0; i < n; i++)
		primeArr[i] = true; //初始化每个数都为素数
	primeArr[0] = primeArr[1] = false;	//初始化0和1不是素数
	for (int i = 2; i <= sqrt(n); i++)	//从i从2到根号n作为质数,依次剔除它的2、3、4……倍
		if (primeArr[i])	//若这个数是素数则剔除它的整数倍
			for (int j = i*2; j <= n; j += i)	//剔除整数倍,上限是n(这里会有重复剔除的数,欧拉筛可优化)
				primeArr[j] = false;
}

//欧拉筛
void EulerSieve(int n)
{
	for (int i = 2; i <= n; i++)
	{
		//从2开始,若i不为0,则表示i为质数,eg:2、3(质数)、4(循环体内被2筛除,minP[4]=2)
		if (!minP[i])prime[cnt++] = minP[i] = i;
		for (int j = 0; j < cnt&&i*prime[j] <= n; j++)
		{
			minP[i*prime[j]] = prime[j]; //分别筛除 i*prime[j],即将当前已知质数乘以i倍筛去
			if (i%prime[j] == 0)break;	//等价于//if(prime[j]>minP[i])break;
			//这一步是保证每个合数仅被其最小因子筛选一次的关键
		}
	}
}

//区间筛
void SegmentSieve(ll a, ll b)	//对区间[a,b)进行筛选 
{
	for (ll i = 0; i * i < b; i++) //初始化[2,sqrt(b))的全为质数
		isSqrtPrime[i] = true;
	for (ll i = 0; i < b - a; i++) //初始化偏移后的[a,b)全为质数
		isPrime[i] = true;

	for (ll i = 2; i * i < b; i++) //埃氏筛选
	{ 
		if (isSqrtPrime[i]) {
			for (ll j = 2 * i; j  <= sqrt(b); j += i)//筛选上限变为sqrt(b)
				isSqrtPrime[j] = false;
			//筛除[a,b)里的非素数
			//(a+i-1)/i 得到最接近(a/i)这个数,也即(a的i倍数),最低是2LL是表示最低2倍
			for (ll j = max(2LL, (a + i - 1) / i) * i; j < b; j += i)
				isPrime[j - a] = false;
		}
	}
}

//打印primeArr[],即打印朴素筛和埃氏筛的结果
void PrintPrime(int n)
{
	int count = 0;
	for (int i = 0; i <= n; i++)
	{
		if (primeArr[i]) {
			count++;
			printf("%5d ", i);
			if (count % 10 == 0)printf("\n");
		}	
	}
}

//主函数,测试
int main()
{
	printf("\n=============================朴素筛=========================\n");
	PurifiedSieve(200);
	PrintPrime(200);
	printf("\n=============================埃氏筛=========================\n");
	EraSieve(200);
	PrintPrime(200);
	printf("\n=============================欧拉筛=========================\n");
	EulerSieve(200);
	for (int i = 0; i < cnt; i++)
	{
		printf("%5d ", prime[i]);
		if ((i + 1) % 10 == 0)printf("\n");
	}
	printf("\n=============================区间筛=========================\n");
	SegmentSieve(100, 200);
	for (ll i = 0,cnt = 0; i < 200 - 100; i++)
		if (isPrime[i])
		{
			printf("%5lld ", i + 100);
			if ((cnt++ + 1) % 10 == 0)printf("\n");
		}   
	return 0;
}

参考:
https://zh.wikipedia.org/wiki/%E5%9F%83%E6%8B%89%E6%89%98%E6%96%AF%E7%89%B9%E5%B0%BC%E7%AD%9B%E6%B3%95
https://www.cnblogs.com/Jackpei/p/10372392.html
https://blog.csdn.net/nixinyis/article/details/64131351
https://www.wmathor.com/index.php/archives/1182/

  • 6
    点赞
  • 18
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值