Add fast prime sieve (min_25)

content/maths/number-theory/fast-sieve.cpp

@@ -0,0 +1,94 @@
+/*
+ *Author:* min\_25
+ *Description:* Fast prime sieve up to $N$.
+ Returns a sorted vector of all primes $<= N$.
+ Uses a segmented sieve with wheel factorization (mod 30).
+ *Time:* ~0.5s for $N = 10^9$.
+ *Status:* Tested on Yosupo (enumerate\_primes)
+*/
+vector<int> sieve(int N, int Q = 17, int L = 1 << 15) {
+  static const int rs[] = {1, 7, 11, 13, 17, 19, 23, 29};
+  struct P {
+    P(int p) : p(p) {}
+    int p, pos[8];
+  };
+
+  auto prime_count = [](int N) -> int {
+    return N > 60184 ? N / (log(N) - 1.1) : max(1., N / (log(N) - 1.11)) + 1;
+  };
+
+  int v = sqrt(N), vv = sqrt(v);
+  vector<bool> isp(v + 1, true);
+  for (int i = 2; i <= vv; i++)
+    if (isp[i])
+      for (int j = i * i; j <= v; j += i)
+        isp[j] = false;
+
+  int rsize = prime_count(N + 30);
+  vector<int> primes = {2, 3, 5};
+  int psize = 3;
+  primes.resize(rsize);
+
+  vector<P> sprimes;
+  size_t pbeg = 0;
+  int prod = 1;
+
+  for (int p = 7; p <= v; p++) {
+    if (!isp[p])
+      continue;
+    if (p <= Q)
+      prod *= p, pbeg++, primes[psize++] = p;
+    auto pp = P(p);
+    for (int t = 0; t < 8; t++) {
+      int j = (p <= Q) ? p : p * p;
+      while (j % 30 != rs[t])
+        j += p << 1;
+      pp.pos[t] = j / 30;
+    }
+    sprimes.push_back(pp);
+  }
+
+  vector<unsigned char> pre(prod, 0xFF);
+  for (size_t pi = 0; pi < pbeg; pi++) {
+    auto pp = sprimes[pi];
+    int p = pp.p;
+    for (int t = 0; t < 8; t++) {
+      unsigned char m = ~(1 << t);
+      for (int i = pp.pos[t]; i < prod; i += p)
+        pre[i] &= m;
+    }
+  }
+
+  int block_size = (L + prod - 1) / prod * prod;
+  vector<unsigned char> block(block_size);
+  unsigned char *pblock = block.data();
+  int M = (N + 29) / 30;
+
+  for (int beg = 0; beg < M; beg += block_size, pblock -= block_size) {
+    int end = min(M, beg + block_size);
+    for (int i = beg; i < end; i += prod)
+      copy(pre.begin(), pre.end(), pblock + i);
+    if (beg == 0)
+      pblock[0] &= 0xFE;
+    for (size_t pi = pbeg; pi < sprimes.size(); pi++) {
+      auto &pp = sprimes[pi];
+      int p = pp.p;
+      for (int t = 0; t < 8; t++) {
+        int i = pp.pos[t];
+        unsigned char m = ~(1 << t);
+        for (; i < end; i += p)
+          pblock[i] &= m;
+        pp.pos[t] = i;
+      }
+    }
+    for (int i = beg; i < end; i++)
+      for (int m = pblock[i]; m > 0; m &= m - 1)
+        primes[psize++] = i * 30 + rs[__builtin_ctz(m)];
+  }
+
+  assert(psize <= rsize);
+  while (psize > 0 && primes[psize - 1] > N)
+    psize--;
+  primes.resize(psize);
+  return primes;
+}

tracker-wf.yaml

@@ -99,6 +99,7 @@ content:
     - counting-divisors.cpp
     - discrete-log.cpp
     - erathostenes-sieve.cpp
+    - fast-sieve.cpp
     - euler-phi.cpp
     - euler-phi-sieve.cpp
     - mobius-sieve.cpp

tracker.yaml

@@ -76,6 +76,7 @@ content:
     - divisors.cpp
     - discrete-log.cpp
     - erathostenes-sieve.cpp
+    - fast-sieve.cpp
     - euler-phi.cpp
     - euler-phi-sieve.cpp
     - mobius-sieve.cpp