（初学者）素数计算法

最新推荐文章于 2022-06-10 15:17:51 发布

GGN_2015

最新推荐文章于 2022-06-10 15:17:51 发布

阅读量8.2k

点赞数 6

分类专栏：算法导论文章标签：算法数论基础素数欧拉线筛

本文链接：https://blog.csdn.net/GGN_2015/article/details/68059353

版权

算法导论专栏收录该内容

76 篇文章 0 订阅

订阅专栏

前言

记得我刚开始学习OI的时候接触到的最早的算法莫过于素数算法。当然了，素数这个知识小学的时候我们就已经有所学习了。

质数(prime number)又称素数，有无限个。一个大于1的自然数，除了1和它本身外，不能被其他自然数整除，换句话说就是该数除了1和它本身以外不再有其他的因数;否则称为合数。

当时我一直有一个疯狂的设想，就是编写一个程序把我能算出来的素数都算出来。因为当时水平有限，用C语言打印了一个一百万以内的质数表都用了好几秒的时间，所以说我一直希望能有一种方法又快而又准确地求出很大范围之内的素数。

1.朴素的质数验证法

我们当时学的最早的那种方法就是“朴素算法”，这是一个用来验证一个数是否是质数的算法。它的原理就是和小学数学书上所写是一模一样的：如果这个数是2，那么它是质数；如果这个数大于2，那么就判断它是否能整除从2到这个数减一中的所有数，如果能则是质数。

代码也是非常简洁的（simple and stupid）：

bool IsPrime(int NumToCheck)//判断一个数是否是质数
{
	if(NumToCheck<=1)//不考虑小于等于1的情况
		return 0;
	if(NumToCheck==2)//如果是2,则是质数
		return 1;
	for(int i=2;i<NumToCheck;i++)
		if(NumToCheck%i==0)//对于从2到n-1的每一个数,一旦整除,则说明不是质数
			return 0;
	return 1;//循环完毕,说明都不能整除,则是质数
}

有了这个朴素算法，我们就可以用最暴力的方法打印出一个质数表来。对于每一个数，判断它是否是质数，如果是，则输出这个数。

看代码：

void OutputPrimeTable(int NumMax)
{
	for(int i=2;i<=NumMax;i++)
		if(IsPrime(i))
			cout<<i<<" ";
	cout<<endl;
}

这个算法虽然很慢，但是当NumMax<=一万时，它的效率还是可以令我们接受的。

程序运行的效果图

因为，朴素算法的素数验证的时间复杂度为O(n)。所以，这个“朴素素数表”的时间复杂度为O(n^2)。

2.对朴素算法的优化

小学的时候，我们学过一个东西叫“分解质因数”。再分解质因数的过程中，你会发现一个数的因数可以理解为是成对出现的，例如： $12=1\times 12=2\times 6=3\times 4$ ，再比如： $16=1\times 16=2\times 8=4\times 4$ 。显然，如果一个数n有一个大于等于sqrt(n)的因数，那么它一定会有一个小于等于sqrt(n)的因数与之对应。所以说，我们在验证素数时可以只枚举小于等于sqrt(n)的数就可以得到正确的答案。

但是一定要注意是“小于等于”根号n，而不是“小于”根号n。就比如说：9=3*3，如果你只验证了“小于”根号n的数（也就是只验证了2），那么你的程序会把9当成一个素数（而实际上它不是）。

这就是一种“优化的朴素算法”，请看代码：

#include<cmath>//引用数学头文件
bool IsPrimeSqrt(int NumToCheck)
{
	if(NumToCheck<=1)
		return 0;
	if(NumToCheck==2)
		return 1;
	int SqrtN=(int)sqrt(NumToCheck);
	for(int i=2;i<=SqrtN;i++)//只判断小于等于sqtr(n)的数
		if(NumToCheck%i==0)
			return 0;
	return 1;
}

