找到 40 亿以下所有素数的最快方法

问题描述 投票:0回答:7

我正在尝试打印 2**32 以下的每个质数。现在我正在使用布尔向量构建一个筛子,然后在制作筛子后打印出素数。仅打印 10 亿以内的素数就需要 4 分钟。有没有更快的方法来做到这一点?这是我的代码

#include <iostream>
#include <cstdlib>
#include <vector>
#include <math.h>

using namespace std;

int main(int argc, char **argv){
  long long limit = atoll(argv[1]);
  //cin >> limit;
  long long sqrtlimit = sqrt(limit);

  vector<bool> sieve(limit+1, false);

  for(long long n = 4; n <= limit; n += 2)
    sieve[n] = true;

  for(long long n=3; n <= sqrtlimit; n = n+2){
    if(!sieve[n]){
      for(long long m = n*n; m<=limit; m=m+(2*n))
        sieve[m] = true;
    }
  }

  long long last;
  for(long long i=limit; i >= 0; i--){
    if(sieve[i] == false){
      last = i;
      break;
    }
  }
  cout << last << endl;

  for(long long i=2;i<=limit;i++)
  {
    if(!sieve[i])
      if(i != last)
        cout<<i<<",";
      else
        cout<<i;
  }
  cout<<endl;
c++ algorithm primes sieve-of-eratosthenes
7个回答
4
投票

我在我的博客中讨论了生成大量素数的问题,我发现前十亿个素数的总和是11138479445180240497。我描述了四种不同的方法:

  1. 暴力破解,使用试除法测试从 2 开始的每个数字。

  2. 使用 2、3、5、7 轮生成候选项,然后使用对基数 2、7 和 61 的强伪素数测试来测试素数;这个方法最多只能计算 2^32,这不足以让我对前十亿个素数求和,但对你来说已经足够了。

  3. Melissa O'Neill 提出的一种算法,使用嵌入优先级队列的筛子,速度相当慢。

  4. 埃拉托色尼分段筛,速度非常快,但需要空间来存储筛分素数和筛子本身。


0
投票

这可能会加快速度:

#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 开始并忽略直到限制的倍数)。


0
投票

最快的方法可能是获取预先生成的列表。

http://www.bigprimes.net/ 有前 14 亿个素数可供下载,其中应包括 300 亿左右以下的每个素数。

我认为当二进制文件大小为几 GB 时,加载二进制文件可能需要很长时间。


0
投票

您是否对最耗时的事情进行了基准测试?是筛子本身,还是输出的书写?

加快筛子速度的一种快速方法是不再担心所有偶数。只有一个偶数是素数,您可以对其进行硬编码。这会将数组的大小减少一半,如果您遇到物理内存的限制,这将非常有帮助。

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
一起使用。


0
投票

我用 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;
}

0
投票

我用 python 编写了代码,可以在 8.7 秒内输出所有小于 10 亿的素数。但我不确定你是指前 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]

0
投票

我在 C++20 中使用 countr_zero 实现了埃拉斯托斯特尼筛,以搜索位图中的下一个素数:

#include <iostream>
#include <vector>
#include <charconv>
#include <string_view>
#include <algorithm>
#include <fstream>
#include <cctype>
#include <cstring>
#include <bit>
#include <cmath>

using namespace std;

