快速查找计算9百万万亿整数内全部素数（质数）的C++代码

用素数筛查找小于等于某个给定整数的全部素数，是一种较为高效的方法，具体的原理网上很多，这里就不赘述了。

但即使都是运用的素数筛原理，不同的算法设计，也可以带来巨大的效率差异。最近从网上搜索学习了相关的原理，验证相关的算法，顺便测试、完善了相关的C++代码，在这里贴出来，分享给有兴趣的读者。

以下两段代码，都可以计算查找9223372036854775807（9百万万亿）内的所有素数。

第一段代码的效率中等，但算法简单易理解。第二段代码据说是目前地球上发现的最高效的方法，但逻辑稍复杂，理解起来有些难度。

这两段代码我均在gcc version 9.2.0+unbutun 20.04 linux上测试通过。

使用第一种算法，在一台2核（CPU主频1.5GHz）+8GB内存的虚拟机(最小化模式安装unbutun 20.04操作系统)上计算查找十亿以内的所有素数，耗时15秒；在同样的硬件水平下，第二种算法耗时80毫秒，显然效率高很多。用第二种算法查找计算千亿范围内的全部素数，耗时4024毫秒；查找计算万亿范围内的全部素数，则耗时毫秒39秒，效率确实比第一种算法高很多。

当计算大量的素数时，需要大量的内存来保存所找到的素数。第一种算法使用动态数组来保存搜索到的素数，第二种算法在这方面的处理更为精准，细化。注意，如果你的计算机内存不大，建议不要尝试去计算超过万亿范围的素数。

第一种算法代码：

/*---------------------------------------------------------------------------------------
 功能：查找输出所有小于指定值的素数,上限是输出小于等于9223372036854775807的全部质数，并输出保存到CSV文件
---------------------------------------------------------------------------------------*/
#include <bits/stdc++.h>
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

//long long型对应的最大整数为9223372036854775807
#define ll long long

//CalPrime函数计算列出小于等于MaxInt的所有素数，并把找到的素数保存到PrimeList数组
void CalPrime(ll MaxInt,ll &CntPrime,ll PrimeList[])
{
    //MaxInt 要查找素数的范围上限;
    //CntPrime 小于等于MaxInt的素数个数（不包含1）;
    //PrimeList 找到的素数存储到这个数组列表;
    ll i,j;
    clock_t start, end;

    bool *valid=new bool[MaxInt];
    CntPrime=0;

    //开始计算，获取当前系统时间，以便计算结束时统计耗时
    start = clock();

    for(i=2; i<=MaxInt; i++)    //将要查找素数范围内的valid全部赋值为true;
        valid[i]=true;
    for(i=2; i<=MaxInt; i++)
    {
        if(valid[i])
        {
            if(MaxInt/i<i)
                break;
            for(j=i*i; j<=MaxInt; j+=i)    //将不是素数的数所对应的valid[i]赋值为false;
                valid[j]=false;
        }
    }
    //将[2至MaxInt]范围内的所有素数存入PrimeList数组中;
    for(i=2; i<=MaxInt; i++)
    {
        if(valid[i])
            PrimeList[CntPrime++]=i;
    }
    delete valid;   //释放bool数组，回收内存
    //显示计算所用的时间,我计算999999999(9亿)内的素数，发现有50847534个,计算用时[15391.218000]毫秒
    //这个算法和目前最快的算法比起来不算突出，我用同样的机器，用另一种更快但也更复杂的算法计算查找999999999(9亿)内的
    //素数耗时是80.968000毫秒
    end = clock() - start;
    printf("%lld(%lld亿)内的素数个数为%lld,计算用时[%lf]毫秒\n",MaxInt,MaxInt/100000000,CntPrime,(double)end/1000);
}
//---------------------------------------------------------------------------------------
int main()
{
    ll  MaxInt,CntPrime;
    ll  PrimeListLen;
    int tmpval;
    CntPrime=0;
    //输入范围，即查找所有小于n的素数
    printf("\n请输入需计算素数的范围上限，输入0则结束运行：");
    cin>>MaxInt;
    if(MaxInt<1) return 0;

    PrimeListLen=MaxInt/2+1; //其实可以用Pi(x)计算更准确的个数，以节约运行所需内存，这里懒得写了
    //存放找到素数的动态数组
    ll *PrimeList =new ll[PrimeListLen];
    //调用getPrime函数搜索小于n的所有素数
    CalPrime(MaxInt,CntPrime,PrimeList);
    //打印输出找到的质数个数，提问是否保存到本地文件
    //999999999内的素数存到CSV文件，大小是478MB
    printf("\n需把小于等于[%lld]的全部[%lld]个素数保存到CSV文件吗？（不保存请输入0）",MaxInt,CntPrime);
    char filename[128];
    cin>>tmpval;
    if(tmpval<1){
        delete PrimeList;
        return 0;
    }
    sprintf(filename,"./prime-in-%lld.csv",MaxInt);
    printf("\n下面开始把所有合计[%lld]个小于等于[%lld]的质数写入文件【%s】.....",CntPrime,MaxInt,filename);
    ofstream out(filename,ios::app);//app表示每次操作前均定位到文件末尾
    if(out.fail()){
        cout<<"error\n";
    }
    for(ll i=0;i<CntPrime;i++)
        out<<PrimeList[i]<<",";
    out.close();
    printf("\n小于【%lld】的[%lld]个质数已全部写入到文件[%s]\n",MaxInt,CntPrime,filename);
    delete PrimeList;   //删除保存素数的数组，回收内存
    return 0;
}
//---------------------------------------------------------------------------------------

