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

dpiehjr4  于 2024-01-09  发布在  其他
关注(0)|答案(7)|浏览(173)

我正在尝试打印2**32以下的所有素数。现在我正在使用bool vector构建一个筛子，然后打印出筛子后的素数。仅打印出10亿个素数就需要4分钟。有没有更快的方法来做到这一点？？这是我的代码

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <vector>
  4. #include <math.h>
  6. using namespace std;
  8. int main(int argc, char **argv){
  9.   long long limit = atoll(argv[1]);
  10.   //cin >> limit;
  11.   long long sqrtlimit = sqrt(limit);
  13.   vector<bool> sieve(limit+1, false);
  15.   for(long long n = 4; n <= limit; n += 2)
  16.     sieve[n] = true;
  18.   for(long long n=3; n <= sqrtlimit; n = n+2){
  19.     if(!sieve[n]){
  20.       for(long long m = n*n; m<=limit; m=m+(2*n))
  21.         sieve[m] = true;
  22.     }
  23.   }
  25.   long long last;
  26.   for(long long i=limit; i >= 0; i--){
  27.     if(sieve[i] == false){
  28.       last = i;
  29.       break;
  30.     }
  31.   }
  32.   cout << last << endl;
  34.   for(long long i=2;i<=limit;i++)
  35.   {
  36.     if(!sieve[i])
  37.       if(i != last)
  38.         cout<<i<<",";
  39.       else
  40.         cout<<i;
  41.   }
  42.   cout<<endl;

字符串

c++

7条答案

按热度按时间
cuxqih21

cuxqih211#

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

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

disho6za2#

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

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <charconv>
  4. #include <cstring>
  5. #include <vector>
  6. #include <algorithm>
  7. #include <cmath>
  8. #include <bit>
  9. #include <fstream>
  10. #include <string_view>
  11. #include <thread>
  12. #include <utility>
  13. #include <new>
  14. #include <span>
  15. #include <array>
  16. #include <cassert>
  17. #if defined(_MSC_VER)
  18.     #include <intrin.h>
  19. #elif defined(__GNUC__) || defined(__clang__)
  20.     #include <cpuid.h>
  21. #endif
  23. #if defined(_MSC_VER)
  24.     #pragma warning(disable: 26495) // uninitialized member
  25. #endif
  27. using namespace std;
  29. #if defined(__cpp_lib_hardware_interference_size)
  30. constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
  31. #else
  32. constexpr ptrdiff_t CL_SIZE = 64;
  33. #endif
  35. size_t getL2Size();
  36. bool smt();
  38. int main( int argc, char **argv )
  39. {
  40.     if( argc < 2 )
  41.         return EXIT_FAILURE;
  42.     auto parse = []( char const *str, auto &value )
  43.     {
  44.         bool hex = str[0] == '0' && (str[1] == 'x' || str[1] == 'X');
  45.         str += 2 * hex;
  46.         auto [ptr, ec] = from_chars( str, str + strlen( str ), value, !hex ? 10 : 16 );
  47.         return !(bool)ec && !*ptr;
  48.     };
  49.     size_t end;
  50.     if( !parse( argv[1], end ) )
  51.         return EXIT_FAILURE;
  52.     if( end < 2 )
  53.         return EXIT_FAILURE;
  54.     if( (ptrdiff_t)end++ < 0 )
  55.         throw bad_alloc();
  56.     unsigned
  57.         hc = jthread::hardware_concurrency(),
  58.         nThreads;
  59.     if( argc < 4 || !parse( argv[3], nThreads ) )
  60.         nThreads = hc;
  61.     if( !nThreads )
  62.         return EXIT_FAILURE;
  63.     using word_t = size_t;
  64.     constexpr size_t
  65.         BITS_PER_CL = CL_SIZE * 8,
  66.         BITS = sizeof(word_t) * 8;
  67.     size_t nBits = end / 2;
  68.     union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE / sizeof(word_t)]; ndi_words_cl() {} };
  69.     vector<ndi_words_cl> sieveCachelines( (nBits + BITS_PER_CL - 1) / BITS_PER_CL );
  70.     span<word_t> sieve( &sieveCachelines[0].words[0], (nBits + BITS - 1) / BITS );
  71.     assert(to_address( sieve.end() ) <= &to_address( sieveCachelines.end() )->words[0]);
  72.     fill( sieve.begin(), sieve.end(), -1 );
  73.     auto punch = [&]( size_t start, size_t end, size_t prime, auto )
  74.     {
  75.         size_t
  76.             bit = start / 2,
  77.             bitEnd = end / 2;
  78.         if( bit >= bitEnd )
  79.             return;
  80.         if( prime >= BITS ) [[likely]]
  81.             do [[likely]]
  82.                 sieve[bit / BITS] &= rotl( (word_t)-2, bit % BITS );
  83.             while( (bit += prime) < bitEnd );
  84.         else
  85.         {
  86.             auto pSieve = sieve.begin() + bit / BITS;
  87.             do [[likely]]
  88.             {
  89.                 word_t
  90.                     word = *pSieve,
  91.                     mask = rotl( (word_t)-2, bit % BITS ),
  92.                     oldMask;
  93.                 do
  94.                     word &= mask,
  95.                     oldMask = mask,
  96.                     mask = rotl( mask, prime % BITS ),
  97.                     bit += prime;
  98.                 while( mask < oldMask );
  99.                 *pSieve++ = word;
  100.             } while( bit < bitEnd );
  101.         }
  102.     };
  103.     size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );
  104.     for( size_t bSqrt = sqrt / 2, bit = 1; bit < bSqrt; ++bit ) [[likely]]
  105.     {
  106.         auto pSieve = sieve.begin() + bit / BITS;
  107.         do [[likely]]
  108.         {
  109.             if( word_t w = *pSieve >> bit % BITS; w ) [[likely]]
  110.             {
  111.                 bit += countr_zero( w );
  112.                 break;
  113.             }
  114.             ++pSieve;
  115.             bit = bit + BITS & -(ptrdiff_t)BITS;
  116.         } while( bit < bSqrt );
  117.         if( bit >= bSqrt ) [[unlikely]]
  118.             break;
  119.         size_t prime = 2 * bit + 1;
  120.         punch( prime * prime, sqrt, prime, false_type() );
  121.     }
  122.     auto scan = [&]( size_t start, size_t end, auto fn )
  123.     {
  124.         size_t
  125.             bit = start / 2,
  126.             bitEnd = end / 2;
  127.         if( bit >= bitEnd )
  128.             return;
  129.         auto pSieve = sieve.begin() + bit / BITS;
  130.         word_t word = *pSieve++ >> bit % BITS;
  131.         for( ; ; )
  132.         {
  133.             size_t nextWordBit = bit + BITS & -(ptrdiff_t)BITS;
  134.             for( unsigned gap; word; word >>= gap, word >>= 1 )
  135.             {
  136.                 gap = countr_zero( word );
  137.                 bit += gap;
  138.                 if( bit >= bitEnd ) [[unlikely]]
  139.                     return;
  140.                 fn( bit++ * 2 + 1, true_type() );
  141.             }
  142.             bit = nextWordBit;
  143.             if( bit >= bitEnd ) [[unlikely]]
  144.                 return;
  145.             word = *pSieve++;
  146.         }
  147.     };
  148.     static auto dbl = []( ptrdiff_t x ) { return (double)x; };
  149.     using range_t = pair<size_t, size_t>;
  150.     vector<range_t> ranges;
  151.     ranges.reserve( nThreads );
  152.     vector<jthread> threads;
  153.     threads.reserve( nThreads );
  154.     auto initiate = [&]( size_t start, auto fn )
  155.     {
  156.         ranges.resize( 0 );
  157.         double threadRange = dbl( end - start ) / nThreads;
  158.         for( size_t t = nThreads, trEnd = end, trStart; t--; trEnd = trStart ) [[likely]]
  159.         {
  160.             trStart = start + (ptrdiff_t)((ptrdiff_t)t * threadRange + 0.5);
  161.             size_t trClStart = trStart & -(2 * CL_SIZE * 8);
  162.             trStart = trClStart >= start ? trClStart : start;
  163.             if( trStart < trEnd )
  164.                 ranges.emplace_back( trStart, trEnd );
  165.         }
  166.         for( range_t const &range : ranges )
  167.             if( &range != &ranges.back() )
  168.                 threads.emplace_back( fn, range.first, range.second, true_type() );
  169.             else
  170.                 fn( range.first, range.second, false_type() );
  171.         threads.resize( 0 );
  172.     };
  173.     double maxCacheRange = dbl( getL2Size() * 8 * 2 ) / 3 * (smt() && nThreads > hc / 2 ? 1 : 2);
  174.     initiate( sqrt,
  175.         [&]( size_t rStart, size_t rEnd, auto )
  176.         {
  177.             double thisThreadRange = dbl( rEnd - rStart );
  178.             ptrdiff_t nCachePartitions = (ptrdiff_t)ceil( thisThreadRange / maxCacheRange );
  179.             double cacheRange = thisThreadRange / dbl( nCachePartitions );
  180.             for( size_t p = nCachePartitions, cacheEnd = rEnd, cacheStart; p--; cacheEnd = cacheStart ) [[likely]]
  181.             {
  182.                 cacheStart = rStart + (ptrdiff_t)((double)(ptrdiff_t)p * cacheRange + 0.5);
  183.                 cacheStart &= -(2 * CL_SIZE * 8);
  184.                 cacheStart = cacheStart >= sqrt ? cacheStart : sqrt;
  185.                 scan( 3, sqrt,
  186.                     [&]( size_t prime, auto )
  187.                     {
  188.                         size_t start = (cacheStart + prime - 1) / prime * prime;
  189.                         start = start % 2 ? start : start + prime;
  190.                         punch( start, cacheEnd, prime, true_type() );
  191.                     } );
  192.             }
  193.         } );
  194.     if( bool count = false; argc >= 3 && (!*argv[2] || (count = strcmp( argv[2], "*" ) == 0)) )
  195.     {
  196.         if( !count )
  197.             return EXIT_SUCCESS;
  198.         atomic<size_t> aNPrimes( 1 );
  199.         initiate( 3,
  200.             [&]( size_t rStart, size_t rEnd, auto )
  201.             {
  202.                 size_t nPrimes = 0;
  203.                 scan( rStart, rEnd, [&]( ... ) { ++nPrimes; } );
  204.                 aNPrimes.fetch_add( nPrimes, memory_order::relaxed );
  205.             } );
  206.         cout << aNPrimes.load( memory_order::relaxed ) << " primes" << endl;
  207.         return EXIT_SUCCESS;
  208.     }
  209.     ofstream ofs;
  210.     ofs.exceptions( ofstream::failbit | ofstream::badbit );
  211.     ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary | ofstream::trunc );
  212.     constexpr size_t
  213.         BUF_SIZE = 0x100000,
  214.         AHEAD = 32;
  215.     union ndi_char { char c; ndi_char() {} };
  216.     vector<ndi_char> rawPrintBuf( BUF_SIZE + AHEAD - 1 );
  217.     span printBuf( &rawPrintBuf.front().c, &rawPrintBuf.back().c + 1 );
  218.     auto
  219.         bufBegin = printBuf.begin(),
  220.         bufEnd = bufBegin,
  221.         bufFlush = bufBegin + BUF_SIZE;
  222.     auto print = [&]() { ofs << string_view( bufBegin, bufEnd ); };
  223.     auto printPrime = [&]( size_t prime, auto )
  224.     {
  225.         auto [ptr, ec] = to_chars( &*bufEnd, &bufEnd[AHEAD - 1], prime );
  226.         if( (bool)ec ) [[unlikely]]
  227.             throw system_error( (int)ec, generic_category(), "converson failed" );
  228.         bufEnd = ptr - &*bufBegin + bufBegin; // pointer to iterator - NOP
  229.         *bufEnd++ = '\n';
  230.         if( bufEnd >= bufFlush ) [[unlikely]]
  231.             print(), bufEnd = bufBegin;
  232.     };
  233.     printPrime( 2, false_type() );
  234.     scan( 3, end, printPrime );
  235.     print();
  236. }
  238. array<unsigned, 4> cpuId( unsigned code )
  239. {
  240.     array<unsigned, 4> regs;
  241. #if defined(_MSC_VER)
  242.     __cpuid( (int *)&regs[0], code );
  243. #elif defined(__GNUC__) || defined(__clang__)
  244.     __cpuid(code, regs[0], regs[1], regs[2], regs[3]);
  245. #endif
  246.     return regs;
  247. }
  249. bool smt()
  250. {
  251.     if( cpuId( 0 )[0] < 1 )
  252.         return false;
  253.     return cpuId( 1 )[3] >> 28 & 1;
  254. }
  256. size_t getL2Size()
  257. {
  258.     if( cpuId( 0x80000000u )[0] < 0x80000006u )
  259.         return 512ull * 1024;
  260.     return (cpuId( 0x80000006u )[2] >> 16) * 1024;
  261. }

