我讨论了在我的blog上生成大量素数的问题，我发现第一个十亿个素数的总和是11138479445180240497。我描述了四种不同的方法：
1.蛮力，从2开始试除法。
1.使用2，3，5，7轮生成候选项，然后使用以2、7和61为底的强伪素数测试来测试素数性;这种方法只适用于2^32，这不足以让我对前10亿个素数求和，但对你来说已经足够了。

1. Melissa奥尼尔提出的一种算法，它使用嵌入在优先级队列中的筛子，这是相当慢的。
   1.埃拉托色尼的一种分段筛，速度非常快，但需要空间来存储筛分素数和筛本身。

我已经在C++20中实现了erastosthenes的筛选，并使用多线程完成了非素数的穿孔：

#include <iostream>
#include <cstdlib>
#include <charconv>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <bit>
#include <fstream>
#include <string_view>
#include <thread>
#include <utility>
#include <new>
#include <span>
#include <array>
#include <cassert>
#if defined(_MSC_VER)
    #include <intrin.h>
#elif defined(__GNUC__) || defined(__clang__)
    #include <cpuid.h>
#endif

#if defined(_MSC_VER)
    #pragma warning(disable: 26495) // uninitialized member
#endif

using namespace std;

#if defined(__cpp_lib_hardware_interference_size)
constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
#else
constexpr ptrdiff_t CL_SIZE = 64;
#endif

size_t getL2Size();
bool smt();

