素数判断相关算法

问题引入

输入一个数n，判断n是否为素数。其中（$2\le n\le 10^7$）

解决方案:

#include<iostream>
using namespace std;
int n;
int main(){
	scanf("%d",&n);
	for(int i=2;i<n;i++){
		if(n%i==0){
			cout<<"不是质数\n";
			return 0;
		}
	}
	cout<<"是质数\n";
	return 0;
}

这段代码就是挨个挨个的尝试，从$2$到$n-1$都除一下，如果余数为$0$,那么说明有因子，就不是质数，直接$return0$结束程序，然后我们要知道计算机每秒大概能运行$10^7$到$10^8$次，所以上面的解决方案能够在一秒内完成.

现在我们对这个问题进行小小的修改

输入一个数n，判断n是否为素数。其中（$2\le n\le 10^{16}$）

如果我们还是一个一个的去尝试的话，大概需要$1000$天，在大家大学毕业之前应该能判断出来……所以我们需要改进，我们想一想，因为是相除取余数，假设$n$是一个质数，那么$n$一定能都被表达成$n=x_1\ast x_2$，然后$n$对其中任意一个因子取余，都可以得到$n\%x_1==0$，或者$n\%x_2==0$然后只用找到最小的那个因子就行了，所以$min(x_1,x_2)$最大等于$\sqrt{n}$，所以我们只需要枚举到$\sqrt{n}$就可以了

#include<iostream>
#include<cmath>
using namespace std;
long long n;
int main(){
	cin>>n;
	for(int i=2;i<=sqrt(n);i++){
		if(n%i==0){
			cout<<"不是质数\n";
			return 0;
		}
	}
	cout<<"是质数\n";
	return 0;
}

现在我们再对问题做一个小小的修改

请求出$[10^6,5\ast 10^6]$内的所有质数

解决方案:

#include<iostream>
#include<cmath>
using namespace std;
bool check(int x){
	for(int i=2;i<=sqrt(x);i++){
		if(x%i==0) return false;
	}
	return true;
} 
int main(){
	for(int i=1000000;i<=5000000;i++){
		if(check(i)) cout<<i<<endl;
	}
	return 0;
}

这里我们把判断($check$)函数单独列出来,让程序更有可读性,利用我们之前所说的方法进行一个判断.其实还是蛮快的,我们现在想看看他到底有多快,因为将数据输出到应用窗口是一个较慢的过程,所以我们把$cout$函数给去掉,只进行判断,或者在$for$循环的前面加上一行代码 freopen(“data.out”,”w”,stdout);将数据输出到文件里面,然后我们就能够发现,大概需要$5$秒左右的时间,不再像之前那样刚开始运行就出结果了.

现在我们又对问题做一个小小的修改

请求出$[2,3\ast 10^7]$内的所有质数

如果还是使用之前的算法,大概需要$70$秒的时间,所以我们考虑优化一下操作,使得程序变得更快.

介绍一种筛质数的方法:欧拉筛

我们先思考,所有的$2$的倍数都不是质数,所有$3$的倍数都不是质数,所有$5$的倍数都不是质数……所以我们对于每次找到的一个质数,我们标记它的所有的倍数的数字(小于$3\ast 10^7$)为合数,下次要判断这个数是不是质数的时候,我们就直接跳过它就行了,这样已经足够快了,但是还是有一个问题,就是同一个数字可能会被多个质数标记为合数,比如$30$会被$2,3,5$同时标记为质数,当一个数字很大,因子很多的时候,被标记的次数也会很多,就会大大的降低效率,所以我们要减少重复筛选.现在我们让一个合数只被自己最小的因子给筛掉,假设当前的枚举到的这个数字是$i$,小于$i$的素数集合是$prime[]$,那么我们就用$prime[]$集合里面的数和$i$来除掉“未来”的合数,$i\ast prime[]$一定是合数,然后在$i$能整除当前枚举到的$prime[j]$的时候就$break$掉,因为当前这个$i$可以分解成$prime[j]\ast k$的形式,所以如果再向后枚举,$i\ast prime[j+1]=prime[j]\ast k\ast prime[j+1]$,所以$prime[j+1]$不是筛掉的数的最小的因子,所以会重复筛去,所以我们就直接$break$掉就好了.