int main( int argc, char **argv )
{
    constexpr size_t BUF_SIZE = 0x100000;
    if( argc < 2 )
        return EXIT_FAILURE;
    char const *sizeStr = argv[1];
    bool hex = sizeStr[0] == '0' && (sizeStr[1] == 'x' || sizeStr[1] == 'X');
    sizeStr += 2 * hex;
    size_t end;
    if( from_chars_result fcr = from_chars( sizeStr, sizeStr + strlen( sizeStr ), end, !hex ? 10 : 16 ); fcr.ec != errc()  || *fcr.ptr )
        return EXIT_FAILURE;
    using word_t = size_t;
    constexpr size_t BITS = sizeof(word_t) * 8;
    if( (ptrdiff_t)++end < 0 )
        throw bad_alloc();
    vector<word_t> sieve( (end + BITS - 1) / BITS, (word_t)0xAAAAAAAAAAAAAAAAu );
    ofstream ofs;
    constexpr size_t AHEAD = 32;
    union ndi_char { char c; ndi_char() {} };
    vector<ndi_char> printBuf;
    vector<ndi_char>::iterator bufBegin, bufEnd, bufFlush;
    if( argc < 3 || *argv[2] )
    {
        ofs.exceptions( ofstream::failbit | ofstream::badbit );
        ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );
        printBuf.resize( BUF_SIZE + AHEAD - 1 );
        bufBegin = printBuf.begin();
        bufEnd = bufBegin;
        bufFlush = bufBegin + BUF_SIZE;
    }
    auto print = [&]()
    {
        ofs << string_view( &bufBegin->c, &to_address( bufEnd )->c );
        bufEnd = bufBegin;
    };
    auto const pSieveEnd = sieve.end();
    size_t prime = 2;
    auto primeTime = [&]( size_t end, auto punch )
    {
        if( prime >= end )
            return;
        for( ; ; )
        {
            if( ofs.is_open() )
            {
                auto [ptr, ec] = to_chars( &bufEnd->c, &bufEnd[AHEAD - 1].c, prime );
                if( ec != errc() ) [[unlikely]]
                    throw system_error( (int)ec, generic_category(), "converson failed" );
                bufEnd = ptr - &bufBegin->c + bufBegin; // NOP
                bufEnd++->c = '\n';
                if( bufEnd >= bufFlush ) [[unlikely]]
                    print();
            }
            if( ++prime >= end ) [[unlikely]]
                return;
            for( auto pSieve = sieve.begin() + prime / BITS; ; )
            {
                if( word_t w = *pSieve >> prime % BITS; w ) [[likely]]
                {
                    prime += countr_zero( w );
                    break;
                }
                prime = prime + BITS & -(ptrdiff_t)BITS;
                if( ++pSieve == pSieveEnd ) [[unlikely]]
                    return;
            }
            if( prime >= end ) [[unlikely]]
                return;
            punch();
        }
    };
    primeTime( (ptrdiff_t)ceil( sqrt( (ptrdiff_t)end ) ),
        [&]()
        {
            auto pSieve = sieve.begin();
            size_t multiple = prime * prime;
            if( prime >= BITS ) [[likely]]
                do
                    pSieve[multiple / BITS] &= rotl( (word_t)-2, multiple % BITS ),
                    multiple += prime;
                while( multiple < end );
            else
            {
                word_t mask = rotl( (word_t)-2, multiple % BITS );
                pSieve += multiple / BITS;
                do
                {
                    word_t
                        word = *pSieve,
                        oldMask;
                    do
                        word &= mask,
                        oldMask = mask,
                        mask = rotl( mask, prime % BITS );
                    while( mask < oldMask );
                    *pSieve = word;
                } while( ++pSieve != pSieveEnd );
            }
        } );
    primeTime( end, [] {} );
    if( ofs.is_open() )
        print();
    return EXIT_SUCCESS;
}

在我的带有 g++ 12 的 AMD Zen4-CPU 上,代码运行 17.7 秒来生成最多 2^32 的所有素数,大约是 200E6 个素数,这会产生一个 2.4 GB 的文件。我不使用直接流 I/O,而是使用我自己的缓冲来加速输出。
该代码使用带有 primeTime lambda 的分段方法,其调用运算符被实例化两次,一次用于

ceil( sqrt( max ) )
以下的所有素数,一次用于此后的所有素数。使用 g++ 和 clang++ 函数,只调用一次,并在调用它们的地方实例化,从而产生非常高效的代码。

© www.soinside.com 2019 - 2024. All rights reserved.