今天的网络赛再一次告诉我——数学真的很重要——队友推了两小时公式O(1)一发AC。所以今天就让我们来学习一下一些经典的素数相关算法吧！
素数
在程序设计竞赛中，我们很少碰到单纯关于素数的题目，但经常会在题目的某一部分用到关于素数的知识。而所谓素数，就是值自然数当中，除了1之外，只能被1和该数本身整除的数，下面我们就来一同学习求解素数的相关算法吧。
筛法求素数
朴素筛法求素数
我们先来看一下代码：
#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#define SIZE 10000000
using namespace std;
int check[SIZE] = {0};//元素值为0代表是素数
int prime[SIZE] = {0};
int main() {
int pos = 0;
int flag;
for (int i = 2 ; i < SIZE ; i++) {
if (!check[i]) {
//如果是素数
prime[++pos] = i;
}
for (int j = 2 * i ; j <= SIZE ; j += i) {
check[j] = 1;
}
}
printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
//for(int i = 1; i <= pos; i++)printf("%d ", prime[i]);
return 0;
}
//Output:1.17这种方法比较好理解，初始时，我们假设全部数都是素数，当找到一个素数时，即check[i]==0时，我们就把它加入素数表内，并且筛去掉它的倍数。显然这个素数乘上另外一个大于等于二的正整数之后得到的数——即它的倍数都是合数，把这些合数都筛掉，就首先了素数的筛选，这也是算法名字——“筛选法”的由来。而这种简单粗暴的方法在1e7的数据规模下平均用时约为1.17s。
其实就算是朴素的筛法也有优化：
#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#define SIZE 10000000
using namespace std;
int check[SIZE] = {0};//元素值为0代表是素数
int prime[SIZE] = {0};
int main() {
int pos = 0, num = int(sqrt(SIZE));
int flag;
for (int i = 2 ; i < SIZE ; i++) {
if (!check[i]) {
//如果是素数
prime[++pos] = i;
}
if(i <= num)
for (int j = i * i ; j <= SIZE ; j += i) {
//i*i的优化快了不少
check[j] = 1;
}
}
printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
//for(int i = 1; i <= pos; i++)printf("%d ", prime[i]);
return 0;
}
//Output:0.66我们把删去素数倍数的操作从i2变成了ii，就快了近一倍，为什么可以这样优化呢？你想，我们假定i为大于2的一个素数，那么按照我们之前的筛法i2就在处理到素数2的时候筛了一遍，处理到素数i的时候又处理了一遍，造成了不必要的时间浪费，而改成ii，则可以从大于等于当前素数的倍数开始，从而避免了这个问题。
线性筛法求素数
我们先来看一下代码：
#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#define SIZE 10000000
using namespace std;
int prime[SIZE], isNotPrime[SIZE] = {1, 1};
int n, m, primeNum = 0, num;
int main() {
for(int i = 2; i <= SIZE; i++) {
if(!isNotPrime[i]) prime[++primeNum] = i;
//key1
for(int j = 1; j <= primeNum && i * prime[j] <= SIZE; j++) {
isNotPrime[i * prime[j]] = 1;
if(i % prime[j] == 0) break;//key2
}
}
printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
//for(int i = 1; i <= primeNum; i++)printf("%d ", prime[i]);
return 0;
}
//Output:0.18线性筛的最大优点就是每个数最多只被筛一次，也就是说时间复杂度是O(n)的。但其实线性筛素数的思想很简单，首先，我们要先明确一个条件，任何一个合数都能表示成一系列素数的积。
接着我们来看程序不管 i 是否是素数，都会执行到“key1”，这时我们分情况讨论：
①如果此时 i 是素数的话，我们进入j循环。那很简单，被筛选出来的合数就是一个大的素数 i 乘以不大于 i 的素数 prime[j]，这样筛选出来的数跟之前的是不会重复的。此时筛出的数都是 N=p1p2的形式, p1，p2之间不相等且均为素数；例如35 = 5 7 只会在 i = 7 的时候被筛选出来；
②如果此时 i 是合数的话，我们进入j循环。此时 i 可以表示成递增素数相乘 i=p1p2...pn, pi都是素数（2<=i<=n）， pi<=pj ( i<=j )，其中p1是最小的系数。根据“key2”的定义，当p1==prime[j] 的时候，筛除就终止了，也就是说，只能筛出不大于p1的质数i，这样筛出来的数也是之前没筛出来过的。我们可以直观地举个例子。i=235，此时能筛除 2i ,不能筛除 3i；如果能筛除3i 的话，当 i' 等于 i'=335 时，筛除2i' 就和前面重复了。
所以这样我们就在线性的时间里把素数筛选出来啦。
优化线性筛法求素数
我们先来看代码：
#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#define SIZE 10000000
using namespace std;
const int HALF = SIZE / 2;
int prime[HALF] = {0, 2}, isNotPrime[HALF] = {1};
//isNotPrime[i] = 0 表示 i * 2 + 1 是素数
int n, m, primeNum = 1, num;
int main() {
for(int i = 1; i <= HALF; i++) {
if(!isNotPrime[i]) prime[++primeNum] = 2 * i + 1;
//key1
for(int j = 2; j <= primeNum && (2 * i + 1) * prime[j] <= SIZE; j++) {
isNotPrime[((2 * i + 1) * prime[j] - 1) / 2] = 1;
if((2 * i + 1) % prime[j] == 0) break;//key2
}
}
printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
//for(int i = 1; i <= primeNum; i++)printf("%d ", prime[i]);
return 0;
}
//Output:0.10什么你觉得刚刚的线性筛还不够快？没事这里我想出了一种可以在约O(n/2)筛选出素数的线性筛。
其大致思路和普通的线性筛一样，只不过我们知道，除了2以外的偶数都是合数，所以这次我们只处理奇数。在key1后筛选合数的时候，偶数*质数 = 偶数 = 合数，所以我们直接去掉所有的偶数，只处理奇数也能筛选出全部的素数，这样我们就少处理了一半的数，实测运行时间也比之前的线性筛快了快一半。
素数的判定
朴素算法
素数的判定最原始的方法就是按照素数的定义进行判定，看看要判定的数存不存在除了1和它本身之外的其他因子：
bool isPrime(int x) {
for(int i = 2; i < x; i++) {
if(x % i == 0)return(false);
}
return(true);
}改良算法
刚刚的那种算法时间复杂度是不太可观的。我们知道对于任意一个合数 x ，它一定可以分解成两个大于1的整数 a 和 b 的积，即 x = a * b；我们很容易发现，a 和 b 当中至少有一个数一定小于等于 x 的平方根。所以我们可以对算法进行如下改进：
bool isPrime(int x) {
for(int i = 2; i * i <= x; i++) {
if(x % i == 0)return(false);
}
return(true);
}更优算法
在上面的算法中，其实我们也只需要找到 x 是否存在质因子即可，也就是说，从 2 到 x 的平方根的那些合数其实也是没有必要进行判断的。所以在要多次进行判断素数时，我们可以提前进行线性筛生成一个素数表，再进行判断，这样循环的次数还可以降低：
//int list[1100000] = {0, 2, 3, 5, ………………};
//如果需要判断的质数最大为 10^12, 只需要筛选√10^12 = 10^6 内的素数即可
bool isPrime(int x) {
for(int i = 1; list[i] * list[i] <= x; i++) {
if(x % list[i] == 0)return(false);
}
return(true);
}米勒-勒宾测试
上面的几种方法都是针对需要判定的数较小的情况，如果需要判定的数大到一定规模，上面的方法就不太适用了，这时可以考虑米勒-勒宾测试法，其描述如下：
令大奇数 n = 2^p * m + 1，其中 p 是非负整数，m 是正奇数。若 b^m ≡ 1（mod n） 或 b^（2^j * m） ≡ 1（mod n）,0 ≤ j ≤ p - 1，则称 n 通过以b为基的米勒-勒宾测试。若n为素数，b为整数，且 n 不能 整除 b，则 n 必然通过以b为基的米勒-勒宾测试。故代码如下：
#include<stdio.h>
#include<iostream>
#include<algorithm>
using namespace std;
long long mul(long long a, long long b, long long mod){//伪高精度乘法
a %= mod;
b %= mod;
long long c=(long double)a * b / mod;
long long ans = a * b - c * mod;
return (ans % mod + mod) % mod;
}
long long pow_mod(long long x, long long n, long long mod){//快速幂
long long res=1;
while(n){
if(n & 1)
res = mul(res, x, mod);
x = mul(x, x, mod);
n >>= 1;
}
return (res + mod) % mod;
}
bool Miller_Rabbin(long long a,long long n){
//把n-1 转化成 (2^r)*s
long long s = n - 1, r = 0;
while((s & 1) == 0){
s >>= 1;
r++;
}
//算出 2^d 存在 k 里
long long k = pow_mod(a, s, n);
//二次探测 看变化过程中是不是等于1 或 n-1
if(k == 1)return true;
for(int i = 0; i < r; i++,k = k * k % n){
if(k == n - 1)return true;
}
return false;
}
而且据说按照下表取b的值 在一定范围内 可以使正确率达到百分之百，变成一个确定的算法：
if n < 1 373 653 b = 2 and 3.
if n < 9 080 191 b = 31 and 73.
if n < 4 759 123 141 b = 2, 7, and 61.
if n < 2 152 302 898 747 b = 2, 3, 5, 7, and 11.
结语
通过这篇BLOG，相信你对素数也有了一些初步的了解。最后希望你喜欢这篇BLOG！
非技术的路过。