Prime Sieving

link

---

Tags: Trial Division, Sundaram Sieve, Erastothenes Sieve, Atkin Sieve, Wheel Sieve, Bitwise Sieve, Segment Bitwise Sieve, Segment Sieve, Sieve

---

0. Lời dẫn đầu

Note này tổng hợp lại các thuật toán về sàng nguyên tố, sắp xếp dựa trên tốc độ của thuật, code và tối ưu hỏa bởi SPyofgame.

Note này nhằm tới những code quá dài và phức tạp. Mặc dù có những code rất khó hiểu.

Bạn có thể xem tốc độ chay của các thuật toán tại đây

Code không có bản quyền và cho phép thoải mái sử dụng nhưng xin đừng sửa nguồn (Original) nếu có

Bạn có thể đóng góp ý kiện tại mục comment. Mình rất vui lòng phản hồi và cập nhật blog

Trong tương lai xa mình sẽ:

• Thêm hướng dẫn đầy đủ, cụ thể và chứng minh cho từng thuật sàng
• Thêm comment cho code
• Chia làm 2 note tiếng Việt với tiếng Anh
• Thêm một số loại sàng khác
• Thêm benchmark chó việc lấy số nguyên tố

---

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+\cdots +n)=O(n^2)$
Space:

$O(n)$
Algo: Brute-forces Sieve

#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}\le \frac{d+k}{2}\iff \sqrt{n}\le \frac{d+k}{2}\iff 1<d\le \sqrt{n}\le \frac{d+k}{2}\le k<n$

Từ đó suy ra ta chỉ cần duyêt qua

$1<d\le \sqrt{n}$

Time:

$O(\sqrt{1}+\sqrt{2}+\cdots +\sqrt{n})=O(\frac{1}{6}\lfloor \sqrt{n}\rfloor (5+6n-3\lfloor \sqrt{n}\rfloor -2{\lfloor \sqrt{n}\rfloor }^2))=O(n\sqrt{n})$
Space:

$O(n)$
Algo: Sqrt Sieve

#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\ge 5$ và

$k,r\in \mathbb{N}$ và

$0\le 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)=\lfloor \frac{\lfloor \sqrt{n}\rfloor }{3}\rfloor$
Time:

$O(\stackrel{f(n)}{\underset{x=1}{\mathrm{\Sigma }}}f(x))=O(n\sqrt{n})$ total -

$O(f(x))=O(\sqrt{x})$ check
Space:

$O(n)$
Algo: Sqrt Sieve

#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

#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\le 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(\stackrel{n}{\underset{x=1}{\mathrm{\Sigma }}}\lfloor \frac{n}{x}\rfloor \right)=O(n\log n)$
Space:

$O(n)$
Algo: Divisor Sieve

#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

#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\le n}{\mathrm{\Sigma }}\lfloor \frac{n}{p}\rfloor \right)=O(n\log \log n)$
Space:

$O(n)$
Algo: Eratosthenes Sieve

#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(\underset{\text{prime }p\le n}{\mathrm{\Sigma }}\lfloor \frac{n}{p}\rfloor -p)=O(n\log \log \sqrt{n})$
Space:

$O(n)$
Algo: Eratosthenes Sieve

#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]\le 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\le 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

#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

#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 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

#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

#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

#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

#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

#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

#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

#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

#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

#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

#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

#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

#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

#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

#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;
}