fork download
copy 
  1. #include <iostream>
  2. #include <utility>
  3. #include <vector>
  4. #include <list>
  5. #include <cassert>
  6. #include <fstream>
  7. #include <cstring>
  8.  
  9. template <class T>
  10. class Sieve {
  11.   uint64_t buffsize;
  12.   T *buff;
  13.   static std::pair<uint64_t, uint64_t> index(uint64_t n) {
  14.     std::pair<uint64_t, uint64_t> retv;
  15.     retv.first = (n - 1) / (8 * sizeof(T)); /* zero origin */
  16.     uint64_t bit = (n - 1) % (8 * sizeof(T));    /* zero origin */
  17.     T mask = 1;
  18.     for (uint64_t i = 0; i < bit; i++) mask = (mask << 1);
  19.     retv.second = mask;
  20.     return retv;
  21.   }
  22.  
  23. public:
  24.   Sieve(uint64_t n) {
  25.     std::pair<uint64_t, uint64_t> retindexf = index(n);
  26.     buffsize = retindexf.first + 1; /* zero origin */
  27.     buff = new T [buffsize];
  28.     for (uint64_t i = 0; i < buffsize; i++) buff[i] = 0;
  29.   }
  30.   ~Sieve() {
  31.     delete buff;
  32.   }
  33.   void setSieve(uint64_t n) {
  34.     std::pair<uint64_t, uint64_t> maskpos = index(n);
  35.     buff[maskpos.first] = buff[maskpos.first] | maskpos.second;
  36.   }
  37.   bool isSieveTrue(uint64_t n) {
  38.     std::pair<uint64_t, uint64_t> maskpos = index(n);    
  39.     return buff[maskpos.first] & maskpos.second;
  40.   }
  41. };
  42.  
  43. template <class T>
  44. class SieveDerived : private Sieve<T> {
  45.   static int index(uint64_t n) {
  46.     int r1 = n / 30;
  47.     int r2 = n % 30;
  48.     switch (r2) {
  49.     case 0:
  50.       /* 1 */
  51.     case 2:
  52.     case 3:
  53.     case 4:
  54.     case 5:
  55.     case 6:
  56.       /* 7 */
  57.     case 8:
  58.     case 9:
  59.     case 10:
  60.       /* 11 */
  61.     case 12:
  62.       /* 13 */
  63.     case 14:
  64.     case 15:
  65.     case 16:
  66.       /* 17 */
  67.     case 18:
  68.       /* 19 */
  69.     case 20:
  70.     case 21:
  71.     case 22:
  72.       /* 23 */
  73.     case 24:
  74.     case 25:
  75.     case 26:
  76.     case 27:
  77.     case 28:
  78.       /* 29 */
  79.       return 0;
  80.  
  81.     case 1:
  82.       return 8 * r1 + 1;
  83.     case 7:
  84.       return 8 * r1 + 2;
  85.     case 11:
  86.       return 8 * r1 + 3;
  87.     case 13:
  88.       return 8 * r1 + 4;
  89.     case 17:
  90.       return 8 * r1 + 5;
  91.     case 19:
  92.       return 8 * r1 + 6;
  93.     case 23:
  94.       return 8 * r1 + 7;
  95.     case 29:
  96.       return 8 * r1 + 8;
  97.  
  98.     default:
  99.       return 0;
  100.     }
  101.   }
  102.  
  103. public:
  104. //  SieveDerived(int n) : Sieve<T>(({int r; while((r = index(n)) == 0) n--; r;})) { } /* HERE!! */
  105.   SieveDerived(uint64_t n) : Sieve<T>([](uint64_t n)->uint64_t{uint64_t r; while ((r = index(n)) == 0) n--; return r;}(n)) { } /* HERE!! */
  106.   ~SieveDerived() { }
  107.   void setSieve(uint64_t n) {
  108.     uint64_t r = index(n);
  109.     if (r == 0) return;
  110.     return ((Sieve<T> *)this)->setSieve(r);
  111.   }
  112.   bool isSieveTrue(int n) {
  113.     uint64_t r = index(n);
  114.     if (r == 0) return true;
  115.     return ((Sieve<T> *)this)->isSieveTrue(index(n));
  116.   }
  117. };
  118.  
  119. template <class T>
  120. class Array : public std::vector<T> {
  121. public:
  122.   Array() : std::vector<T>() { }
  123.   ~Array() { }
  124.   T &operator[](std::size_t idx) {
  125.     if (idx >= this->size()) {
  126.       while (this->size() <= idx) this->push_back(0);
  127.     }
  128.     return (*(std::vector<T> *)this)[idx];
  129.   }
  130. };
  131. /*------------------------------------------------------------------*/
  132. //                                   123456789AB
  133. constexpr uint64_t MAXNUM_Default = 100000000000;
  134. constexpr unsigned int windowSize_Default =   30 * 5000000;
  135. char const filename[] = "prime.bin";
  136.  
  137. Array<uint64_t> array;
  138. uint64_t maxSpanLast = 0;
  139.  
  140. void spancheck(uint64_t prime, uint64_t n) {
  141.   uint64_t maxSpanNowFreq = 0; uint64_t maxSpanNow = 0;
  142.   for (uint64_t k = 1; k < array.size(); k++)
  143.     if (array[k] > maxSpanNowFreq) {
  144.       maxSpanNow = k;
  145.       maxSpanNowFreq = array[k];
  146.     }
  147.   if (maxSpanNow != maxSpanLast) {
  148.     std::cout << "span = " << maxSpanNow * 2 << ", p = " << prime << ", n = " << n << std::endl;
  149.     maxSpanLast = maxSpanNow;
  150.   }
  151. }
  152.  
  153. void registspan(uint64_t &n, uint64_t prime, uint64_t &primeNow, uint64_t &primeLast) {
  154.   n++;
  155. //  std::cout << "size() = " << primeTable.size() << " " << std::flush;
  156.   primeNow = prime;
  157.   array[(primeNow - primeLast) / 2]++;
  158.   primeLast = primeNow;
  159. }
  160.  
  161. static bool flag_continue = 0;
  162. static uint64_t MAXNUM = 0;
  163. static unsigned int windowSize = 0;
  164.  
  165. int main(int argc, char *argv[]) {
  166.   if (argc > 4) {
  167.   label_usage:
  168.     std::cerr << argv[0] << " : search prime-to-prime span." << std::endl;
  169.     std::cerr << "usage: " << argv[0] << " <search-max-range> <sieve-range>" << std::endl;
  170.     std::cerr << "option: /continue : continuing search using generated table : 'prime.bin' already :" << filename << "." << std::endl;
  171.     return 0;
  172.   }
  173.  
  174.   int nowstep = 0;
  175.   for (int i = 1; i < argc; i++) {
  176.     if (strstr(argv[i], "/continue") == argv[i]) {
  177.       flag_continue = true;
  178.       continue;
  179.     }
  180.     switch (nowstep) {
  181.     case 0:
  182.       MAXNUM = atoi(argv[i]);
  183.       nowstep++;
  184.       continue;
  185.     case 1:
  186.       windowSize = atoi(argv[i]);
  187.       nowstep++;
  188.       continue;
  189.     default:
  190.       goto label_usage;
  191.     }
  192.   }
  193.   if (MAXNUM == 0) MAXNUM = MAXNUM_Default;
  194.   if (windowSize == 0) windowSize = windowSize_Default;
  195.   if (windowSize % 30 != 0) {
  196.     std::cerr << "bad windowSize, aborted." << std::endl;
  197.     return 0;
  198.   }
  199.  
  200.   uint64_t n = 0;
  201.   uint64_t primeLast, primeNow;
  202.   uint64_t offset = 0;
  203.  
  204.   std::fstream f;
  205.   std::fstream::pos_type hasReadPosition;
  206.   std::fstream::pos_type hasWrittenPosition;
  207.  
  208.  
  209.   if (!flag_continue) {
  210.      f.open(filename, std::ios::in | std::ios::out | std::ios::binary | std::ios::trunc );
  211.   } else {
  212.      f.open(filename, std::ios::in | std::ios::out | std::ios::binary);
  213.   }
  214.   if (!f) { std::cerr << "cannnot open the file, aborted." << std::endl; exit(0); }
  215.  
  216.   if (!flag_continue) {
  217.     uint64_t w;
  218.  
  219.     f.seekp(0, std::ios_base::beg); /* seekput1-1 */
  220.  
  221.     SieveDerived<uint64_t> sieve(windowSize);
  222.     sieve.setSieve(1);
  223. //      sieve.setSieve(2); std::cout << 2 << ", ";
  224.     w = 2; f.write((char *)&w, sizeof(w));
  225.     n++;
  226. //      sieve.setSieve(3); std::cout << 3 << ", ";
  227.     w = 3; f.write((char *)&w, sizeof(w));
  228.     n++;
  229.     primeLast = 3;
  230.     sieve.setSieve(4);
  231. //      sieve.setSieve(5); std::cout << 5 << ", ";
  232.     w = 5; f.write((char *)&w, sizeof(w));
  233.     hasWrittenPosition = f.tellp();
  234.  
  235.     n++;
  236.     primeNow = 5;
  237.     array[(primeNow - primeLast) / 2]++;
  238.     primeLast = primeNow;
  239.     sieve.setSieve(6);
  240.     for (uint64_t i = 7; i <= windowSize; i++) {
  241.       if (i > MAXNUM)
  242.         break;
  243.       if (!sieve.isSieveTrue(i)) {
  244.         std::cout << i << ", " << std::flush;
  245.  
  246.         f.seekp(hasWrittenPosition, std::ios_base::beg);
  247.         f.write(reinterpret_cast<char *>(&i), sizeof(i));
  248.         hasWrittenPosition = f.tellp();
  249.  
  250.         registspan(n, i, primeNow, primeLast);
  251.         for (uint64_t j = i * 2; j <= windowSize; j += i) sieve.setSieve(j);
  252.         spancheck(i, n);
  253.       }
  254.     }
  255.     offset = windowSize;
  256.  
  257.   } else { /* continue */
  258.     uint64_t data, last;
  259.     std::cout << "continue." << std::endl;
  260.  
  261.     /* step.1 */
  262.     f.seekg(0, std::ios_base::beg);
  263.     data = 0;
  264.     for (;;) {
  265.       hasReadPosition = f.tellg();
  266.       f.read((char *)&data, sizeof(data));
  267.  
  268.       if (!f.good())
  269.         break;
  270.       if (data == 2) { n++; continue; }
  271.       if (data == 3) { n++; primeLast = data; continue; }
  272.       if (data <= primeLast) {
  273.         std::cerr << "data = " << data << ", primeLast = " << primeLast << std::endl;
  274.         std::cerr << "data corupted, aborted.(1)" << std::endl;
  275.         std::cerr << "prime.bin is chunked. pos = " << hasReadPosition << std::endl;
  276.         data = primeLast;
  277.         break;
  278.       }
  279.       registspan(n, data, primeNow, primeLast);
  280.       if (n % 100000000 == 0) std::cout << "n = " << n << std::endl;
  281.     }
  282.     if (data == 0) { std::cerr << "data corupted, aborted.(2)" << std::endl; exit(0); }
  283.     last = data;
  284.     std::cout << "last = " << last << std::endl;
  285.     last++;
  286.     f.clear(); /* needed */
  287.     f.seekp(hasReadPosition, std::ios_base::beg);  /* seekput2-1 */
  288.     hasWrittenPosition = f.tellp();
  289.  
  290.     /* step.2 */
  291.     int flag;
  292.     uint64_t primeStep2;
  293.     uint64_t near_offset = ((last + (windowSize - 1))/ windowSize) * windowSize;
  294.  
  295.     for (primeStep2 = last; primeStep2 < near_offset; primeStep2++) {
  296.  
  297.       f.seekg(0, std::ios_base::beg);
  298.       hasReadPosition = f.tellg();
  299.  
  300.       for (;;) {
  301.         uint64_t p;
  302.  
  303.         f.seekg(hasReadPosition, std::ios_base::beg);
  304.         f.read((char *)&p, sizeof(p));
  305.         hasReadPosition = f.tellg();
  306.  
  307.         if (primeStep2 % p == 0) { flag = 1; break; }
  308.         if (primeStep2 / p < p) { flag = 2; break; }
  309.       }
  310.       if (flag == 1) continue;
  311.       assert(flag == 2);
  312.  
  313.       f.seekp(hasWrittenPosition, std::ios_base::beg);
  314.       f.write(reinterpret_cast<char *>(&primeStep2), sizeof(primeStep2));
  315.       hasWrittenPosition = f.tellp();
  316.  
  317.       std::cout << "prime = " << primeStep2 << std::endl;
  318.       registspan(n, primeStep2, primeNow, primeLast);
  319.     }
  320.  
  321.     offset = primeStep2;
  322.     std::cout << "offset = " << offset << std::endl;
  323.     std::cout << "n = " << n << std::endl;
  324.   }
  325.  
  326.   for (; offset < MAXNUM; offset += windowSize) {
  327.     std::cout << "offset = " << offset << std::endl;
  328.     SieveDerived<uint64_t> sieve(windowSize);
  329.  
  330.     f.seekg(0, std::ios_base::beg);
  331.     hasReadPosition = f.tellg();
  332.  
  333.     for (;;) {
  334.       /* loop out 1 */
  335.       if ((uint64_t)f.tellg() >= n * sizeof(uint64_t))
  336.         break;
  337.       uint64_t p;
  338.  
  339.       f.seekg(hasReadPosition, std::ios_base::beg);      
  340.       f.read((char *)&p, sizeof(p));
  341.       hasReadPosition = f.tellg();
  342.  
  343.       /* loop out 2 */
  344.       if (p > (offset + windowSize) / 2)
  345.         break;
  346.       /* loop body */
  347.       for (uint64_t i = offset - (offset % p) + p; i < offset + windowSize; i += p)
  348.         sieve.setSieve(i - offset);
  349.     }
  350.  
  351.     for (uint64_t index = 0; index < windowSize; index++) {
  352.       if (offset + index > MAXNUM)
  353.         break;
  354.       if (!sieve.isSieveTrue(index)) {
  355.         std::cout << offset + index << ", " << std::flush;
  356.         uint64_t prime = offset + index;
  357.  
  358.         f.seekp(hasWrittenPosition, std::ios_base::beg);
  359.         f.write(reinterpret_cast<char *>(&prime), sizeof(prime));
  360.         hasWrittenPosition = f.tellp();
  361.  
  362.         registspan(n, prime, primeNow, primeLast);
  363.         spancheck(prime, n);
  364.  
  365.       }
  366.     }
  367.   } /* for (offset = windowSize; ...) */
  368.   f.close();
  369.   std::cout << std::endl << "n = " << n << std::endl;
  370.   for (uint64_t k = 0; k < array.size(); k++) {
  371.     if (array[k] > 0)
  372.       std::cout << k * 2 << ":" << array[k] << std::endl;
  373.   }
  374.   return 0;
  375. }
  376. /* end */
  377.  
Success #stdin #stdout #stderr 0s 4320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
cannnot open the file, aborted.
https://ideone.com/6Ww9nq
language:
C++14 (gcc 8.3)
created:
4 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!