Sieve in O(N)

Author: hitonanode

number/README.md

@@ -2,8 +2,14 @@
 
 ## Contents
 
-- `eratosthenes.hpp` Sieve of eratosthenes
+- `count_primes.hpp` Prime-counting function, O(N^(2/3))
 - `cyclotomic_polynomials.hpp` Cyclotomic polynomial （円分多項式, $x^k - 1$の因数分解など）
+- `discrete_logarithm.hpp` Discrete logarithm by the baby-step giant-step algorithm
+- `enumerate_partitions.hpp` 自然数の分割の列挙
+- `euler_toptient_phi.hpp` Generate table of Euler's totient function
+- `factorize.hpp` 大きい整数の素因数分解・素数判定
+- `modint_runtime.hpp` 実行時 modint （`modint.hpp` と共通のインターフェース）
+- `sieve.hpp` 線形篩・素因数分解・メビウス関数
 
 ## Cyclotomic Polynomials
 

number/count_primes.hpp

@@ -5,7 +5,7 @@
 
 // CUT begin
 struct CountPrimes {
-    // Count Primes less than x (\pi(x)) for each x = N / i (i = 1, ..., N) in O(N^(2/3)) time
+    // Count Primes less than or equal to x (\pi(x)) for each x = N / i (i = 1, ..., N) in O(N^(2/3)) time
     // Learned this algorihtm from <https://old.yosupo.jp/submission/14650>
     // Reference: <https://min-25.hatenablog.com/entry/2018/11/11/172216>
     using Int = long long;
@@ -36,7 +36,8 @@ struct CountPrimes {
                 for (size_t j = p * 2; j < is_prime.size(); j += p) is_prime[j] = 0;
             }
         }
-        for (Int now = n; now; now = n / (n / now + 1)) vs.push_back(now); // [N, N / 2, ..., 1], Relevant integers (decreasing) length ~= 2sqrt(N)
+        for (Int now = n; now; now = n / (n / now + 1))
+            vs.push_back(now); // [N, N / 2, ..., 1], Relevant integers (decreasing) length ~= 2sqrt(N)
         s = vs.size();
 
         // pi[i] = (# of integers x s.t. x <= vs[i],  (x is prime or all factors of x >= p))
