素数的计算: 从试除到筛法(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
  1. 去掉素数2的所有倍数:
  • 2 3 5 7 9 11 13 15 17 19 21 23 25
  1. 去掉素数3的所有倍数:
  • 2 3 5 7 11 13 17 19 23 25
  1. 去掉素数5的所有倍数:
  • 2 3 7 11 17 19 23
  1. 因为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开始筛选的,且剩余序列中从2i^2-1的所有数都是素数。证明如下:
如果i=2i^2之前的序列为: 2,3,都是素数。
如果i>2,因为1<i^2-1 < i^2,所以i-1<\sqrt{i^2-1} < ii-1i的前一个数,即\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都是有预分配空间的,当预分配的空间的容量不够时,便会重新分配更大的空间,并把所有的旧数据复制到新的空间中,所以这里的开销也是难以忽视的(其它语言也存在类似的问题)。所以,我们需要知道我们应该分配多少的空间来存储求得的素数,也就是说我们得估计素数的个数以便预分配存储空间。
而如何估计素数的数量?前面我们介绍的素数定理就排上用场啦,这里就不具体说啦。

参考资料

  1. 质素 - 维基百科,自由的百科全书
  2. 算术基本定理 - 维基百科,自由的百科全书
  3. 素数计数函数 - 维基百科,自由的百科全书
  4. Sieve of Eratosthenes - Wikipedia, the free encyclopedia

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,688评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,559评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,749评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,581评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,741评论 3 271
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,684评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,122评论 2 292
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,847评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,441评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,939评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,333评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,783评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,275评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,830评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,444评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,553评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,618评论 2 249

推荐阅读更多精彩内容

  • 漆黑的夜晚好安静 太阳已经落进大海里 月亮好像一个香蕉啊 树叶在唱着儿歌 星星为她伴舞 (四周大的女儿在某一天和妈...
    看客不说话阅读 139评论 2 1
  • 第一组 第一组 第二组 第二组 第三组 第三组 第四组 第四组 第五组 第五组 第六组 第六组
    246任晶璐阅读 182评论 0 0
  • 我不勇敢,但是还是要飞 生活就像一根刺,让你痛不欲生,让你偶尔不痛不痒,但却始终都拔不出来的刺。 从大一刚入学到现...
    昙若惊鸿阅读 311评论 1 1
  • 文/三月鱼 1. 前几天,和一个许久不见的女友聊天,她说最近比较烦。 问及原因,原来是她妹妹正在闹离婚,妹妹有事,...
    653e0adfb5bf阅读 1,009评论 5 16