Thuật toán tìm số nguyên tố với sàng Atkin

def sieve_of_atkin(limit):
    # Mảng đánh dấu các số nguyên tố, ban đầu tất cả đều là False
    primes = [False] * (limit + 1)
    
    # Sàng các số nguyên tố bắt đầu từ 2 và 3
    if limit > 2:
        primes[2] = True
    if limit > 3:
        primes[3] = True

    # Lặp qua tất cả các cặp số (x, y)
    for x in range(1, int(limit**0.5) + 1):
        for y in range(1, int(limit**0.5) + 1):
            # Tính n = 4x^2 + y^2
            n = 4*x*x + y*y
            if n <= limit and (n % 12 == 1 or n % 12 == 5):
                primes[n] = not primes[n]

            # Tính n = 3x^2 + y^2
            n = 3*x*x + y*y
            if n <= limit and n % 12 == 7:
                primes[n] = not primes[n]

            # Tính n = 3x^2 - y^2
            if x > y:
                n = 3*x*x - y*y
                if n <= limit and n % 12 == 11:
                    primes[n] = not primes[n]

    # Loại bỏ bội số của các số nguyên tố lên đến căn bậc hai của limit
    for r in range(5, int(limit**0.5) + 1):
        if primes[r]:
            for i in range(r*r, limit + 1, r*r):
                primes[i] = False

    # In ra các số nguyên tố
    return [p for p, prime in enumerate(primes) if prime]

# Ví dụ sử dụng
limit = 20
print(sieve_of_atkin(limit))

#include <iostream>
#include <vector>
#include <cmath>

std::vector<int> sieve_of_atkin(int limit) {
    std::vector<bool> primes(limit + 1, false);

    if (limit > 2)
        primes[2] = true;
    if (limit > 3)
        primes[3] = true;

    for (int x = 1; x * x < limit; x++) {
        for (int y = 1; y * y < limit; y++) {
            int n = 4 * x * x + y * y;
            if (n <= limit && (n % 12 == 1 || n % 12 == 5))
                primes[n] = !primes[n];

            n = 3 * x * x + y * y;
            if (n <= limit && n % 12 == 7)
                primes[n] = !primes[n];

            if (x > y) {
                n = 3 * x * x - y * y;
                if (n <= limit && n % 12 == 11)
                    primes[n] = !primes[n];
            }
        }
    }

    for (int r = 5; r * r < limit; r++) {
        if (primes[r]) {
            for (int i = r * r; i < limit; i += r * r)
                primes[i] = false;
        }
    }

    std::vector<int> primeNumbers;
    for (int a = 2; a <= limit; a++) {
        if (primes[a])
            primeNumbers.push_back(a);
    }

    return primeNumbers;
}

int main() {
    int limit = 20;
    std::vector<int> primeNumbers = sieve_of_atkin(limit);

    for (int prime : primeNumbers)
        std::cout << prime << " ";

    return 0;
}