而这种方法的效率其实已经可以满足我们"日常生活"的需要了。因为在n小于等于一百万的时候，用它打印质数表的效率还是非常可观的。我带领同学写一个程序验证一下这个观点。

#include<iostream>
#include<cstdlib>
#include<cmath>
bool IsPrimeSqrt(int NumToCheck)//升级版朴素质数判断函数
{
	if(NumToCheck<=1)
		return 0;
	if(NumToCheck==2)
		return 1;
	int SqrtN=(int)sqrt(NumToCheck);
	for(int i=2;i<=SqrtN;i++)
		if(NumToCheck%i==0)
			return 0;
	return 1;
}
void OutputPrimeTableSqrt(int NumMax)//（去掉了输出的）打印质数表函数
{
	for(int i=2;i<=NumMax;i++)
		if(IsPrimeSqrt(i));//在这里我们不对计算出的质数进行输出,只是计算
		//我们只统计计算所用的时间,因为输出是非常耗时的
		//我们在可以把这个计算出的素数表储存到内存中,这个O(1)的复杂度我们可以忽略不记
		
		//我为什么不把输入加进去,实际上是有一个故事的...
		//请耐心看到文章最后
}
#include<ctime>
void CheckTime(int NumMax)//下一个用于检测函数性能的程序
{
	clock_t start_t=clock();//计时开始
	OutputPrimeTableSqrt(NumMax);//运行函数
	clock_t end_t=clock();//计时结束
	double timeLast=(double)(end_t-start_t)/CLOCKS_PER_SEC;//计算时间
	cout<<"Check For NumMax="<<NumMax<<" LastTime="<<timeLast<<endl;//输出结论
}
int main()
{
	CheckTime(10000);
	CheckTime(100000);
	CheckTime(1000000);
	CheckTime(10000000);//分别对以上的四组数据进行时间测试
	system("pause");
	return 0;
}

当数据范围超过一百万时，这个算法就显得“捉襟见肘”了。有图有真相：

对算法时间的统计

所以，我们也许需要一些更高端的方法来弥补这个算法的缺点。而这个算法就是著名的：

3.（埃拉托斯特尼筛法）埃氏筛法

它并不像朴素一样，对每个数进行判断，因为它是一种“筛法”。就像是一个“筛子”中有一大堆数，每次把其中我确定出来的不是素数的数从“筛子”里筛出，最后剩下的就是我想要的素数了。首先，二是素数。然后每当我确定出一个素数，就把所有它的倍数从这个筛子里筛出。代码如下：

bool IsPrimeTable[100000001];//IsPrimeTable[i]表示i是否在筛子中
void PrimeSelect(int NumMax)//筛法
{
	IsPrimeTable[0]=IsPrimeTable[1]=0;//0和1排除
	for(int i=2;i<=NumMax;i++)
		IsPrimeTable[i]=1;//把所有数装入筛子
	for(int i=2;i<=NumMax;i++)
		if(IsPrimeTable[i])//如果已经循环到这个数,而它仍未被删除,那么它一定是质数
			for(int j=i*2;j<=NumMax;j+=i)//把除了它自身以外所有它的倍数从筛子中筛出
				IsPrimeTable[j]=0;
}
void OutputSelectPrime(int NumMax)//输出筛子中剩余的数
{
	PrimeSelect(NumMax);//筛法
	for(int i=1;i<=NumMax;i++)
		if(IsPrimeTable[i])//在筛子中则输出
			cout<<i<<" ";
	cout<<endl;
}

这种方法的效率还是很高的，NumMax<=10,000,000时它可以从容面对。但当NumMax=100,000,000时他可能就得需要几秒钟的时间来完成。这个算法的时间复杂度一定是大于O(n)的，这在朴素算法的基础上已经是一个飞跃了。

为什么大于O(n)呢？因为我们每找到一个素数就筛去所有它的倍数，而这时每一个合数会被它所有的质因子都筛一遍。假设每一个合数都只被筛去一遍，那么复杂度为O（n），所以这种筛法的复杂度一定是大于O（n）的。

