素数的计算 - 从试除到筛法
昨天搜索素数的问题的时候，找到一篇很棒的文章，转载一下，并加上一些自己的理解。
文章链接: 素数的计算: 从试除到筛法(C++实现)
素数定理：
素数的个数是有规律的，对于正实数x，定义T(x)为素数计数函数：不大于x的素数的个数。
那么会有：T(x) 约等于 x / ln(x) .
其中 T(x) / (x/ln(x)) 总是小于1.17。这个性质能够帮我们确定x以内的最大素数个数，以分配存储空间(即数组开多大)。
对于找不大于x的所有素数，有很多方法，下面从最简单的到最高效的介绍。
试除法
思想：
判断方法：对于一个数n，将其分别除以[2,n]的每一个整数，如果有可以整除的，则不是素数。
循环所有不大于x的数，确定数xi是否为素数以求出不大于x的所有素数。
算法实现:
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
// 判断n是否为素数
bool is_prime(int n)
{
if (n < 2)
return false;
for (int i = 2; i < n; ++i)
if (n % i == 0)
return false;
return true;
}
// 计算所有不大于n的素数
void get_prime(vector<int>& prime, int n)
{
for(int i = 2; i <= n; ++i)
if(is_prime(i)) // 判断i是否是素数
prime.push_back(i);
}
int main()
{
int n = 100000;
vector<int> prime;
get_prime(prime, n);
return 0;
}
算法复杂度: O(n)
试除法优化
思想：
对上面思想的分析，我们发现，若有数n = x*y，那么y与x中必有一个满足：k <= sqrt(n)。即在sqrt(n)之前没有找到能够整数n的数，那么n之后也不会有，所以对于判断n是否为素数的方法中的循环，只需要判断到sqrt(n)即可。
算法实现
更改is_prime方法：
bool is_prime(int n)
{
if (n < 2)
return false;
for (int i = 2; i * i <= n; ++i)
if (n % i == 0)
return false;
return true;
}
算法复杂度： O(n * sqrt(n))
上面的都是暴力手法，下面介绍科学的艺术手法。
埃氏筛法
埃拉托斯特尼筛法，是古希腊数学家发明的计算素数的方法(妈耶，古希腊就研究出来了)。
思想：
对于求解不大于n的素数：
找出不大于sqrt(n)内的素数 p1、p2、p3 .... pn (1 <= n <= sqrt(k))
依次剔除不大于n的pi的倍数。
剩下的都是素数。
这里有一张埃氏筛法的工作原理图：
算法实现：
void get_prime (vector<bool> &isPrime, int n) {
isPrime.assign(n + 1, true);
if (n < 2) return;
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
// 若计算规模较大，加上 i < sqrt(n) 的判断条件
for (int j = i * i; j <= n && i < sqrt(n); j += i) {
isPrime[j] = false;
}
}
}
}
这里有一个原文没有提到的点，在我的实践中，当 n > 46000 时，会出现 i*i 超出int范围的问题，此时j会是负值，造成非法数组访问的错误。所以如果计算规模较大，加上 i < sqrt(n) 的判断条件。
算法时间复杂度： O(nloglogn)
欧拉筛法
思想
埃氏筛法会对两个素数的公倍数多次剔除，根据这一问题优化后，便出现了欧拉筛法：对于一个没有被筛过的数，只要被第一个数筛了就行了。
算法实现
void get_prime(vector<int> &prime, vector<bool> &isPrime, int n) {
isPrime.assign(n + 1, true);
int max_num = (n / log(n)) * 1.17 + 1;
if (n < 2) return;
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
prime.push_back(i);
}
for (int j = 0; j < isPrime.size() && i * prime[j] <= n; j++) {
isPrime[i*prime[j]] = false;
if (i % prime[j] == 0) break;
}
}
}
这个算法的精髓主要在于：
if (i % prime[j] == 0) break;
证明：对于一个数i，若能整除prime[j]，那么有 i = a * prime[j]。当判断prime[j+1]能否被i整除时，有 i * prime[j + 1] = a * prime[j] * prime[j+1]，即此数(i*prime[j+1])已经被prime[j]剔除过了，所以跳出。
此算法时间复杂度：O(n)
数组开多大
最前面已经说了，数组开多大可以用T(x) / (x/ln(x)) <= 1.17 来决定，那么：
int max_num = (n / log(n)) * 1.17 + 1;
