---
tags: Template, Trial Division, Sundaram Sieve, Erastothenes Sieve, Atkin Sieve, Wheel Sieve, Bitwise Sieve, Segment Bitwise Sieve, Segment Sieve, Sieve, SPyofgame
---
Prime Sieving
===
[link](link)
-----
###### Tags: `Trial Division`, `Sundaram Sieve`, `Erastothenes Sieve`, `Atkin Sieve`, `Wheel Sieve`, `Bitwise Sieve`, `Segment Bitwise Sieve`, `Segment Sieve`, `Sieve`
-----
### 1. Brute-forces
Với mọi số nguyên $x$ trong đoạn $[0, n]$. Ta đếm số ước $d$ của $x$ trong đoạn $[0, x]$. Nếu $x$ có đúng $2$ ước thì $x$ là số nguyên tố
> **Time:** $O(1 + 2 + \dots + n) = O(n^2)$
> **Space:** $O(n)$
> **Algo:** Brute-forces Sieve
```cpp=
#include <iostream>
using namespace std;
const int LIM = 1e6;
bool checkPrime(int n)
{
int cnt = 0;
for (int d = 1; d <= n; ++d)
if (n % d == 0)
++cnt;
return cnt == 2;
}
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
cntPrime = 0;
for (int x = 0; x <= n; ++x)
{
if (checkPrime(x))
{
++cntPrime;
isPrime[x] = true;
}
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 2.Trial Division
Nhận thấy là $n$ là số nguyên tố, khi và chỉ khi nó có đúng $2$ ước nguyên dương là $1$ và $n$
Nếu $n > 1$ và tồn tai một giá trị $d$ trong đoạn $[2, n]$ mà $n$ chia hết thì $n$ sẽ là hợp số
Thấy rằng khi tồn tại giá trị $d$ như thế, thì tồn tại số nguyên $k$ để $n = d \times k$
Không mất tính tổng quát ta giả sử $1 < d < k < n$
Áp dụng bất đẳng thức AM-GM cho 2 số nguyên dương $d$ và $k$ ta có $\sqrt{dk} \leq \frac{d + k}{2} \Leftrightarrow \sqrt{n} \leq \frac{d + k}{2} \Leftrightarrow 1 < d \leq \sqrt{n} \leq \frac{d + k}{2} \leq k < n$
Từ đó suy ra ta chỉ cần duyêt qua $1 < d \leq \sqrt{n}$
> **Time:** $O(\sqrt{1} + \sqrt{2} + \dots + \sqrt{n}) = O(\frac{1}{6} \left \lfloor \sqrt{n} \right \rfloor (5 + 6 n - 3\left \lfloor \sqrt{n} \right \rfloor - 2 \left \lfloor \sqrt{n} \right \rfloor^2)) = O(n \sqrt{n})$
> **Space:** $O(n)$
> **Algo:** Sqrt Sieve
```cpp=
#include <iostream>
using namespace std;
const int LIM = 1e6;
bool checkPrime(int n)
{
if (n <= 1) return false;
for (int d = 2; d * d <= n; ++d)
if (n % d == 0)
return false;
return true;
}
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
cntPrime = 0;
for (int x = 0; x <= n; ++x)
{
if (checkPrime(x))
{
++cntPrime;
isPrime[x] = true;
}
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 3.Improved Trial Division
Nhận xét các số có dạng $n = 6k + r$ với $n \geq 5$ và $k, r \in \mathbb{N}$ và $0 \leq r < 6$
- Khi $r \in \{0, 2, 4\}$ thì $n$ chia hết cho $2 \Rightarrow n$ không phải số nguyên tố.
- Khi $r \in \{0, 3\}$ thì $n$ chia hết cho $3 \Rightarrow n$ không phải số nguyên tố.
- Vậy ta chỉ cần kiểm tra các số thuộc dạng $n = 6k + 1$ và $n = 6k + 5$
Ta áp dụng tính chất trên cho cả việc sàng lẫn việc kiểm tra
> For $f(n) = \left \lfloor \frac{\left \lfloor \sqrt{n} \right \rfloor}{3} \right \rfloor$
> **Time:** $O\left(\overset{f(n)}{\underset{x = 1}{\Large \Sigma}} f(x) \right) = O(n \sqrt{n})$ total - $O(f(x)) = O(\sqrt{x})$ check
> **Space:** $O(n)$
> **Algo:** Sqrt Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
bool checkPrime(int n)
{
if (n <= 1) return false;
if (n == 2 || n == 3 || n == 5) return true;
if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;
for (int d = 5; d * d <= n; d += 6)
{
if (n % d == 0) return false;
if (n % (d + 2) == 0) return false;
}
return true;
}
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, false, sizeof(isPrime));
isPrime[2] = isPrime[3] = true;
cntPrime = (n >= 2) + (n >= 3);
for (int x = 5; x <= n; x += 6)
{
if (checkPrime(x))
{
++cntPrime;
isPrime[x] = true;
}
if (checkPrime(x + 2))
{
++cntPrime;
isPrime[x + 2] = true;
}
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 4.Second Improved Trial Division
Ta nhận thấy nếu $x$ không phải số nguyên tố thì tồn tại một số $1 < d < X$ mà $X$ chia hết.
Mà mọi ước của $d$ đều là ước của $x$
Nên ta chỉ cần kiểm tra xem $x$ có ước nguyên tố nào hay không
> **Time:** $O\left(\frac{n\sqrt{n}}{\ln n}\right)$
> **Space:** $O(n)$
> **Algo:** Sqrt Sieve
```cpp
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
typedef unsigned int uint;
const int LIM = 1e6;
int cntPrime;
int prime[LIM + 1];
void sieve(int n)
{
cntPrime = 0;
if (n < 2) return ;
prime[cntPrime++] = 2;
if (n < 3) return ;
prime[cntPrime++] = 3;
int sqrtn = 3, tickn = 1;
for (int x = 5, t = 1; x <= n; x += 2 + 2 * (t ^= 1))
{
bool ok = true;
for (int i = 0; i < cntPrime && prime[i] <= sqrtn; ++i)
{
if (x % prime[i] == 0)
{
ok = false;
break;
}
}
if (ok) prime[cntPrime++] = x;
if (--tickn == 0) tickn = sqrtn++;
if (t && --tickn == 0) tickn = sqrtn++;
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 5.Divisor Sieve
Với mỗi số $x$. Đánh dấu các bội của $x$ khác $x$ gồm các số $2x, 3x, \dots, kx$ ($kx \leq n$) không phải số nguyên tố. Chứng minh:
- Vì các số có dạng $kx$ và $k > 1$ nên có ít nhất $3$ ước là $1, x, k$ nên không phải số nguyên tố
> **Time:** $O\left(\overset{n}{\underset{x = 1}{\Large \Sigma}} \left \lfloor \frac{n}{x} \right \rfloor \right) = O(n\log n)$
> **Space:** $O(n)$
> **Algo:** Divisor Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int p = 2; p * p <= n; ++p)
for (int v = 2; p * v <= n; ++v)
isPrime[p * v] = false;
cntPrime = 0;
for (int p = 0; p <= n; ++p) cntPrime += isPrime[p];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 6.Sundaram Sieve
> **Time:** $O(n\log n)$
> **Space:** $O(n)$
> **Algo:** Divisor Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM / 2 + 1];
void sieve(int n)
{
int k = (n - 2) / 2;
memset(isPrime, true, sizeof(isPrime[0]) * (k + 1));
isPrime[0] = false;
for (int i = 1; i * i <= n; ++i)
for (int j = i; i + j + 2 * i * j <= k; ++j)
isPrime[i + j + 2 * i * j] = false;
cntPrime = 1;
for (int i = 1; i <= k; ++i)
if (isPrime[i] == true) /// prime 2 * i + 1
++cntPrime;
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 7.Eratosthenes Sieve
Nhận thấy rằng các số có nhiều ước thì sẽ bị đánh dấu nhiều lần. Vì khi $x$ là hợp số, thì với mọi ước $v$ của $x$, thì bội $kx$ của $x$ cũng là bội của $v$. Nên để hạn chế điều đó, ta chỉ đánh dấu các bội của $x$ khi $x$ là số nguyên tố.
> **Time:** $O\left(\underset{\text{prime }p \leq n}{\Large \Sigma} \left \lfloor \frac{n}{p} \right \rfloor \right) = O(n\log\log n)$
> **Space:** $O(n)$
> **Algo:** Eratosthenes Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int p = 2; p * p <= n; ++p) if (isPrime[p])
for (int x = p + p; x <= n; x += p)
isPrime[x] = false;
cntPrime = 0;
for (int p = 0; p <= n; ++p) cntPrime += isPrime[p];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 8.Improved Eratosthenes Sieve
Mình có thể cải thiện hơn nữa. Bằng cách đánh dấu các số chẵn khác $2$ không phải số nguyên tố. Rồi xét các số lẻ
Thay vì đánh dấu các bội $kx$ với $k > 1$, thì nhận thấy mọi $kx$ với $k < x$ đều đã bị đánh dấu bởi một số nguyên tố nào đó trước $x$. Nên ta đánh dấu từ $k = x$ luôn.
> **Time:** $O\left(\underset{\text{prime }p \leq n}{\Large \Sigma} \left \lfloor \frac{n}{p} \right \rfloor - p \right) = O(n\log\log \sqrt{n})$
> **Space:** $O(n)$
> **Algo:** Eratosthenes Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int p = 4; p <= n; p += 2) isPrime[p] = false;
for (int p = 3; p * p <= n; p += 2) if (isPrime[p])
for (int x = p * p; x <= n; x += p)
isPrime[x] = false;
cntPrime = (n >= 2);
for (int p = 3; p <= n; p += 2) cntPrime += isPrime[p];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 9.Linear Eratosthenes Sieve
Nhưng, với cách trên thì mình vẫn phải đánh dấu nhiều lần các số có nhiều ước nguyên tố.
Ta gọi `lpf[x]` là ước nguyên tố nhỏ nhất của $x$ (viết tắt của `lowest_prime_factor`)
Khi này mọi số $x > 1$ đều có thể được biểu diễn bằng $x = k \times lpf[x]$. Mà từ định nghĩa ta có $lpf[x] \leq lpf[k]$
Giờ nếu ta chỉ đánh dấu các số có dạng $kx$ với $k$ nguyên tố và $k \leq lpf[x]$. Thì mỗi số chỉ bị đánh dấu bởi một số nguyên tố mà thôi
> **Time:** $O(n)$
> **Space:** $O(n)$ lpf[] + $O\left(\frac{n}{ln(n)} \right)$ prime[]
> **Algo:** Eratosthenes Sieve
```cpp=
#include <iostream>
#include <vector>
using namespace std;
const int LIM = 1e6;
vector<int> prime;
vector<int> lpf;
void sieve(int n)
{
prime.assign(1, 2);
lpf.assign(n + 1, 2);
for (int x = 3; x <= n; x += 2)
{
if (lpf[x] == 2) prime.push_back(lpf[x] = x);
for (int i = 0; i < prime.size() && prime[i] <= lpf[x] && prime[i] * x <= n; ++i)
lpf[prime[i] * x] = prime[i];
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << prime.size();
return 0;
}
```
-----
### 10.Improved Linear Eratosthenes Sieve
Mình biết trước được số lượng phần tử rồi thì có thể sài mảng tĩnh để xử lí nhanh hơn
> **Time:** $O(n)$
> **Space:** $O(n)$ lpf[] + $O\left(\frac{n}{ln(n)} \right)$ prime[]
> **Algo:** Eratosthenes Sieve
```cpp=
#include <iostream>
#include <vector>
using namespace std;
const int LIM = 1e6;
const int LIM_P = 664579 + 66;
int cntPrime;
int prime[LIM_P];
int lpf[LIM + 1];
void sieve(int n)
{
cntPrime = 0;
prime[cntPrime++] = 2;
fill_n(lpf, n + 1, 2);
for (int x = 3; x <= n; x += 2)
{
if (lpf[x] == 2) prime[cntPrime++] = (lpf[x] = x);
for (int i = 0; i < cntPrime && prime[i] <= lpf[x] && prime[i] * x <= n; ++i)
lpf[prime[i] * x] = prime[i];
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
---
-----
### 11.Atkin Sieve
Cái thuật sàng này thì phức tạp hơn nhiều và cũng sử dụng nhiều phép toán nhân chia nên khó tối ưu hơn là sàng Erastothenes.
Nhưng mà thôi kệ, cứ viết về nó vậy.
[Ý tưởng](https://www.ams.org/journals/mcom/2004-73-246/S0025-5718-03-01501-1/S0025-5718-03-01501-1.pdf) của **A. O. L. Atkin** và **Daniel J. Bernstein** là như sau:
- Một số nguyên dương không chính phương $p \in 1 + 4k$ là số nguyên tố khi và chỉ khi $4x^2 + y^2 = p$ có lẻ cặp nghiệm nguyên dương $(x, y)$
Để phủ hết các số nguyên tố ta kiểm tra thêm các dạng
- $p \in 1 + 4k$ xét phương trình $4x^2 + y^2 = p$ với $x, y > 0$
- $p \in 7 + 12k$ xét phương trình $3x^2 + y^2 = p$ với $x, y > 0$
- $p \in 11 + 12k$ xét phương trình $3x^2 - y^2 = p$ với $x > y > 0$
- $\dots$
> **Time:** $O(n)$
> **Space:** $O(n)$
> **Algo:** Atkin Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
const int LIM = 1e6;
int cntPrime = 0;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, false, sizeof(isPrime));
if (n < 2) return ;
isPrime[2] = true; cntPrime = 1;
if (n < 3) return ;
isPrime[3] = true; cntPrime = 2;
if (n < 4) return ;
for (int x = 1; x * x <= n; ++x)
{
int A = x * x;
for (int y = 1; y * y <= n - A; ++y)
{
int B = y * y, p;
if (p = 4 * A + B, p <= n && (p % 12 == 1 || p % 12 == 5)) isPrime[p] = !isPrime[p];
if (p = 3 * A + B, p <= n && p % 12 == 7) isPrime[p] = !isPrime[p];
if (p = 3 * A - B, p <= n && p % 12 == 11 && x > y) isPrime[p] = !isPrime[p];
}
}
int rn = sqrt(n);
for (int x = 5; x <= rn; ++x) if (isPrime[x])
for (int k = 1; k * x * x <= n; ++k)
isPrime[k * x * x] = false;
cntPrime = count(isPrime, isPrime + n + 1, 1);
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 12.Improved Atkin Sieve
Ta chia làm các vòng for độc lập để có thể có những tối ưu riêng từng cái
> **Time:** $O(n)$
> **Space:** $O(n)$
> **Algo:** Atkin Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, false, sizeof(isPrime));
if (n < 2) return ;
isPrime[2] = true; cntPrime = 1;
if (n < 3) return ;
isPrime[3] = true; cntPrime = 2;
if (n < 4) return ;
for (int upx = sqrt(n / 4) + 1, x = 1; x <= upx; ++x)
{
int X = 4 * x * x;
for (int upy = sqrt(n - X), y = 1; y <= upy; ++y)
{
int p = X + y * y;
if (p % 12 == 1 || p % 12 == 5)
isPrime[p] = !isPrime[p];
}
}
for (int upx = sqrt(n / 3) + 1, x = 1; x <= upx; x += 2)
{
int t = 1;
int X = 3 * x * x;
for (int upy = sqrt(n - X), y = 2; y <= upy; y += 6)
{
int p = X + y * y;
isPrime[p] = !isPrime[p];
}
for (int upy = sqrt(n - X), y = 4; y <= upy; y += 6)
{
int p = X + y * y;
isPrime[p] = !isPrime[p];
}
}
for (int rn = sqrt(n), x = 1; x <= rn; ++x)
{
int X = 3 * x * x;
for (int y = max(int(sqrt(X - n)) + 1, 1); y < x; ++y)
{
int p = X - y * y;
if (p % 12 == 11)
isPrime[p] = !isPrime[p];
}
}
int rn = sqrt(n);
for (int i = 5; i <= rn; i += 2) if (isPrime[i])
for (int t = i * i, j = t; j < n; j += t)
isPrime[j] = false;
for (int i = 5; i <= n; i += 2)
if (isPrime[i])
++cntPrime;
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 13.Non-modulo Atkin Sieve
Áp dụng toán học và tìm những trường hợp nào nghiệm thỏa cả phương trình lẫn điều kiện modulo và đơn giản hóa lại thành phép cộng trừ ta thu được như sau
> **Time:** $O(n)$
> **Space:** $O(n)$
> **Algo:** Atkin Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
constexpr int LIM = 1e6;
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, false, sizeof(isPrime));
if (n < 2) return ;
isPrime[2] = true; cntPrime = 1;
if (n < 3) return ;
isPrime[3] = true; cntPrime = 2;
if (n < 4) return ;
for (int upx = sqrt(n / 4) + 1, x = 1; x <= upx; ++x)
{
int X = 4 * x * x;
int t = 0;
int b = (x % 3 == 0);
for (int y = 1, p = X + y * y; p <= n; y += 2 + (b ? 2 * (t ^= 1) : 0), p = X + y * y)
isPrime[p] = !isPrime[p];
}
for (int upx = sqrt(n / 3) + 1, x = 1; x <= upx; x += 2)
{
int t = 1;
int X = 3 * x * x;
for (int y = 2, p = X + y * y; p <= n; y += 2 + 2 * (t ^= 1), p = X + y * y)
isPrime[p] = !isPrime[p];
}
for (int rn = sqrt(n), x = 1; x <= rn; ++x)
{
int X = 3 * x * x;
int t = x & 1;
for (int y = 1 + t, p = X - y * y; y < x; y += 2 + 2 * (t ^= 1))
{
p = X - y * y;
if (p <= n) isPrime[p] = !isPrime[p];
}
}
int rn = sqrt(n);
for (int i = 5; i <= rn; i += 2) if (isPrime[i])
for (int t = i * i, j = t; j < n; j += t)
isPrime[j] = false;
for (int i = 5; i <= n; i += 2)
if (isPrime[i])
++cntPrime;
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 14.Improved Non-modulo Atkin Sieve
> **Time:** $O(n)$
> **Space:** $O(n)$
> **Algo:** Atkin Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
constexpr int LIM = 1e6;
constexpr int offset[2] = {2, 4};
constexpr int n3[3] = {1, 2, 0};
int cntPrime;
bool isPrime[LIM + 1];
void sieve(int n)
{
memset(isPrime, false, sizeof(isPrime));
if (n < 2) return ;
isPrime[2] = true; cntPrime = 1;
if (n < 3) return ;
isPrime[3] = true; cntPrime = 2;
if (n < 4) return ;
for (int upx = sqrt(n / 4) + 1, x = 1, x3 = 1, X = 4; x <= upx; ++x, x3 = n3[x3], X = 4 * x * x)
{
bool t = false;
for (int y = 1, p = X + y * y; p <= n; y += (x3 ? 2 : offset[t ^= 1]), p = X + y * y)
isPrime[p] = !isPrime[p];
}
for (int upx = sqrt(n / 3) + 1, x = 1, X = 3; x <= upx; x += 2, X = 3 * x * x)
{
bool t = true;
for (int y = 2, p = X + y * y; p <= n; y += offset[t ^= 1], p = X + y * y)
isPrime[p] = !isPrime[p];
}
// This code is slightly faster
for (int rn = sqrt(n), x = 1, X = 3, v = 2, T = 1; x <= rn; ++x, X = 3 * x * x)
{
while (T * T <= X - n) ++T;
bool t = x & 1;
int y = 1 + t;
if (X - y * y > n) y += (T - y) / 6 * 6;
while (X - y * y > n) y += offset[t ^= 1];
for (int p = X - y * y; y < x; y += offset[t ^= 1], p = X - y * y)
isPrime[p] = !isPrime[p];
}
bool t = 1;
int rn = sqrt(n);
for (int x = 5; x <= rn; x += offset[t ^= 1]) if (isPrime[x])
for (int t = x * x, p = t; p <= n; p += t)
isPrime[p] = false;
for (int x = 5; x <= n; x += 6)
cntPrime += isPrime[x] + isPrime[x + 2]; /// x is Prime or/and (x + 2) is Prime
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 15.Eratosthenes Wheel-2 Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(n)$
> **Algo:** Erastothenes Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM / 2 + 1];
void sieve(int n)
{
memset(isPrime, true, sizeof(isPrime[0]) * (n / 2 + 1));
isPrime[0 / 2] = isPrime[1 / 2] = false;
for (int p = 2; p * p <= n; ++p) if (isPrime[p / 2])
for (int x = p + p; x <= n; x += p) if (x % 2 == 1)
isPrime[x / 2] = false;
cntPrime = (n >= 2);
for (int p = 3; p <= n; p += 2) if (isPrime[p / 2]) cntPrime += isPrime[p / 2];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 16.Improved Eratosthenes Wheel-2 Sieve
> **Time:** $O(n \log \log \sqrt{n})$
> **Space:** $O(n)$
> **Algo:** Erastothenes Sieve
```cpp=
#include <iostream>
#include <cstring>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bool isPrime[LIM / 2 + 1];
void sieve(int n)
{
memset(isPrime, true, sizeof(isPrime[0]) * (n / 2 + 1));
isPrime[0] = false;
for (int p = 3; p * p <= n; p += 2) if (isPrime[p >> 1])
for (int x = p + p + p; x <= n; x += 2 * p)
isPrime[x >> 1] = false;
cntPrime = (n >= 2) + (n >= 3);
for (int p = 5; p <= n; p += 6) cntPrime += isPrime[p >> 1] + isPrime[(p + 2) >> 1];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 17.Eratosthenes Wheel-2 Bitwise Sieve
> **Time:** $O(n \log \log \sqrt{n})$
> **Space:** $O(n)$
> **Algo:** Erastothenes Sieve
```cpp=
#include <iostream>
#include <bitset>
using namespace std;
const int LIM = 1e6;
int cntPrime;
bitset<LIM / 2 + 1> isNotPrime;
void sieve(int n)
{
isNotPrime[0] = true;
for (int p = 3; p * p <= n; p += 2) if (!isNotPrime[p >> 1])
for (int x = p * p; x <= n; x += 2 * p)
isNotPrime[x >> 1] = true;
cntPrime = (n >= 2);
for (int p = 3; p <= n; p += 2) cntPrime += !isNotPrime[p >> 1];
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 18.Improved Eratosthenes Wheel-2 Bitwise Sieve
> **Time:** $O(n \log \log \sqrt{n})$
> **Space:** $O(n)$
> **Algo:** Erastothenes Sieve
```cpp=
#include <iostream>
#include <bitset>
using namespace std;
typedef unsigned int uint;
const uint LIM = 1e6;
uint cntPrime;
uint isNotPrime[(LIM >> 5) / 2 + 3];
uint get(uint p) { return (isNotPrime[p >> 5] >> (p & 31)) & 1; }
void set(uint p) { isNotPrime[p >> 5] |= 1 << (p & 31); }
void sieve(uint n)
{
isNotPrime[0] = true;
for (uint p = 3; p * p <= n; p += 2) if (!get(p >> 1))
for (uint x = p * p; x <= n; x += 2 * p)
set(x >> 1);
cntPrime = (n >= 2) + (n >= 3);
for (uint p = 5; p <= n; p += 6) cntPrime += !get(p >> 1) + !get((p + 2) >> 1);
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 19.Improved Erastothenes Wheel-2-3-5-7 Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
#define all(x) (x).begin(), (x).end()
const int LIM = 100000000;
bool isPrime[(LIM + 1 - 7) / 2];
typedef unsigned int uint;
constexpr uint WHEEL_1[] = {2, 1, 2, 1, 2, 3, 1, 3};
constexpr uint WHEEL_0[] = {1, 0, 1, 0, 1, 2, 0, 2};
constexpr uint NEXT8[] = {1, 2, 3, 4, 5, 6, 7, 0};
uint cntPrime;
void sieve(uint n)
{
cntPrime = 0u; if (n < 2u) return ;
cntPrime = 1u; if (n < 3u) return ;
cntPrime = 2u; if (n < 5u) return ;
cntPrime = 3u; if (n < 7u) return ;
uint length = (n - 7u) / 2u;
uint sqrt_lim = (sqrt(n) - 7u) / 2u;
memset(isPrime, true, sizeof(isPrime[0]) * (length + 1));
for (uint i = 0u, w = 0u; i <= sqrt_lim; i += WHEEL_1[w], w = NEXT8[w]) if (isPrime[i])
{
uint p = 7u + i + i;
uint pa[] = {p, p + p, p + p + p};
for (uint j = (p * p - 7) / 2, m = w; j <= length; j += pa[WHEEL_0[m]], m = NEXT8[m])
isPrime[j] = false;
}
for (uint i = 0, w = 0; i <= length; i += WHEEL_1[w], w = NEXT8[w]) if (isPrime[i]) ++cntPrime; /// 7 + i + i is Prime
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 20.Improved Erastothenes Wheel-2-3-5-7 Segment Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
typedef unsigned int uint;
constexpr uint LIM = 1e6;
constexpr uint NEXT8[] = {1, 2, 3, 4, 5, 6, 7, 0};
constexpr uint WHEEL_1[] = {2, 1, 2, 1, 2, 3, 1, 3};
constexpr uint WHEEL_0[] = {1, 0, 1, 0, 1, 2, 0, 2};
constexpr uint WHEEL_POS[] = {0, 2, 3, 5, 6, 8, 11, 12};
constexpr uint WHEEL_IDX[] = {0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7, 7, 7, 0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7};
constexpr uint WHEEL_DUP[] = {0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, 15, 15, 17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27};
constexpr uint L1CACHEPOW = 14u + 3u;
constexpr uint L1CACHESZ = (1u << (int)L1CACHEPOW);
constexpr uint SEGSZ = L1CACHESZ / 15u * 15u;
constexpr uint PRESZ = (65535u - 7u) / 2u;
bool isPrime[PRESZ + 1];
bool segPrime[SEGSZ];
uint cntPrime;
void sieve(uint n)
{
cntPrime = 0u; if (n < 2u) return;
cntPrime = 1u; if (n < 3u) return;
cntPrime = 2u; if (n < 5u) return;
cntPrime = 3u; if (n < 7u) return;
const uint upper = (n - 7u) / 2u;
uint sqrt_upper = (sqrt(n) - 7u) / 2u;
memset(isPrime, true, sizeof(isPrime));
for (uint i = 0u; i <= 124u; ++i) if (isPrime[i])
{
uint p = 7u + i + i;
for (uint j = (p * p - 7u) / 2u; j <= PRESZ; j += p)
isPrime[j] = false;
}
for (uint i = 0u, w = 0u, si = 0u; i <= upper; i += WHEEL_1[w], si += WHEEL_1[w], si = (si >= SEGSZ) ? 0u : si, w = NEXT8[w])
{
if (si == 0u)
{
memset(segPrime, true, sizeof(segPrime));
for (uint j = 0u, bw = 0u; j <= PRESZ; j += WHEEL_1[bw], bw = NEXT8[bw]) if (isPrime[j])
{
uint p = 7u + j + j;
uint k = p * (j + 3u) + j;
if (k >= i + SEGSZ) break;
uint pa[3] = {p, p + p, p + p + p};
uint sw = bw;
if (k < i)
{
k = (i - k);
k %= (15u * p);
if (k != 0)
{
uint os = WHEEL_POS[bw];
sw = os + ((k + p - 1u) / p);
sw = WHEEL_DUP[sw];
k = (sw - os) * p - k;
sw = WHEEL_IDX[sw];
}
}
else k -= i;
for (; k < SEGSZ; k += pa[WHEEL_0[sw]], sw = NEXT8[sw])
segPrime[k] = false;
}
}
if (segPrime[si]) ++cntPrime; /// 7 + i + i is Prime
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 21.Erastothenes Segment Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
```cpp=
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
const int LIM = 1e6;
const int SQRT_LIM = sqrt(LIM);
bool isPrime[SQRT_LIM + 1];
int cntPrime;
void sieve(int n)
{
int LIM = floor(sqrt(n));
vector<int> prime;
memset(isPrime, true, sizeof(isPrime[0]) * LIM);
for (int p = 2; p * p <= LIM; ++p) if (isPrime[p] == true)
for (int x = p + p; x <= LIM; x += p)
isPrime[x] = false;
cntPrime = 0;
for (int p = 2; p <= LIM; ++p) if (isPrime[p])
{
++cntPrime;
prime.push_back(p);
}
for (int low = LIM; low <= n; low += LIM)
{
int high = min(n, low + LIM);
memset(isPrime, true, sizeof(isPrime[0]) * LIM);
for (int i = 0; i < prime.size(); ++i)
{
int p = (low + prime[i] - 1) / prime[i] * prime[i];
for (int x = p; x <= high; x += prime[i])
isPrime[x - low] = false;
}
for (int n = low; n < high; ++n)
if (isPrime[n - low] == true)
++cntPrime;
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 22.Improved Erastothenes Segment Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
#define all(x) (x).begin(), (x).end()
const int LIM = 1e6;
const int OPTIMAL_RANGE = 1 << 15;
const int offset[2] = {2, 4};
int cntPrime;
bool mark[OPTIMAL_RANGE];
bool isPrime[OPTIMAL_RANGE];
int seg_prime[OPTIMAL_RANGE];
int seg_multi[OPTIMAL_RANGE];
void sieve(int lim)
{
if (lim < 2) return ;
int sqrt = std::sqrt(lim);
int segment_size = max(sqrt, OPTIMAL_RANGE);
memset(isPrime, true, sizeof(isPrime));
cntPrime = (lim >= 2) + (lim >= 3);
int qdrt = std::sqrt(sqrt);
for (int x = 5, t = 0; x <= qdrt; x += offset[t ^= 1]) if (isPrime[x])
for (int p = x * x; p <= sqrt; p += x)
isPrime[p] = false;
int n_seg = 0;
int s = 5, n = 5;
bool tn = 1, ts = 1;
for (int low = 0; low <= lim; low += segment_size)
{
int high = min(low + segment_size - 1, lim);
memset(mark, true, sizeof(mark[0]) * (high - low + 1));
for (; s * s <= high; s += offset[ts ^= 1]) if (isPrime[s])
{
seg_prime[n_seg] = s;
seg_multi[n_seg] = s * s - low;
++n_seg;
}
for (int x = 0; x < n_seg; ++x)
{
int p = seg_multi[x];
for (int k = seg_prime[x] * 2; p < segment_size; p += k) mark[p] = false;
seg_multi[x] = p - segment_size;
}
for (; n <= high; n += offset[tn ^= 1])
if (mark[n - low])
++cntPrime;
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 23.Erastothenes Segment Bitwise Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
> **Original:** [primesieve](https://github.com/kimwalisch/primesieve)
```cpp=
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
typedef unsigned int uint;
constexpr uint LIM = 1e6;
constexpr uint SQRT_LIM = sqrt(LIM);
uint isPrime[(SQRT_LIM >> 5) + 1];
uint is_prime(uint p) { return (isPrime[p >> 5] >> (p & 31)) & 1; }
void unset(uint p) { isPrime[p >> 5] &= ~(1 << (p & 31)); }
int cntPrime;
void sieve(uint n)
{
const uint LIM = floor(sqrt(n));
const uint RANGE = sizeof(isPrime[0]) * ((LIM >> 5) + 1);
memset(isPrime, -1, RANGE);
for (uint p = 2; p * p <= LIM; ++p) if (is_prime(p))
for (int x = p + p; x <= LIM; x += p)
unset(x);
vector<uint> prime;
for (int p = 2; p <= LIM; ++p)
if (is_prime(p))
prime.push_back(p);
cntPrime = prime.size();
for (uint low = LIM; low <= n; low += LIM)
{
uint high = min(n + 1, low + LIM);
memset(isPrime, -1, RANGE);
for (uint i = 0; i < prime.size(); ++i)
for (uint x = prime[i], p = (low + x - 1) / x * x; p <= high; p += x)
unset(p - low);
for (uint n = low; n < high; ++n)
if (is_prime(n - low))
++cntPrime;
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
### 24.Improved Erastothenes Segment Bitwise Sieve
> **Time:** $O(n \log \log n)$
> **Space:** $O(\sqrt{n})$
> **Algo:** Erastothenes Sieve
> **Original:** [primesieve](https://github.com/kimwalisch/primesieve)
```cpp=
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
typedef unsigned int uint;
constexpr uint unset_bit[8] = { uint(-2), uint(-2), uint(-3), uint(-3), uint(-5), uint(-5), uint(-9), uint(-9)};
constexpr uint popcnt[128] = {
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
};
constexpr uint LIM = 1e6;
constexpr uint SQRT_LIM = sqrt(LIM);
constexpr uint OPTIMAL_RANGE = 1 << 15;
constexpr uint offset[2] = {2, 4};
uint8_t mark[max(SQRT_LIM, OPTIMAL_RANGE)];
uint8_t is_prime[SQRT_LIM + 1];
uint cntPrime;
void sieve(uint lim)
{
cntPrime = 0;
if (lim <= 1) return ;
uint sqrt = std::sqrt(lim);
uint mark_size = std::max(sqrt, OPTIMAL_RANGE);
uint segment_size = mark_size * 8;
uint n = 3, s = 3;
bool tn = 1, ts = 1;
std::vector<uint> seg_prime;
std::vector<uint> seg_multi;
memset(is_prime, true, sizeof(is_prime));
for (uint low = 0; low <= lim; low += segment_size)
{
uint high = min(lim, low + segment_size - 1);
mark_size = ((high - low) >> 3) + 1;
fill_n(mark, mark_size, 0xf);
for (; n * n <= high; n += offset[tn ^= 1]) if (is_prime[n])
for (uint j = n * n; j <= sqrt; j += n)
is_prime[j] = false;
for (; s * s <= high; s += offset[ts ^= 1]) if (is_prime[s])
{
seg_prime.push_back(s);
seg_multi.push_back(s * s - low);
}
for (uint n = 0; n < seg_prime.size(); ++n)
{
uint j = seg_multi[n];
for (uint k = seg_prime[n] * 2; j < segment_size; j += k) mark[j >> 3] &= unset_bit[j & 7];
seg_multi[n] = j - segment_size;
}
if (high == lim) mark[mark_size - 1] &= ~(0xf << (((lim & 7) + 1) >> 1));
for (uint n = 0; n < mark_size; n++) cntPrime += popcnt[mark[n]]; /// prime at bit[n] = low + n * 16 + i * 2 + 1
}
}
int main()
{
int n;
cin >> n;
sieve(n);
cout << "Number of primes up to " << n << " is " << cntPrime;
return 0;
}
```
-----
# Bench-mark
[GNU G++17 7.3.0 CF 23/9/2021](https://codeforces.com/problemset/customtest)
| Num | Algorithm | Time Complexity | $N = 10^4$ | $N = 10^5$ | $N = 10^6$ | $N = 10^7$ | $N = 10^8$ | $N = 10^9$ |
| --: | ------------------------------------------------: | :-------------------------------------- | :--------: | :--------: | :--------: |----------: | :--------: | :--------: |
| 1. | Brute-forces | $O(n^2)$ | 109ms | 11215ms | >15000ms | >15000ms | >15000ms | MLE |
| 2. | Trial Division | $O(n\sqrt{n})$ | 15ms | 15ms | 171ms | 3993ms | >15000ms | MLE |
| 3. | Improved Trial Division | $O(n\sqrt{n})$ | 0ms | 0ms | 61ms | 1325ms | >15000ms | MLE |
| 4. | Second Improved Trial Division | $O\left(\frac{n\sqrt{n}}{\ln n}\right)$ | 0ms | 0ms | 61ms | 1325ms | >15000ms | MLE |
| 5. | Divisor Sieve | $O(n\log n)$ | 15ms | 31ms | 31ms | 576ms | 14617ms | MLE |
| 6. | Sundaram Sieve | $O(n\log n)$ | 15ms | 31ms | 31ms | 31ms | 1419ms | >15000ms |
| 7. | Eratosthenes Sieve | $O(n\log\log n)$ | 0ms | 15ms | 15ms | 78ms | 1169ms | MLE |
| 8. | Improved Eratosthenes Sieve | $O(n\log\log \sqrt{n})$ | 0ms | 15ms | 15ms | 93ms | 1014ms | MLE |
| 9. | Linear Eratosthenes Sieve | $O(n)$ | 0ms | 15ms | 15ms | 61ms | 795ms | MLE |
| 10. | Improved Linear Eratosthenes Sieve | $O(n)$ | 0ms | 0ms | 15ms | 77ms | 732ms | MLE |
| 11. | Atkin Sieve | $O(n)$ | 0ms | 15ms | 15ms | 93ms | 1122ms | MLE |
| 12. | Improved Atkin Sieve | $O(n)$ | 0ms | 0ms | 0ms | 61ms | 733ms | MLE |
| 13. | Non-modulo Atkin Sieve | $O(n)$ | 0ms | 15ms | 15ms | 46ms | 561ms | MLE |
| 14. | Improved Non-modulo Atkin Sieve | $O(n)$ | 0ms | 0ms | 0ms | 46ms | 402ms | MLE |
| 15. | Eratosthenes Wheel-2 Sieve | $O(n\log\log n)$ | 0ms | 0ms | 15ms | 62ms | 857ms | 9000ms |
| 16. | Improved Eratosthenes Wheel-2 Sieve | $O(n\log\log \sqrt{n})$ | 0ms | 0ms | 0ms | 15ms | 546ms | 6666ms |
| 17. | Eratosthenes Wheel-2 Bitwise Sieve | $O(n\log\log \sqrt{n})$ | 0ms | 0ms | 0ms | 30ms | 234ms | 4648ms |
| 18. | Improved Eratosthenes Wheel-2 Bitwise Sieve | $O(n\log\log \sqrt{n})$ | 0ms | 0ms | 0ms | 30ms | 218ms | 4071ms |
| 19. | Improved Erastothenes Wheel-2-3-5-7 Sieve | $O(n \log \log n)$ | 0ms | 15ms | 15ms | 30ms | 358ms | 4242ms |
| 20. | Improved Erastothenes Wheel-2-3-5-7 Segment Sieve | $O(n \log \log n)$ | 0ms | 0ms | 0ms | 30ms | 186ms | 1980ms |
| 21. | Erastothenes Segment Sieve | $O(n \log \log n)$ | 0ms | 15ms | 31ms | 46ms | 312ms | 2776ms |
| 22. | Improved Erastothenes Segment Sieve | $O(n \log \log n)$ | 0ms | 15ms | 15ms | 30ms | 109ms | 1466ms |
| 23. | Erastothenes Segment Bitwise Sieve | $O(n \log \log n)$ | 0ms | 0ms | 0ms | 62ms | 546ms ? | 5366ms ? |
| 24. | Improved Erastothenes Segment Bitwise Sieve | $O(n \log \log n)$ | 0ms | 0ms | 0ms | 15ms | 46ms | 592ms |
Sign in with Wallet
Connect another wallet