Prime Sieving

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+\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;
}

---

Bench-mark

GNU G++17 7.3.0 CF 23/9/2021

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