Cover

素数的计算 - 从试除到筛法

昨天搜索素数的问题的时候，找到一篇很棒的文章，转载一下，并加上一些自己的理解。

文章链接: 素数的计算: 从试除到筛法(C++实现)

素数定理：

素数的个数是有规律的，对于正实数x，定义T(x)为素数计数函数：不大于x的素数的个数。

那么会有：T(x) 约等于 x / ln(x) .

其中 T(x) / (x/ln(x)) 总是小于1.17。这个性质能够帮我们确定x以内的最大素数个数，以分配存储空间(即数组开多大)。

对于找不大于x的所有素数，有很多方法，下面从最简单的到最高效的介绍。

试除法

思想：

1. 判断方法：对于一个数n，将其分别除以[2,n]的每一个整数，如果有可以整除的，则不是素数。

2. 循环所有不大于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的素数：

1. 找出不大于sqrt(n)内的素数 p1、p2、p3 .... pn (1 <= n <= sqrt(k))

2. 依次剔除不大于n的pi的倍数。

3. 剩下的都是素数。

这里有一张埃氏筛法的工作原理图：

埃氏筛法

算法实现：

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;