那有没有O（n）的质数筛法呢？这便是神奇的：

4.欧拉筛法

那么，欧拉筛法又是如何实现O（n）筛选的呢？只让每个合数被筛出一次？只让每个合数被它的最小质因子筛出! 而又如何让一个数只被它的最小质因子筛出呢？它的思想是非常“深邃”的，先给同学们看代码，然后再进行解读：

bool vis[1000000001];//记录一个数当前是否被"筛出"过
int PrimeCount=0;//记录当前质数表的大小
int PrimeTable[50847537];//质数表，你可能会感觉我定义的内存可能有点少，实则不然
//因为我是预先知道10^9以内的素数的个数的
void AddPrime(int thePrime)//向质数表中添加一个质数
{
	PrimeTable[++PrimeCount]=thePrime;
}
void PrimeSelect(int maxNum)//欧拉筛法
{
	vis[0]=vis[1]=vis[2]=1;
	AddPrime(2);//把2加到质数表
	for(int i=3;i<=maxNum;i++)
	{
		if(!vis[i])//如果这个结点仍未被"筛出"
			AddPrime(i);//把它加入质数表
		vis[i]=1;//把这个数筛出
		for(int j=1;j<=PrimeCount && i*PrimeTable[j]<=maxNum;j++)
		{
			vis[i*PrimeTable[j]]=1;//把 这个数 与 所有已知质数 的积 "筛出"
			if(i%PrimeTable[j]==0)//关键在这里,好好想想为什么要这么写
				break;
		}
	}
	
}

欧拉筛法就是用这句话实现了“每个合数”只被筛出一次。

if(i%PrimeTable[j]==0)break;//其中i为任意数，PrimeTable[j]一定是一个质数。

（请立刻打开你的脑洞，因为下面的这段话不是很好理解。）
对于每一个大于1的数，我在它身上乘一个质数，得到的结果一定是一个合数。对于一个大于1的数，我在它身上乘上一个它自己的最小质因子，得到的合数的最小质因子，仍为原数的最小质因子。同理，让一个大于1的数乘上一个小于它最小质因子的质数，那么得到的新的合数的最小质因子一定是我乘上去的这个新的质数。因为PrimeTable中的质数是从小到大储存的，所以第一次满足关系式(i%PrimeTable[j]==0)的时候，PrimeTable[j]恰好为i的最小质因子。这样，我们就保证了每一个合数只能被它的最小质因子删去，故复杂度为O（n）。

再让我们探究一下欧拉筛法的运算效率。为了计算它的运行速度，我写了这样的一个程序：

void CheckTime(int maxNum)//计算时间
{
	for(int i=0;i<=maxNum;i++)
		vis[i]=0;
	PrimeCount=0;//初始化不计时
	clock_t start_t=clock();//开始计时
	PrimeSelect(maxNum);//运行筛法
	clock_t end_t=clock();//结束计时
	double lastTime=(double)(end_t-start_t)/CLOCKS_PER_SEC;//计算用时
	cout<<"Check time for maxNum="<<maxNum<<" time="<<lastTime<<"s"<<endl;
	//输出测试结果
}
int main()
{
	CheckTime(10000);
	CheckTime(100000);
	CheckTime(1000000);
	CheckTime(10000000);
	CheckTime(100000000);
	CheckTime(1000000000);
	//对于以上的六个数据范围分别计时
	system("pause");
	return 0;
}

最后的效果是这样的：

欧拉筛法的运行测试

欧拉筛法的效率就非常的可观了，当maxNum $\leq 10^8$ 时都属于日常可以接受的范畴，尽管等于 $10^8$ 时时间略微超出了1秒。这也就是为什么我们总把“ $10^8$ 数量级运算”当做一个OI的一个“高压线”，也因此O（n）算法一般的安全边界是在 $n=10^7$ 。

5.后记

以上就是我要给同学们分享的全部内容了。

