fork(1) download
copy 
  1. #include <iostream>
  2. #include <fstream>
  3. #include <algorithm>
  4. #include <list>
  5. #include <vector>
  6. #include <cmath>
  7. #include <cassert>
  8.  
  9. using namespace std;
  10.  
  11. class PrimeGenerator {
  12.   public:
  13.     PrimeGenerator() {
  14.         normalSolution();
  15.     }
  16.  
  17.     int nextPrime() {
  18.         int p = m_primes.front();
  19.         m_primes.pop_front();
  20.         // https://e...content-available-to-author-only...a.org/wiki/Prime-counting_function
  21.         // pi(m_size) > m_size / ln(m_size) for m_size >= 17 (which starts at 100 so its true)
  22.         // pi(m_size / 4) <= 1.25506 * (m_size / 4) / ln(m_size / 4) <= 1.25506 * (m_size / 4) / ln(m_size ^ 0.5) = 1.25506 * (m_size / 2) / ln(m_size)
  23.         // so the number of primes bigger than m_size / 4 but lower than or equal to m_size is
  24.         // > m_size / ln(m_size) - 1.25506 * (m_size / 2) / ln(m_size) = m_size / ln(m_size) * 0.37247
  25.         // so doSomeSteps is called at least m_size / ln(m_size) * 0.37247
  26.         // and because doSomeSteps does O(ln(m_size)) operations and p > m_size / 4 it does O(ln p) operations
  27.         if (p > m_size / 4)
  28.             doSomeSteps();
  29.         return p;
  30.     }
  31.  
  32.   private:
  33.     // iterator, meant to hold a state in normalSolution
  34.     struct Iterator {
  35.         int step; // 0 -> filling with 0, 1 -> doing the sieve, 2 -> removing primes
  36.         int i;
  37.         list<int>::iterator p;
  38.     } m_iterator;
  39.  
  40.     void normalSolution() {
  41.         m_size = START_SIZE;
  42.         m_calls = ceil(log(m_size) * 43);
  43.         m_lowest_prime = vector<int>(m_size + 1, 0);
  44.         // Sieve of erathostenes in linear time,
  45.         // this algorithm makes N + pi(N) + N operations, upper bounded by 3N
  46.         // the if is true in pi(N) cases <= N
  47.         // we update at most once each element => N operations
  48.         // the for condition will have to break one for each element => N operations
  49.         for (int i = 2; i <= m_size; ++i) {
  50.             if (m_lowest_prime[i] == 0) {
  51.                 m_lowest_prime[i] = i;
  52.                 m_primes.push_back(i);
  53.             }
  54.             for (auto &prime : m_primes) {
  55.                 if (prime > m_lowest_prime[i])
  56.                     break;
  57.                 if (i * prime > m_size)
  58.                     break;
  59.                 m_lowest_prime[prime * i] = prime;
  60.             }
  61.         }
  62.         m_lowest_prime_second = vector<int>(m_size * 4 + 1);
  63.         m_iterator = {0, 0, m_primes_second.begin()};
  64.     }
  65.  
  66.     // to do the sieve on 4 * m_size we need 16 * m_size operations
  67.     // 4 * m_size to fill the array with 0
  68.     // 8 * m_size to simulate the sieve in linear time
  69.     // pi(m_size * 4) to remove the already printed primes, <= 4 * m_size
  70.     // and we have at least m_size / ln(m_size) * 0.37247 calls, so to distribute them all
  71.     // we need to do at least 16ln(m_size) / 0.37247 ~ 43 ln(m_size) which is O(ln(m_size))
  72.     void doSomeSteps() {
  73.         int operations = m_calls;
  74.         for (int i = 0; i < operations; ++i)
  75.             doStep();
  76.         if (done())
  77.             swap();
  78.     }
  79.  
  80.     // complexity of this is O(1)
  81.     void doStep() {
  82.         if (m_iterator.step == 0) { // filling with 0, this will happen exactly 4 * m_size times
  83.             m_lowest_prime_second[m_iterator.i++] = 0;
  84.             if (m_iterator.i > m_size * 4)
  85.                 m_iterator = {1, 2, m_primes_second.begin()};
  86.             return;
  87.         }
  88.         if (m_iterator.step == 1) { // doing the actual sieve
  89.             // when its prime, this will be at most pi(4 * m_size) <= 4 * m_size
  90.             if (m_lowest_prime_second[m_iterator.i] == 0) {
  91.                 m_lowest_prime_second[m_iterator.i] = m_iterator.i;
  92.                 m_primes_second.push_back(m_iterator.i);
  93.                 m_iterator.p = m_primes_second.begin();
  94.                 return;
  95.             }
  96.             // break condition, this will happen exactly 4 * m_size
  97.             if (m_iterator.p == m_primes_second.end()
  98.                 || *m_iterator.p > m_lowest_prime_second[m_iterator.i]
  99.                 || *m_iterator.p * m_iterator.i > m_size * 4) {
  100.                 m_iterator = {1, m_iterator.i + 1, m_primes_second.begin()};
  101.                 if (m_iterator.i > 4 * m_size)
  102.                     m_iterator = {2, 0, m_primes_second.end()};
  103.             } else {
  104.                 m_lowest_prime_second[*m_iterator.p * m_iterator.i] = *m_iterator.p;
  105.                 m_iterator.p++;
  106.             }
  107.             return;
  108.         }
  109.         if (m_iterator.step == 2) { // removing primes
  110.             if (m_primes_second.front() <= m_size)
  111.                 m_primes_second.pop_front();
  112.             else
  113.                 m_iterator.step = 3;
  114.             return;
  115.         }
  116.     }
  117.  
  118.     // when we have finished creating the second sieve (which is 4 times bigger) and exhausted the primes from the first
  119.     bool done() {
  120.         return m_iterator.step == 3 && m_primes.begin() == m_primes.end();
  121.     }
  122.  
  123.     void swap() {
  124.         ::swap(m_lowest_prime, m_lowest_prime_second);
  125.         ::swap(m_primes, m_primes_second);
  126.         m_size *= 4;
  127.         m_calls = ceil(43 * log(m_size));
  128.         m_lowest_prime_second = vector<int>(4 * m_size + 1);
  129.         m_iterator = {0, 0, m_primes.begin()};
  130.     }
  131.     list<int> m_primes, m_primes_second;
  132.     vector<int> m_lowest_prime, m_lowest_prime_second;
  133.     int m_size, m_calls;
  134.  
  135.     static const int START_SIZE = 100;
  136. };
  137.  
  138. vector<int> primes(int N) {
  139.     vector<bool> sieve(N + 1, true);
  140.     vector<int> primes;
  141.     for (int i = 2; i * i <= N; ++i)
  142.         if (sieve[i] == true) {
  143.             for (int j = i * i; j <= N; j += i)
  144.                 sieve[j] = false;
  145.         }
  146.     for (int i = 2; i <= N; ++i)
  147.         if (sieve[i] == true)
  148.             primes.push_back(i);
  149.     return primes;
  150. }
  151.  
  152. vector<int> generator(int P) {
  153.     vector<int> primes;
  154.     PrimeGenerator G;
  155.     while (P--)
  156.         primes.push_back(G.nextPrime());
  157.     return primes;
  158. }
  159.  
  160. int main() {
  161.     ios::sync_with_stdio(false);
  162.  
  163.     int N = 10 * 1000 * 1000;
  164.     auto normalSieve = primes(N);
  165.     auto smartGenerator = generator(normalSieve.size());
  166.     assert(normalSieve == smartGenerator);
  167. }
  168.  
Success #stdin #stdout 4.49s 32240KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/NnOU63
language:
C++14 (gcc 8.3)
created:
8 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!