@@ -74,8 +75,9 @@ struct CountPrimes {
         };
         for (; primes[ip] <= n3; ip++, pre++) {
             const auto &p = primes[ip];
-            for (int i = 0; i < n3 and p * p <= vs[i]; i++) trans2(i, p); // upto n3, total trans2 times: O(N^(2/3) / logN)
-            _fenwick_rec_update(ip, primes[ip], true);                    // total update times: O(N^(2/3) / logN)
+            for (int i = 0; i < n3 and p * p <= vs[i]; i++)
+                trans2(i, p);                          // upto n3, total trans2 times: O(N^(2/3) / logN)
+            _fenwick_rec_update(ip, primes[ip], true); // total update times: O(N^(2/3) / logN)
         }
         for (int i = s - n3 - 1; i >= 0; i--) {
             int j = i + ((i + 1) & (-i - 1));

number/cyclotomic_polynomials.hpp

@@ -1,12 +1,12 @@
 #pragma once
-#include "number/eratosthenes.hpp"
+#include "sieve.hpp"
 #include <vector>
 
 // CUT begin
 // cyclotomic polynominals : 円分多項式
 // ret[i][j] = [x^j]\Phi_i(x) for i <= Nmax, which are known to be integers
 std::vector<std::vector<int>> cyclotomic_polynomials(int Nmax) {
-    std::vector<int> moebius = SieveOfEratosthenes(Nmax).GenerateMoebiusFunctionTable();
+    std::vector<int> moebius = Sieve(Nmax).GenerateMoebiusFunctionTable();
     std::vector<std::vector<int>> ret(Nmax + 1);
     auto multiply = [](const std::vector<int> &a, const std::vector<int> &b) { // a * b
         std::vector<int> ret(int(a.size() + b.size()) - 1);

number/enumerate_partitions.hpp

@@ -4,11 +4,12 @@
 
 // CUT begin
 // Enumerate Partitions of number <= N （自然数の分割の列挙）
-// - can be used for N less than ~50
+// N = 50: 221 MB, 274 ms (length = 1295970) (g++, C++17, AtCoder)
 // - Example:
 //   - 1 => [[1,],]
 //   - 2 => [[1,],[1,1,],[2,],]
 //   - 3 => [[1,],[1,1,],[1,1,1,],[2,],[2,1,],[3,],]
+// 分割数の一覧：<https://oeis.org/A000041/list>
 struct {
     std::vector<std::vector<int>> ret;
     std::vector<int> now;

number/eratosthenes.hpp

@@ -9,7 +9,7 @@
 // (*this)[i] = (divisor of i, greater than 1)
 // Example: [0, 1, 2, 3, 2, 5, 3, 7, 2, 3, 2, 11, ...]
 // Complexity: Space O(MAXN), Time (construction) O(MAXNloglogMAXN)
-struct SieveOfEratosthenes : std::vector<int> {
+struct [[deprecated("Use sieve.hpp (linear time complexity)")]] SieveOfEratosthenes : std::vector<int> {
     std::vector<int> primes;
     SieveOfEratosthenes(int MAXN) : std::vector<int>(MAXN + 1) {
         std::iota(begin(), end(), 0);

number/sieve.hpp

@@ -0,0 +1,76 @@
+#pragma once
+#include <cassert>
+#include <map>
+#include <vector>
+
+// CUT begin
+// Linear sieve algorithm for fast prime factorization
+// Complexity: O(N) time, O(N) space:
+// - MAXN = 10^7:  ~44 MB,  80~100 ms (Codeforces / AtCoder GCC, C++17)
+// - MAXN = 10^8: ~435 MB, 810~980 ms (Codeforces / AtCoder GCC, C++17)
+// Reference:
+// [1] D. Gries, J. Misra, "A Linear Sieve Algorithm for Finding Prime Numbers,"
+//     Communications of the ACM, 21(12), 999-1003, 1978.
+// - <https://cp-algorithms.com/algebra/prime-sieve-linear.html>
+// - <https://37zigen.com/linear-sieve/>
+struct Sieve {
+    std::vector<int> min_factor;
+    std::vector<int> primes;
+    Sieve(int MAXN) : min_factor(MAXN + 1) {
+        for (int d = 2; d <= MAXN; d++) {
+            if (!min_factor[d]) {
+                min_factor[d] = d;
+                primes.emplace_back(d);
+            }
+            for (const auto &p : primes) {
+                if (p > min_factor[d] or d * p > MAXN) break;
+                min_factor[d * p] = p;
+            }
+        }
+    }
+    // Prime factorization for 1 <= x <= MAXN^2
+    // Complexity: O(log x)           (x <= MAXN)
+    //             O(MAXN / log MAXN) (MAXN < x <= MAXN^2)
+    template <typename T> std::map<T, int> factorize(T x) {
+        std::map<T, int> ret;
+        assert(x > 0 and x <= ((long long)min_factor.size() - 1) * ((long long)min_factor.size() - 1));
+        for (const auto &p : primes) {
+            if (x < T(min_factor.size())) break;
+            while (!(x % p)) x /= p, ret[p]++;
+        }
+        if (x >= T(min_factor.size())) ret[x]++, x = 1;
+        while (x > 1) ret[min_factor[x]]++, x /= min_factor[x];
+        return ret;
+    }
+    // Enumerate divisors of 1 <= x <= MAXN^2
+    // Be careful of highly composite numbers <https://oeis.org/A002182/list> <https://gist.github.com/dario2994/fb4713f252ca86c1254d#file-list-txt>:
+    // (n, (# of div. of n)): 45360->100, 735134400(<1e9)->1344, 963761198400(<1e12)->6720
+    template <typename T> std::vector<T> divisors(T x) {
+        std::vector<T> ret{1};
+        for (const auto p : factorize(x)) {
+            int n = ret.size();
+            for (int i = 0; i < n; i++) {
+                for (T a = 1, d = 1; d <= p.second; d++) {
+                    a *= p.first;
+                    ret.push_back(ret[i] * a);
+                }
+            }
+        }
+        return ret; // NOT sorted
+    }
+    // Moebius function Table, (-1)^{# of different prime factors} for square-free x
+    // return: [0=>0, 1=>1, 2=>-1, 3=>-1, 4=>0, 5=>-1, 6=>1, 7=>-1, 8=>0, ...] <https://oeis.org/A008683>
+    std::vector<int> GenerateMoebiusFunctionTable() {
+        std::vector<int> ret(min_factor.size());
+        for (unsigned i = 1; i < min_factor.size(); i++) {
+            if (i == 1)
+                ret[i] = 1;
+            else if ((i / min_factor[i]) % min_factor[i] == 0)
+                ret[i] = 0;
+            else
+                ret[i] = -ret[i / min_factor[i]];
+        }
+        return ret;
+    }
+};
+// Sieve sieve(1 << 15);  // (can factorize n <= 10^9)

number/test/enumerate_primes.test.cpp

@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/enumerate_primes"
+#include "../sieve.hpp"
+#include <iostream>
+#include <vector>
+using namespace std;
+
+int main() {
+    int N, A, B;
+    cin >> N >> A >> B;
+    Sieve sieve(N);
+    vector<int> ret;
+    for (unsigned i = 0; A * i + B < sieve.primes.size(); i++) ret.push_back(sieve.primes[A * i + B]);
+    cout << sieve.primes.size() << ' ' << ret.size() << '\n';
+    for (auto p : ret) cout << p << ' ';
+    cout << '\n';
+}

number/test/gen_primes.test.cpp

@@ -1,16 +1,15 @@
-#include "number/eratosthenes.hpp"
-#include <algorithm>
+#include "../sieve.hpp"
 #include <iostream>
 #include <vector>
 #define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ALDS1_1_C"
 using namespace std;
 
 int main() {
-    SieveOfEratosthenes sieve(10000);
+    Sieve sieve(10000);
     int N;
     cin >> N;
     int ret = 0;
-    for (int i = 0; i < N; i++) {
+    while (N--) {
         int x;
         cin >> x;
         bool flg = true;
@@ -21,5 +20,5 @@ int main() {
             }
         ret += flg;
     }
-    cout << ret << endl;
+    cout << ret << '\n';
 }

number/test/prime_factorization.test.cpp

@@ -1,16 +1,16 @@
-#include "number/eratosthenes.hpp"
+#include "../sieve.hpp"
 #include <iostream>
 #define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=NTL_1_A"
 using namespace std;
 
 int main() {
-    SieveOfEratosthenes sieve(40000);
+    Sieve sieve(1 << 15);
     int N;
     cin >> N;
-    map<long long int, int> factors = sieve.Factorize(N);
-    cout << N << ":";
+    map<long long int, int> factors = sieve.factorize<long long>(N);
+    cout << N << ':';
     for (auto p : factors) {
-        while (p.second--) cout << " " << p.first;
+        while (p.second--) cout << ' ' << p.first;
     }
-    cout << endl;
+    cout << '\n';
 }

utilities/segtree_relevant_points.hpp

@@ -1,3 +1,4 @@
+#pragma once
 #include <algorithm>
 #include <vector>