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

作者: 程点 | 来源:发表于2018-09-10 20:19 被阅读2次

什么是素数

素数又称质数，是指大于1的自然数中，除了1和它本身，不能被其它自然数整除的数字。1被定义为非素数。大于1的自然数，如果不是素数则为合数。上古大神欧几里得已经证明了有无穷多个素数存在（欧几里得定理）。

素数定理

引入维基百科的描述

在数论中，素数定理描述素数在自然数中分布的渐进情况，给出随着数字的增大，素数的密度逐渐降低的直觉的形式化描述

素数的个数是有规律可循的，对于正实数x，定义 $\pi (x)$ 为素数计数函数，即不大于x的素数的个数。
$\pi(x)$ 的估计方法有很多，这里介绍比较简单的一种:
$\pi(x) \approx \frac{x}{\ln x}$
其中 $\ln x$ 是x的自然对数。
维基百科中给出了 $x$ 从10到 $10^{25}$ ， $\pi(x)$ 与 $\frac{x}{\ln x}$ 的数值，可以看出随着 $x$ 的增长， $\frac{\pi(x)}{\frac{x}{\ln x}}$ 始终小于 $1.17$ （这里的比值很关键，我们后面需要用来估计用于存储素数的数组的预分配空间的大小）。

常规计算方法（试除法）及实现

基本思想

根据定义，给定一个自然数，先判断它是否是素数。若要计算出所有不大于n的素数，则从2到n，依次判断是否是素数。根据如何判断一个数是否是素数又有以下几种常规方法。

常规方法一

为了判断自然数x是否是素数，我们可以从2开始，分别除以不大于x的每个数，如果存在能整除x的数，则x不是素数。
常规方法一的c++实现

#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(\sum_{i=2}^{n}i^2) =O(n^2)$
实际上我们没必要比较这么多次，下面我们介绍方法二。

常规方法二

如果 $n$ 为合数，已知 $n=\sqrt{n} \times \sqrt{n}$ ，假设 $n=xy$ ，如果 $x>\sqrt{n}$ ，那么 $y$ 必然小于 $\sqrt{n}$ 。即我们只需要判断到 $\sqrt{n}$ 即可。
常规方法二的C++实现
实际上只需要根据方法一更改is_prime函数即可，其余代码完全一样。
is_prime函数

bool is_prime(int n)
{
    if (n < 2)
        return false;
    /**
     * 只需要判断到根号n即可，这里必须要使用小于等于符号
     * （不能仅仅使用小于符号，例如判断9）
     **/
    for (int i = 2; i * i <= n; ++i)
        if (n % i == 0)
            return false;
    return true;
}

该算法的时间复杂度为
$O(\sum_{i=2}^{n}\sqrt{i}) = O(n\sqrt{n})$
实际上我们的算法还有进步的空间，下面介绍常用的两种性能不错的筛法。

埃拉托斯特尼筛法及C++实现

基本思想

埃拉托斯特尼筛法（sieve of Eratosthenes ） 是古希腊数学家埃拉托斯特尼发明的计算素数的方法。对于求解不大于 $n$ 的所有素数，我们先找出 $\sqrt{n}$ 内的所有素数 $p_1,p_2,...,p_{k-1}, p_k$ ，其中 $k = \lfloor \sqrt{n}\rfloor$ ，依次剔除 $p_i(1 <= i <=k)$ 的倍数，剩下的所有数都是素数。
维基百科中有一张关于使用埃拉托斯特尼筛法求解素数的GIF图

埃拉托斯特尼筛法求解素数
下面我们来看一个来自维基百科的例子，如果求25内的所有素数
1.列出大于等于2小于等于25的所有数组成的序列：

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

去掉素数2的所有倍数：

2 3 5 7 9 11 13 15 17 19 21 23 25

去掉素数3的所有倍数：

2 3 5 7 11 13 17 19 23 25

去掉素数5的所有倍数：

2 3 7 11 17 19 23

因为 $5 = \sqrt{25}$ ，所以算法结束。剩下的所有数都是素数。

埃拉托斯特尼筛法的证明

证明之前我们先粗略介绍算术基本定理。算法基本定理的描述是：每个大于1的自然数均可写为素数的积，而且这些素因子按大小排列之后，写法仅有一种方式。然后开始我们的证明（反证法）：
假设正整数 $n$ 为通过埃拉托斯特尼筛法得到序列中的合数且 $n$ 不能被小于等于 $\sqrt{n}$ 的素数整除。因为 $n$ 是合数，则必然有 $xy=n$ ，如果 $1<x<\sqrt{n}$ ，那么 $\sqrt{n} < y < n$ ，又因为 $n$ 不能被小于等于 $\sqrt{n}$ 的素数整除，故 $x$ 为合数，根据算术基本定理,合数 $x$ 一定写成关于某个素数 $p$ 的积，因为 $1 < x < \sqrt{n}$ ，则 $1<p<\sqrt{n}$ 。与假设矛盾，证毕。

C++实现

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

void get_prime(vector<int> &prime, int n)
{
    // 初始化一个布尔数组，用于判断下标i是否是素数
    vector<bool> is_prime(n + 1, true);
    if (n < 2)
        return;
    for (int i = 2; i <= n; ++i)
    {
        if (is_prime[i])
        {
            prime.push_back(i); // 把素数加入到数组中
            // 这里只需要从 i * i 开始筛选即可，证明见下文
            for (int j = i * i; j <= n; j += i)
                is_prime[j] = false;
        }
    }
}
int main()
{
    int n = 100;
    vector<int> prime;
    get_prime(prime, n);
    // 打印计算出的素数
    for_each(prime.begin(), prime.end(), [](int &n) {
        cout << n << " ";
    });
    return 0;
}