#include<iostream>
#define maxn 100000100
using namespace std;
int Max,n,prime[maxn],tot;
bool vis[maxn];
void getprime(int Max){
	vis[1]=true;
	for(int i=1;i<=Max;i++){
		if(!vis[i]) prime[++tot]=i;
		for(int j=1;j<=tot&&prime[j]*i<=Max;j++){
			vis[prime[j]*i]=true;
			if(i%prime[j]==0) break;
		}
	}
}
signed main(){
	freopen("data.out","w",stdout);
	Max=30000000; 
	getprime(Max);
	for(int i=1;i<=tot;i++) printf("%d\n",prime[i]);
	return 0;
}

然后我们加上freopen之后看看时间效率,还是挺快的,只需要3秒钟,而且是在数据范围扩大了的情况下只需要三秒钟.如果用cout进行输出的话,则需要20秒钟,因为cout在输出的时候需要电脑去判断数据类型,所以比较慢.用上一个方法,并且用printf输出会花50秒钟左右的时间,所以在筛素数方面就相差了47秒,数据范围扩大,这个差距将会更加的明显.

问题引入:

输入一个数字$n$,接下来输入$n$个数字$a[i]$,判断$a[i]$是否为质数

其中$a[i]\le 10^{18},n<10^5$

显然,不可能去用欧拉筛,因为数组开不下,因为需要$9\ast 10^5$T的储存空间,而且运算$10^{18}$次也需要$300$年左右的时间,用$\sqrt{a_i}$的算法也需要$11$天.

所以最后介绍一种叫做Miller_Rabin的算法:

• 费马小定理:如果$p$是一个质数并且$a$不是$p$的倍数,则$a^{p-1}\equiv 1(modp)$

但是这个定理是有缺陷的,卡迈克尔数就是这个定理的反例,但是卡麦尔数毕竟只是少部分的数字大部分的数字还是能够通过费马小定理的

• 二次探测定理:若$p$为质数,$a^2\equiv 1(modp)$

所以$a=1$或者$a=p-1$

如果$p$是质数,我们就能写成$2^k\ast t$的形式,所以有$a^{2^k\ast t}\equiv 1(modp)$,然后我们依次对$t,2t,4t.....t\ast 2^k$进行二次探测看看是不是全部都是$1$或者$p-1$,只要中间有一个不满足,那就说明是合数,否则就是质数,但是这样下来正确率是$0.75$,所以我们多取几个$a$,就能使正确率大大提高,当我们取$a=2,3,5,7,11,13$的时候,正确率就已经是$0.999755859375$了,对于计算机来说,就是一个正确并且优秀的算法.

#include<cstdio>
#include<algorithm>
#include<cmath>
#define int long long
using namespace std;
int a[6]={2, 3, 5, 7, 11, 13},n;
inline int pow_add(int x, int y, int p){
    return ((x*y-(int)(((long double)x*y+0.5)/p)*p)%p+p)%p;
} 
inline int pow_mul(int a, int d, int p){
    int ans = 1;
    while(d){
        if(d % 2) ans = pow_add(ans, a, p);
        a = pow_add(a, a, p);
        d >>= 1;
    }
    return ans % p;
}
inline bool test(int a, int p, int d){
    if(!(p % 2)) return false;
    while(!(d & 1)) d >>= 1;
    register int t = pow_mul(a, d, p);
    while(d != p - 1){
        int y = pow_add(t, t, p);
        if(y == 1 && t != 1 && t != p - 1) return false;
        t = y; d <<= 1;
    }
    return t == 1;
}
inline bool Miller_rabin(int x){
    for(int i = 0; i < 6; i ++){
        if(a[i] == x) return true;
        if(!test(a[i], x, x - 1))
           return false;
    } 
    return true;       
}
signed main(){
	scanf("%lld",&n);
	if(Miller_rabin(n)) printf("Yes");
	else printf("No");
    return 0;
} 

其中$pow\_mul$是快速幂,大家可以了解一下,$pow\_add$是避免溢出.