字符串
几乎所有的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的文件大约需要额外的两秒钟。

展开查看全部
赞(0）分享  回复(0）举报  2024-01-09
deyfvvtc

deyfvvtc3#

这可能会加快一点：

  1. #include <algorithm>
  2. #include <iostream>
  3. #include <iterator>
  4. #include <vector>
  6. int main()
  7. {
  8.     std::vector<unsigned long long> numbers;
  9.     unsigned long long maximum = 4294967296;
  10.     for (unsigned long long i = 2; i <= maximum; ++i)
  11.     {
  12.         if (numbers.empty())
  13.         {
  14.             numbers.push_back(i);
  15.             continue;
  16.         }
  18.         if (std::none_of(numbers.begin(), numbers.end(), [&](unsigned long long p)
  19.         {
  20.             return i % p == 0;
  21.         }))
  22.         {
  23.             numbers.push_back(i);
  24.         }
  26.     }
  28.     std::cout << "Primes:  " << std::endl;
  29.     std::copy(numbers.begin(), numbers.end(), std::ostream_iterator<int>(std::cout, " "));
  31.     return 0;
  32. }

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

展开查看全部
赞(0）分享  回复(0）举报  2024-01-09
gtlvzcf8

gtlvzcf84#

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

赞(0）分享  回复(0）举报  2024-01-09
fv2wmkja

fv2wmkja5#

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

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

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

赞(0）分享  回复(0）举报  2024-01-09
w8rqjzmb

w8rqjzmb6#

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

  1. #include <stdlib.h>
  2. #include <stdio.h>
  4. #define ARRSIZE 4000000000
  5. #define MAXCALC ARRSIZE/2
  7. int main() {
  8.     long int z;
  9.     long int *arr = (long int*) malloc((ARRSIZE) * sizeof(long int));
  10.     for (long int x=3;x <= MAXCALC; x++) {
  11.         if (x % 10 == 3 || x % 10 == 7 || x % 10 == 9) {
  12.             for (long int y=3; y < MAXCALC; y++){
  13.                     z = x * y;
  14.                     if (z < ARRSIZE)
  15.                         arr[z] = 1;
  16.                     else
  17.                         break;
  18.             }
  19.         }
  20.     }
  21.     printf("2 3 5 ");
  22.     for (long int x=7; x < ARRSIZE; x++) {
  23.         if (x % 2 != 0 && x % 10 != 5)
  24.             if (arr[x] != 1)
  25.                 printf("%ld ", x);
  26.     }
  27.     printf("\n");
  29.     free(arr);
  30.     return EXIT_SUCCESS;
  31. }

字符串

展开查看全部
赞(0）分享  回复(0）举报  2024-01-09
uplii1fm

uplii1fm7#

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

  1. import numpy as np
  2. from math import isqrt
  4. def all_primes_less_than(n):
  5.     is_prime = np.full(n,True)
  6.     is_prime[:2] = False
  7.     is_prime[4::2] = False
  8.     for p in range(3,isqrt(n),2):
  9.         if is_prime[p]: is_prime[p*p::p] = False
  10.     return is_prime.nonzero()[0]

字符串

赞(0）分享  回复(0）举报  2024-01-09
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com