int main( int argc, char **argv )
{
    if( argc < 2 )
        return EXIT_FAILURE;
    auto parse = []( char const *str, auto &value )
    {
        bool hex = str[0] == '0' && (str[1] == 'x' || str[1] == 'X');
        str += 2 * hex;
        auto [ptr, ec] = from_chars( str, str + strlen( str ), value, !hex ? 10 : 16 );
        return !(bool)ec && !*ptr;
    };
    size_t end;
    if( !parse( argv[1], end ) )
        return EXIT_FAILURE;
    if( end < 2 )
        return EXIT_FAILURE;
    if( (ptrdiff_t)end++ < 0 )
        throw bad_alloc();
    unsigned
        hc = jthread::hardware_concurrency(),
        nThreads;
    if( argc < 4 || !parse( argv[3], nThreads ) )
        nThreads = hc;
    if( !nThreads )
        return EXIT_FAILURE;
    using word_t = size_t;
    constexpr size_t
        BITS_PER_CL = CL_SIZE * 8,
        BITS = sizeof(word_t) * 8;
    size_t nBits = end / 2;
    union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE / sizeof(word_t)]; ndi_words_cl() {} };
    vector<ndi_words_cl> sieveCachelines( (nBits + BITS_PER_CL - 1) / BITS_PER_CL );
    span<word_t> sieve( &sieveCachelines[0].words[0], (nBits + BITS - 1) / BITS );
    assert(to_address( sieve.end() ) <= &to_address( sieveCachelines.end() )->words[0]);
    fill( sieve.begin(), sieve.end(), -1 );
    auto punch = [&]( size_t start, size_t end, size_t prime, auto )
    {
        size_t
            bit = start / 2,
            bitEnd = end / 2;
        if( bit >= bitEnd )
            return;
        if( prime >= BITS ) [[likely]]
            do [[likely]]
                sieve[bit / BITS] &= rotl( (word_t)-2, bit % BITS );
            while( (bit += prime) < bitEnd );
        else
        {
            auto pSieve = sieve.begin() + bit / BITS;
            do [[likely]]
            {
                word_t
                    word = *pSieve,
                    mask = rotl( (word_t)-2, bit % BITS ),
                    oldMask;
                do
                    word &= mask,
                    oldMask = mask,
                    mask = rotl( mask, prime % BITS ),
                    bit += prime;
                while( mask < oldMask );
                *pSieve++ = word;
            } while( bit < bitEnd );
        }
    };
    size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );
    for( size_t bSqrt = sqrt / 2, bit = 1; bit < bSqrt; ++bit ) [[likely]]
    {
        auto pSieve = sieve.begin() + bit / BITS;
        do [[likely]]
        {
            if( word_t w = *pSieve >> bit % BITS; w ) [[likely]]
            {
                bit += countr_zero( w );
                break;
            }
            ++pSieve;
            bit = bit + BITS & -(ptrdiff_t)BITS;
        } while( bit < bSqrt );
        if( bit >= bSqrt ) [[unlikely]]
            break;
        size_t prime = 2 * bit + 1;
        punch( prime * prime, sqrt, prime, false_type() );
    }
    auto scan = [&]( size_t start, size_t end, auto fn )
    {
        size_t
            bit = start / 2,
            bitEnd = end / 2;
        if( bit >= bitEnd )
            return;
        auto pSieve = sieve.begin() + bit / BITS;
        word_t word = *pSieve++ >> bit % BITS;
        for( ; ; )
        {
            size_t nextWordBit = bit + BITS & -(ptrdiff_t)BITS;
            for( unsigned gap; word; word >>= gap, word >>= 1 )
            {
                gap = countr_zero( word );
                bit += gap;
                if( bit >= bitEnd ) [[unlikely]]
                    return;
                fn( bit++ * 2 + 1, true_type() );
            }
            bit = nextWordBit;
            if( bit >= bitEnd ) [[unlikely]]
                return;
            word = *pSieve++;
        }
    };
    static auto dbl = []( ptrdiff_t x ) { return (double)x; };
    using range_t = pair<size_t, size_t>;
    vector<range_t> ranges;
    ranges.reserve( nThreads );
    vector<jthread> threads;
    threads.reserve( nThreads );
    auto initiate = [&]( size_t start, auto fn )
    {
        ranges.resize( 0 );
        double threadRange = dbl( end - start ) / nThreads;
        for( size_t t = nThreads, trEnd = end, trStart; t--; trEnd = trStart ) [[likely]]
        {
            trStart = start + (ptrdiff_t)((ptrdiff_t)t * threadRange + 0.5);
            size_t trClStart = trStart & -(2 * CL_SIZE * 8);
            trStart = trClStart >= start ? trClStart : start;
            if( trStart < trEnd )
                ranges.emplace_back( trStart, trEnd );
        }
        for( range_t const &range : ranges )
            if( &range != &ranges.back() )
                threads.emplace_back( fn, range.first, range.second, true_type() );
            else
                fn( range.first, range.second, false_type() );
        threads.resize( 0 );
    };
    double maxCacheRange = dbl( getL2Size() * 8 * 2 ) / 3 * (smt() && nThreads > hc / 2 ? 1 : 2);
    initiate( sqrt,
        [&]( size_t rStart, size_t rEnd, auto )
        {
            double thisThreadRange = dbl( rEnd - rStart );
            ptrdiff_t nCachePartitions = (ptrdiff_t)ceil( thisThreadRange / maxCacheRange );
            double cacheRange = thisThreadRange / dbl( nCachePartitions );
            for( size_t p = nCachePartitions, cacheEnd = rEnd, cacheStart; p--; cacheEnd = cacheStart ) [[likely]]
            {
                cacheStart = rStart + (ptrdiff_t)((double)(ptrdiff_t)p * cacheRange + 0.5);
                cacheStart &= -(2 * CL_SIZE * 8);
                cacheStart = cacheStart >= sqrt ? cacheStart : sqrt;
                scan( 3, sqrt,
                    [&]( size_t prime, auto )
                    {
                        size_t start = (cacheStart + prime - 1) / prime * prime;
                        start = start % 2 ? start : start + prime;
                        punch( start, cacheEnd, prime, true_type() );
                    } );
            }
        } );
    if( bool count = false; argc >= 3 && (!*argv[2] || (count = strcmp( argv[2], "*" ) == 0)) )
    {
        if( !count )
            return EXIT_SUCCESS;
        atomic<size_t> aNPrimes( 1 );
        initiate( 3,
            [&]( size_t rStart, size_t rEnd, auto )
            {
                size_t nPrimes = 0;
                scan( rStart, rEnd, [&]( ... ) { ++nPrimes; } );
                aNPrimes.fetch_add( nPrimes, memory_order::relaxed );
            } );
        cout << aNPrimes.load( memory_order::relaxed ) << " primes" << endl;
        return EXIT_SUCCESS;
    }
    ofstream ofs;
    ofs.exceptions( ofstream::failbit | ofstream::badbit );
    ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary | ofstream::trunc );
    constexpr size_t
        BUF_SIZE = 0x100000,
        AHEAD = 32;
    union ndi_char { char c; ndi_char() {} };
    vector<ndi_char> rawPrintBuf( BUF_SIZE + AHEAD - 1 );
    span printBuf( &rawPrintBuf.front().c, &rawPrintBuf.back().c + 1 );
    auto
        bufBegin = printBuf.begin(),
        bufEnd = bufBegin,
        bufFlush = bufBegin + BUF_SIZE;
    auto print = [&]() { ofs << string_view( bufBegin, bufEnd ); };
    auto printPrime = [&]( size_t prime, auto )
    {
        auto [ptr, ec] = to_chars( &*bufEnd, &bufEnd[AHEAD - 1], prime );
        if( (bool)ec ) [[unlikely]]
            throw system_error( (int)ec, generic_category(), "converson failed" );
        bufEnd = ptr - &*bufBegin + bufBegin; // pointer to iterator - NOP
        *bufEnd++ = '\n';
        if( bufEnd >= bufFlush ) [[unlikely]]
            print(), bufEnd = bufBegin;
    };
    printPrime( 2, false_type() );
    scan( 3, end, printPrime );
    print();
}

array<unsigned, 4> cpuId( unsigned code )
{
    array<unsigned, 4> regs;
#if defined(_MSC_VER)
    __cpuid( (int *)&regs[0], code );
#elif defined(__GNUC__) || defined(__clang__)
    __cpuid(code, regs[0], regs[1], regs[2], regs[3]);
#endif
    return regs;
}