在上述代码的get_prime函数中

for (int j = i * i; j <= n; j += i)

这一行中的的 $j$ 是从 $i^2$ 开始的，也就是说对于当前素数 $i$ 是从 $i^2$ 开始筛选的，且剩余序列中从 $2$ 到 $i^2-1$ 的所有数都是素数。证明如下：
如果 $i=2$ ， $i^2$ 之前的序列为: $2,3$ ，都是素数。
如果 $i>2$ ，因为 $1<i^2-1 < i^2$ ，所以 $i-1<\sqrt{i^2-1} < i$ ， $i-1$ 为 $i$ 的前一个数，即 $\sqrt{i^2-1}$ 大于 $i-1$ 及之前的所有数。因为 $i$ 之前的所有素数都筛选过了，且均小于 $\sqrt{i^2-1}$ ，根据埃拉托斯特尼筛法可知，序列中 $i^2$ 之前的所有数都是素数，所以 $i$ 只需要从 $i^2$ 开始筛选即可，证毕。
埃拉托斯特尼筛法的时间复杂度为 $O(n\log \log n)$ ，证明俺暂时不会（╯－＿－）╯╧╧

欧拉筛法及C++实现

基本思想

回顾上文的埃拉托斯特尼筛法，核心思想就是依次剔除素数的倍数。这里我们很容易就能想到此算法存在的缺点：对于两个素数的公倍数，我们可能会进行多次剔除。
于是便出现了欧拉筛法，核心思想便是对于一个没被筛选过的数来说，它只需要被第一个整除它的素数筛掉就行。直接上代码。

C++实现

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

void get_prime(vector<int> &prime, int n)
{
    if (n < 2)
        return;
    // 初始化一个布尔数组，用于记录每个数是否是素数
    vector<bool> is_prime(n + 1, true);
    for (int i = 2; i <= n; ++i)
    {
        // 如果是素数则将该数字加入到素数序列中
        if (is_prime[i])
            prime.push_back(i);
        // 从第一个素数开始，将它的倍数设置为非素数
        for (int j = 0; j < prime.size() && i * prime[j] <= n; ++j)
        {
            is_prime[i * prime[j]] = false;
            /**
             * 如果 i%prime[j] == 0, 则 i = prime[j] * a
             * i * prime[j+1] = prime[j] * a * prime[j+1]
             * prime[j]已经筛过i * prime[j+1]这个数了，不需要用prime[j+1]再筛选一次
             **/
            if (i % prime[j] == 0)
                break;
        }
    }
}

int main()
{

    int n = 100;
    vector<int> prime;
    get_prime(prime, n);
    for_each(prime.begin(), prime.end(), [](int &n){
        cout << n << " ";
    });
    return 0;
}

在上述代码的get_prime函数中

    if (i % prime[j] == 0)
        break;

这里就是欧拉筛法的核心所在啦，对于数字 $i$ 来说，如果 $prime[j]$ 能整除 $i$ ，则说明 $i=prime[j] \times a$ ，其中， $1<a<i$ ，而对于比 $prime[j]$ 大的素数 $prime[j+1]$ ， $i \times prime[j+1] = prime[j] \times a \times prime[j+1]$ ，也就是说 $i \times prime[j+1]$ 已经被 $prime[j]$ 筛选过了，不需要再筛选了。
欧拉筛法的时间复杂度为 $O(n)$ ，证明俺暂时也不会（╯－＿－）╯╧╧

总结

本文一共介绍了四种计算素数的方法。两种试除法，试除法一时间复杂度为 $O(n^2)$ ，试除法二时间复杂度为 $O(n \sqrt{n})$ ，两者的空间复杂度为 $O(1)$ ；两种筛法，埃拉托斯特尼筛法时间复杂度为 $O(n\log \log n)$ ，欧拉筛法时间复杂度为 $O(n)$ ，两者空间复杂度均为 $O(n)$ 。
在我们上面的实现中还存在一个问题值得我们思考，那就是如何存储素数，上面的代码中我们是用C++的vector存储素数，vector虽为动态数组，能够动态添加和删除数据，但是在大部分STL的vector的实现版本中，vector都是有预分配空间的，当预分配的空间的容量不够时，便会重新分配更大的空间，并把所有的旧数据复制到新的空间中，所以这里的开销也是难以忽视的（其它语言也存在类似的问题）。所以，我们需要知道我们应该分配多少的空间来存储求得的素数，也就是说我们得估计素数的个数以便预分配存储空间。
而如何估计素数的数量？前面我们介绍的素数定理就排上用场啦，这里就不具体说啦。

参考资料

本作品采用知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议进行许可。转载请注明: 作者staneyffer，首发于我的博客，原文链接: https://chengfy.com/post/10

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

什么是素数

素数定理

常规计算方法（试除法）及实现

基本思想

常规方法一

常规方法二

埃拉托斯特尼筛法及C++实现

基本思想

埃拉托斯特尼筛法的证明

C++实现

欧拉筛法及C++实现

基本思想

C++实现

总结

参考资料

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

数据结构和算法分析

C++