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

什么是素数

素数又称质数，是指大于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

最后编辑于：2018.11.16 00:26:02

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 206,126评论 6赞 481
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 88,254评论 2赞 382
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 152,445评论 0赞 341
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 55,185评论 1赞 278
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 64,178评论 5赞 371
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,970评论 1赞 284
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,276评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,927评论 0赞 259
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 43,400评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,883评论 2赞 323
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,997评论 1赞 333
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,646评论 4赞 322
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,213评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,204评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,423评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,423评论 2赞 352
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,722评论 2赞 345

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

什么是素数

素数定理

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

基本思想

常规方法一

常规方法二

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

基本思想

埃拉托斯特尼筛法的证明

C++实现

欧拉筛法及C++实现

基本思想

C++实现

总结

参考资料

推荐阅读更多精彩内容