bool smt()
{
    if( cpuId( 0 )[0] < 1 )
        return false;
    return cpuId( 1 )[3] >> 28 & 1;
}

size_t getL2Size()
{
    if( cpuId( 0x80000000u )[0] < 0x80000006u )
        return 512ull * 1024;
    return (cpuId( 0x80000006u )[2] >> 16) * 1024;
}

字符串
几乎所有的CPU时间都花在了多线程代码中，用于对非质数进行打孔。在打孔线程中，位范围被分区以适应L2缓存。如果线程数量超过SMT系统上逻辑核数量的一半，则线程子分区的大小为一半。
我只是增加了只存储奇素数的功能，即不存储2的倍数，从而提高了大约87%的性能。
在我的Ryzen 9 7950 X Zen 4 16核CPU上，产生所有素数高达2^32，这几乎是210 E6素数，在所有核心上大约需要0.15秒，抑制文件输出;在一个核心上，算法大约需要1.73秒。程序的第三个参数是线程数。在我的16个核心Zen 4上有16个线程，在一个线程上，我得到了70%的扩展。占用更多的SMT兄弟线程几乎不会带来进一步的加速，因为代码依赖于L2缓存的吞吐量。
文件输出是通过to_chars和我自己的缓冲来完成的，以加快I/O。在我的计算机（64 GB，PCIe 4.0 SSD）上，从上述范围生成2.2GB的文件大约需要额外的两秒钟。

这可能会加快一点：

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

int main()
{
    std::vector<unsigned long long> numbers;
    unsigned long long maximum = 4294967296;
    for (unsigned long long i = 2; i <= maximum; ++i)
    {
        if (numbers.empty())
        {
            numbers.push_back(i);
            continue;
        }

        if (std::none_of(numbers.begin(), numbers.end(), [&](unsigned long long p)
        {
            return i % p == 0;
        }))
        {
            numbers.push_back(i);
        }

    }

    std::cout << "Primes:  " << std::endl;
    std::copy(numbers.begin(), numbers.end(), std::ostream_iterator<int>(std::cout, " "));

    return 0;
}

字符串
它有点像埃拉托塞尼筛的逆（不是从极限以下的每个数字开始并消除倍数，而是从2开始并忽略极限以下的倍数）。

最快的方法可能是使用预先生成的列表。
http://www.bigprimes.net/有前14亿个质数可供下载，其中应该包括300亿以下的所有质数。
我想加载二进制文件可能需要太长的时间，当它是几个千兆字节的大小。

你有没有衡量过什么是最耗时的？是筛子本身，还是输出的写作？
一个快速的方法来加速筛子是停止担心所有的偶数。只有一个偶数是素数，你可以硬编码它。这将把你的数组大小减少一半，这将极大地帮助你克服物理内存的限制。

vector<bool> sieve((limit+1)/2, false);
...
  for(long long m = n*n/2; m<=limit/2; m=m+n)
    sieve[m] = true;

字符串
至于输出本身，cout是出了名的低效。自己调用itoa或其他等效函数，然后使用cout.write输出可能会更高效。您甚至可以采用传统方法，将fwrite与stdout一起使用。

我用C写了一个快速的方法，在我的Ryzen 9 3900x上使用一个线程在三分钟内输出高达40亿个素数。如果你把它输出到一个文件，它最终是2.298GB，我认为它使用了大约20GB的内存来完成。

#include <stdlib.h>
#include <stdio.h>

#define ARRSIZE 4000000000
#define MAXCALC ARRSIZE/2

int main() {
    long int z;
    long int *arr = (long int*) malloc((ARRSIZE) * sizeof(long int));
    for (long int x=3;x <= MAXCALC; x++) {
        if (x % 10 == 3 || x % 10 == 7 || x % 10 == 9) {
            for (long int y=3; y < MAXCALC; y++){
                    z = x * y;
                    if (z < ARRSIZE)
                        arr[z] = 1;
                    else
                        break;
            }
        }
    }
    printf("2 3 5 ");
    for (long int x=7; x < ARRSIZE; x++) {
        if (x % 2 != 0 && x % 10 != 5)
            if (arr[x] != 1)
                printf("%ld ", x);
    }
    printf("\n");

    free(arr);
    return EXIT_SUCCESS;
}

字符串

我写了代码i python输出所有素数小于10亿在8.7秒.但我不确定如果你的第一个40亿素数或艾伦素数小于.无论如何，这是我的解决方案：

import numpy as np
from math import isqrt

def all_primes_less_than(n):
    is_prime = np.full(n,True)
    is_prime[:2] = False
    is_prime[4::2] = False
    for p in range(3,isqrt(n),2):
        if is_prime[p]: is_prime[p*p::p] = False
    return is_prime.nonzero()[0]

字符串