我正在尝试打印 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;
我在我的博客中讨论了生成大量素数的问题,我发现前十亿个素数的总和是11138479445180240497。我描述了四种不同的方法:
暴力破解,使用试除法测试从 2 开始的每个数字。
使用 2、3、5、7 轮生成候选项,然后使用对基数 2、7 和 61 的强伪素数测试来测试素数;这个方法最多只能计算 2^32,这不足以让我对前十亿个素数求和,但对你来说已经足够了。
Melissa O'Neill 提出的一种算法,使用嵌入优先级队列的筛子,速度相当慢。
埃拉托色尼分段筛,速度非常快,但需要空间来存储筛分素数和筛子本身。
这可能会加快速度:
#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 亿左右以下的每个素数。
我认为当二进制文件大小为几 GB 时,加载二进制文件可能需要很长时间。
您是否对最耗时的事情进行了基准测试?是筛子本身,还是输出的书写?
加快筛子速度的一种快速方法是不再担心所有偶数。只有一个偶数是素数,您可以对其进行硬编码。这会将数组的大小减少一半,如果您遇到物理内存的限制,这将非常有帮助。
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;
}
我用 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]
我在 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++ 函数,只调用一次,并在调用它们的地方实例化,从而产生非常高效的代码。