下图是第一种算法运行的效果：

我尝试把10亿以内的所有5084千万多个素数写到CSV文件试了一下，文件的大小是478MB。

第二种算法的代码：

/*
 功能：查找所有小于等于指定整数的素数,上限是输出小于等于9223372036854775807的全部质数
 这是目前能找到的，最高效的查找素数算法
 */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <bits/stdc++.h>
#include <iostream>
#include <vector>
using namespace std;

#define __int64 long long
__int64 *primarr, *v;
__int64 q = 1, p = 1;

//π(n)
__int64 pi(__int64 n, __int64 primarr[], __int64 len)
{
    __int64 i = 0, mark = 0;
    for (i = len - 1; i > 0; i--) {
        if (primarr[i] < n) {
            mark = 1;
            break;
        }
    }
    if (mark)
        return i + 1;
    return 0;
}

//Φ(x,a)
__int64 phi(__int64 x, __int64 a, __int64 m)
{
    if (a == m)
        return (x / q) * p + v[x % q];
    if (x < primarr[a - 1])
        return 1;
    return phi(x, a - 1, m) - phi(x / primarr[a - 1], a - 1, m);
}

__int64 prime(__int64 n)
{
    char *mark;
    __int64 mark_len;
    __int64 count = 0;
    __int64 i, j, m = 7;
    __int64 sum = 0, s = 0;
    __int64 len, len2, len3;

    mark_len = (n < 10000) ? 10002 : ((__int64)exp(2.0 / 3 * log(n)) + 1);

    //筛选n^(2/3)或n内的素数
    mark = (char *)malloc(sizeof(char) * mark_len);
    memset(mark, 0, sizeof(char) * mark_len);
    for (i = 2; i < (__int64)sqrt(mark_len); i++) {
        if (mark[i])
            continue;
        for (j = i + i; j < mark_len; j += i)
            mark[j] = 1;
    }
    mark[0] = mark[1] = 1;

    //统计素数数目
    for (i = 0; i < mark_len; i++)
        if (!mark[i])
            count++;

    //保存素数
    primarr = (__int64 *)malloc(sizeof(__int64) * count);
    j = 0;
    for (i = 0; i < mark_len; i++)
        if (!mark[i])
            primarr[j++] = i;

    if (n < 10000)
        return pi(n, primarr, count);

    //n^(1/3)内的素数数目
    len = pi((__int64)exp(1.0 / 3 * log(n)), primarr, count);
    //n^(1/2)内的素数数目
    len2 = pi((__int64)sqrt(n), primarr, count);
    //n^（2/3)内的素数数目
    len3 = pi(mark_len - 1, primarr, count);

    //乘积个数
    j = mark_len - 2;
    for (i = (__int64)exp(1.0 / 3 * log(n)); i <= (__int64)sqrt(n); i++) {
        if (!mark[i]) {
            while (i * j > n) {
                if (!mark[j])
                    s++;
                j--;
            }
            sum += s;
        }
    }
    free(mark);
    sum = (len2 - len) * len3 - sum;
    sum += (len * (len - 1) - len2 * (len2 - 1)) / 2;

    //欧拉函数
    if (m > len)
        m = len;
    for (i = 0; i < m; i++) {
        q *= primarr[i];
        p *= primarr[i] - 1;
    }
    v = (__int64 *)malloc(sizeof(__int64) * q);
    for (i = 0; i < q; i++)
        v[i] = i;
    for (i = 0; i < m; i++)
        for (j = q - 1; j >= 0; j--)
            v[j] -= v[j / primarr[i]];

    sum = phi(n, len, m) - sum + len - 1;
    free(primarr);
    free(v);
    return sum;
}
//---------------------------------------------------------------------
int main()
{
    __int64 n;
    __int64 count;
    int h;
    clock_t start, end;
    n=10000;

    while(n>1)
    {
        printf("\n请输入需计算的整数上限，输入0则结束运行：");
        //输入范围，即查找所有小于n的素数
        cin>>n;
        if(n==0)  break;
        p=1;
        q=1;
        start = clock();
        count = prime(n);
        end = clock() - start;
        printf("%lld(%lld亿)内的素数个数为%lld\n",n,n/100000000,count);
        printf("用时%lf毫秒\n",(double)end/1000);
    }
    return 0;
}
//---------------------------------------------------------------------

下图是第二种算法计算10亿以内所有素数的输出截图：

计算万亿范围内的所有素数，耗时39秒。

（如果觉得有帮助，请点个赞鼓励下作者吧^_^）