人们对素数的探索如同大海捞针，对素数的了解也可以说是九牛一毛。虽然我们已经知道了素数的很多特性还有定理，但是从古至今一直都没有人能找到素数的确切的表达式。“革命尚未成功，同志仍需努力！”

粘贴一些，先贤对数学的探究作为“文化常识”，有兴趣的同学可以看一下：

（我最佩服的）质数无限定理：
我之所以佩服是因为这个定理是探究所有以下质数定理的前提。

质数的个数是无穷的。欧几里得的《几何原本》中有一个经典的证明。它使用了证明常用的方法：反证法。具体证明如下：假设质数只有有限的n个，从小到大依次排列为p1，p2，……，pn，设N=p1×p2×……×pn，那么，N+1是素数或者不是素数。

如果N+1为素数，则N+1要大于p1，p2，……，pn，所以它不在那些假设的素数集合中。

如果N+1为合数，因为任何一个合数都可以分解为几个素数的积；而N和N+1的最大公约数是1，所以N+1不可能被p1，p2，……，pn整除，所以该合数分解得到的素因数肯定不在假设的素数集合中。

因此无论该数是素数还是合数，都意味着在假设的有限个素数之外还存在着其他素数。所以原先的假设不成立。也就是说，素数有无穷多个。

其它相关内容：

在一个大于1的数a和它2倍之间（即区间(a, 2a]中）必存在至少一个素数。

存在任意长度的素数等差数列。（格林和陶哲轩，2004年）

一个偶数可以写成两个数字之和，其中每一个数字都最多只有9个质因数。（挪威布朗，1920年）

一个偶数必定可以写成一个质数加上一个合成数，其中的因子个数有上界。（瑞尼，1948年）

一个偶数必定可以写成一个质数加上一个最多由5个因子所组成的合成数。
【后来，有人简称这结果为 (1 + 5) (中国，1968年)】

一个充分大偶数必定可以写成一个素数加上一个最多由2个质因子所组成的合成数。
【简称为 (1 + 2) (中国陈景润)】

著名猜想：

哥德巴赫猜想：是否每个大于2的偶数都可写成两个素数之和？

孪生素数猜想：孪生素数就是差为2的素数对，例如11和13。是否存在无穷多的孪生素数？

斐波那契数列内是否存在无穷多的素数？

是否有无穷多个的梅森素数？

在n2与（n+1）2之间是否每隔n就有一个素数？

是否存在无穷个形式如X2+1素数？

黎曼猜想

6.彩蛋——测试算法运行时间

为什么我刚才说我坚决不把输出的部分加到测试时间的部分中，因为输出远比运算要慢。我曾经用“埃氏筛法”计算一亿以内的质数，并把它输出到一个叫“Prime.txt”的文件里。令我非常欣喜的是，它计算仅仅用了十几秒时间，但是输出…输出了半个多小时。程序运行了一会，我就不耐烦地就去吃晚饭了，但是吃完饭回来一看还没输出完。因此，算法用时的计算不能包含输入输出所用的时间（在线算法就无奈了）。

#include<iostream>
#include<cstdlib>
#include<ctime>
using namespace std;

int CheckTime(/*用于测试的数据*/)
{
	//初始化的内容写在这里
	time_t start_t=clock();//取开始时间
	//调用算法程序
	time_t end_t=clock();//去结束时间
	double lastTime=double(end_t-start_t)/CLOCKS_PER_SEC;//求得用时
	//注意:一定要用double变量,这样就能精确到小数点后三位(更精确的方法我也不会了)
	cout<<"lastTime="<<lastTime<<"s"<<endl;
}

下面的这个也是同理：

int CheckTime(const char* ProgramPath)//测试一个外部程序的耗时
{
	time_t start_t=clock();
	system(ProgramPath);//调用外部程序
	time_t end_t=clock();
	double lastTime=double(end_t-start_t)/CLOCKS_PER_SEC;
	cout<<"lastTime="<<lastTime<<"s"<<endl;
}