fork(6) download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. const int N = 1e8;
  4. bool comp[N];
  5. vector<int> primes;
  6.  
  7. void normal_sieve() {
  8.     memset(comp, 0, sizeof comp);
  9.  
  10.     primes.clear();
  11.     bool done = false;
  12.     for (int i = 2; i < N; i++) {
  13.         if (!comp[i]) {
  14.             primes.push_back(i);
  15.  
  16.             if (i * i > N) {
  17.                 done = true;
  18.             }
  19.  
  20.             if (!done)
  21.                 for (int j = i * i; j < N; j += i)
  22.                     comp[j] = 1;
  23.         }
  24.     }
  25. }
  26.  
  27. void normal_optimized_sieve() {
  28.     memset(comp, 0, (sizeof comp) / 2);
  29.     primes.clear();
  30.  
  31.     primes.push_back(2);
  32.     int n = N / 2;
  33.     int sr = round(sqrt(N));
  34.     for (int i = 1; i < n; i++) {
  35.         if (!comp[i]) {
  36.             int p = 2 * i + 1;
  37.             primes.push_back(p);
  38.  
  39.             if (p > sr)
  40.                 continue;
  41.  
  42.             for (int j = (p + 1) * i; j < n; j += p)
  43.                 comp[j] = 1;
  44.         }
  45.     }
  46. }
  47.  
  48. void linear_sieve() {
  49.     memset(comp, 0, sizeof comp);
  50.     primes.clear();
  51.     for (int i = 2; i < N; i++) {
  52.         if (!comp[i])
  53.             primes.push_back(i);
  54.         for (int j = 0, si = primes.size(); j < si && i * primes[j] < N; j++) {
  55.             comp[primes[j] * i] = 1;
  56.             if (i % primes[j] == 0)
  57.                 break;
  58.         }
  59.     }
  60. }
  61.  
  62. void linear_sieve_optimized() {
  63.     memset(comp, 0, (sizeof comp) / 2);
  64.     primes.clear();
  65.     primes.push_back(2);
  66.     for (int i = 1, p = 3; p < N; i++, p+=2) {
  67.         if (!comp[i])
  68.             primes.push_back(p);
  69.         for (int j = 1, si = primes.size(); j < si && p * primes[j] < N; j++) {
  70.             comp[(primes[j] * p) >> 1] = 1;
  71.             if (p % primes[j] == 0)
  72.                 break;
  73.         }
  74.     }
  75. }
  76.  
  77. void segmented_sieve() {
  78.     int n = N;
  79.     const int S = 50000;
  80.  
  81.     primes.clear();
  82.     int nsqrt = round(sqrt(n));
  83.     vector<char> is_prime(nsqrt + 1, true);
  84.  
  85.     ::primes.push_back(2);
  86.     vector<char> block(S);
  87.  
  88.     vector<int> primes, start_indices;
  89.     for (int i = 3; i <= nsqrt; i += 2) {
  90.         if (is_prime[i]) {
  91.             primes.push_back(i);
  92.             start_indices.push_back(i * i);
  93.             for (int j = i * i; j <= nsqrt; j += 2 * i)
  94.                 is_prime[j] = false;
  95.         }
  96.     }
  97.  
  98.     for (int k = 0; k * S <= n; k++) {
  99.         fill(block.begin(), block.end(), true);
  100.         int start = k * S;
  101.         for (auto i = 0u; i < primes.size(); i++) {
  102.             int p = primes[i];
  103.             int idx = start_indices[i];
  104.             for (; idx < S; idx += 2 * p)
  105.                 block[idx] = false;
  106.             start_indices[i] = idx - S;
  107.         }
  108.         if (k == 0)
  109.             block[1] = false;
  110.         for (int i = 1; i < S && start + i <= n; i += 2) {
  111.             if (block[i])
  112.                 ::primes.push_back(start + i);
  113.         }
  114.     };
  115. }
  116.  
  117. void segmented_sieve_optimized() {
  118.     int n = N;
  119.     const int S = 30000;
  120.  
  121.     primes.clear();
  122.     int nsqrt = round(sqrt(n));
  123.  
  124.     vector<char> is_prime(nsqrt + 1, true);
  125.     vector<int> primes, start_indices;
  126.     for (int i = 3; i <= nsqrt; i += 2) {
  127.         if (is_prime[i]) {
  128.             primes.push_back(i);
  129.             start_indices.push_back((i * i - 1) / 2);
  130.             for (int j = i * i; j <= nsqrt; j += 2 * i)
  131.                 is_prime[j] = false;
  132.         }
  133.     }
  134.  
  135.     ::primes.push_back(2);
  136.     vector<char> block(S);
  137.     int high = (n - 1) / 2;
  138.     for (int low = 0; low <= high; low += S) {
  139.         fill(block.begin(), block.end(), true);
  140.         for (auto i = 0u; i < primes.size(); i++) {
  141.             int p = primes[i];
  142.             int idx = start_indices[i];
  143.             for (; idx < S; idx += p)
  144.                 block[idx] = false;
  145.             start_indices[i] = idx - S;
  146.         }
  147.         if (low == 0)
  148.             block[0] = false;
  149.         for (int i = 0; i < S && low + i <= high; i++) {
  150.             if (block[i])
  151.                 ::primes.push_back((low + i) * 2 + 1);
  152.         }
  153.     };
  154. }
  155.  
  156. void aitken() {
  157.     int n = N;
  158.     std::vector<bool> m(n + 1, false);
  159.  
  160.     int root = std::sqrt(n) + 1;
  161.     for (int i = 1; i <= root; i++) {
  162.         for (int j = 1; j <= root; j++) {
  163.             int a = 4 * i * i + j * j;
  164.             int b = 3 * i * i + j * j;
  165.             int c = 3 * i * i - j * j;
  166.  
  167.             if (a <= n && (a % 12 == 1 || a % 12 == 5))
  168.                 m[a].flip();
  169.             if (b <= n && b % 12 == 7)
  170.                 m[b].flip();
  171.             if (i > j && c <= n && c % 12 == 11)
  172.                 m[c].flip();
  173.         }
  174.     }
  175.  
  176.     for (int i = 5; i * i <= n; i++) {
  177.         if (m[i]) {
  178.             for (int j = 1; j * i * i <= n; j++)
  179.                 m[j * i * i] = false;
  180.         }
  181.     }
  182.  
  183.     primes.clear();
  184.     primes.push_back(2);
  185.     primes.push_back(3);
  186.     for (int i = 5; i < n; i++) {
  187.         if (m[i])
  188.             primes.push_back(i);
  189.     }
  190. }
  191.  
  192. template <typename T>
  193. void time_compare(string name, T func, vector<int> comparison) {
  194.     primes.clear();
  195.     auto t0 = clock();
  196.     (*func)();
  197.     auto t1 = clock();
  198.     if (!comparison.empty())
  199.         assert(comparison == primes);
  200.     cout << setw(20) << name << ": " << (1.0 * (t1 - t0)) / CLOCKS_PER_SEC << " seconds"
  201.          << endl;
  202. }
  203.  
  204. int main() {
  205.     ::primes.reserve(6'000'00);
  206.  
  207.     time_compare("normal", normal_sieve, {});
  208.     auto comparison = primes;
  209.     time_compare("normal_optimized", normal_optimized_sieve, comparison);
  210.     time_compare("linear", linear_sieve, comparison);
  211.     time_compare("linear_optimized", linear_sieve_optimized, comparison);
  212.     time_compare("segmented", segmented_sieve, comparison);
  213.     time_compare("segmented_optimized", segmented_sieve_optimized, comparison);
  214.     time_compare("aitken", aitken, comparison);
  215. }
  216.  
Success #stdin #stdout 3.6s 147648KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
              normal: 1.09295 seconds
    normal_optimized: 0.468311 seconds
              linear: 0.546933 seconds
    linear_optimized: 0.321325 seconds
           segmented: 0.188974 seconds
 segmented_optimized: 0.133987 seconds
              aitken: 0.800168 seconds
https://ideone.com/I4piBg
language:
C++14 (gcc 8.3)
created:
6 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!