Решето́ А́ткина — метод нахождения всех простых чисел, не превышающих заданное натуральное limit. В отличие от решета Эратосфена, последовательно исключающего числа, кратные уже найденным простым, данный алгоритм производит предварительное просеивание, а затем из найденных чисел исключает кратные квадратам простых, благодаря чему имеет теоретически лучший показатель асимптотической сложности $\mathop {O} ({\frac {1}{\log \log(limit)}})$ .

Авторы: Артур Оливер Лонсдэйл Аткин, Дэниэл Юлиус Бернштайн (англ. Arthur Oliver Lonsdale Atkin, Daniel Julius Bernstein).

Алгоритм

Числа 2, 3 и 5 рассматриваются, как заведомо простые.
Создаётся «решето» — список, сопоставляющий каждому натуральному числу n из диапазона (5; limit] флаг «простое\составное». Изначально все флаги устанавливаются в положение «составное».
Для очередного n из диапазона (5; limit] находится остаток от деления на 60:
- 1, 13, 17, 29, 37, 41, 49, 53: уравнение — 4x² + y² = n; (1)
- 7, 19, 31, 43: уравнение — 3x² + y² = n; (2)
- 11, 23, 47, 59: уравнение — 3x² − y² = n (при x > y) (3)
Значение остатка, не равное ни одному из вышеуказанных, свидетельствует о том, что число n — составное и таким образом исключается из дальнейшего рассмотрения. Иначе, если соответствующее уравнение имеет нечётное число решений (пар натуральных значений x и y), n помечается в «решете» как «простое».
Процедура повторяется начиная с шага № 3, пока не будут перебраны все n ≤ limit.
В заключение для каждого найденного простого n все значения, кратные его квадрату, помечаются в «решете» как «составные».

Реализации

В приведённых ниже реализациях предприняты следующие шаги, направленные на оптимизацию вычислений:

Флаги «простое\составное» хранятся упакованными в массив битов (vector<bool> для С++, BitArray для C#, BitSet для Java, BitVec для Rust и т. п.), возвращаемый функцией. Значение n-го бита true соответствует простому n, false — составному.
Умножение на степень двойки производится операцией побитового сдвига (<<). Например a << 2 вместо 4 * a.
Во избежание переполнения вычисления ведутся в длинном целом типе, в «решето» значения заносятся в обычном целом типе (unsigned long long / unsigned в C++, long / int в C# и Java, u64 / u32 в Rust).
Квадраты последовательных натуральных чисел вычисляются по формуле (x + 1)² = x² + 2x + 1, т. е. прибавлением последовательных нечётных чисел. Другими словами, если squares — последовательность квадратов натуральных чисел {0² = 0, 1² = 1, 2² = 4…}, а odds — последовательность нечётных чисел {1, 3, 5…}, то square_i+1 = square_i + odds_i (для любого i ≥ 0).
n делится с остатком на 12 (а не на 60). Полученные значения 1 и 5 соответствуют случаю (1), 7 — (2), 11 — (3), прочие — составному n.

C++

// atkin_sieve.h
#include <vector>


namespace AtkinSieve {
    std::vector<bool> getPrimesUpTo (unsigned const limit);
}


// atkin_sieve.cpp
#include "atkin_sieve.h"


std::vector<bool> AtkinSieve::getPrimesUpTo (unsigned const limit) {
    // Предварительное просеивание
    std::vector<bool> sieve(limit);    
    for (unsigned long long x2 = 1ull, dx2 = 3ull; x2 < limit; x2 += dx2, dx2 += 2ull)
        for (unsigned long long y2 = 1ull, dy2 = 3ull, n; y2 < limit; y2 += dy2, dy2 += 2ull) {
            // n = 4x² + y²
            n = (x2 << 2ull) + y2;
            if (n <= limit && (n % 12ull == 1ull || n % 12ull == 5ull))
                sieve[n].flip();
            // n = 3x² + y²
            n -= x2;
            if (n <= limit && n % 12ull == 7ull)
                sieve[n].flip();
            // n = 3x² - y² (при x > y)
            if (x2 > y2) {
                n -= y2 << 1ull;
                if (n <= limit && n % 12ull == 11ull)
                    sieve[n].flip();
            }
        }
    // Все числа, кратные квадратам, помечаются как составные
    unsigned r = 5u;
    for (unsigned long long r2 = r * r, dr2 = (r << 1ull) + 1ull; r2 < limit; ++r, r2 += dr2, dr2 += 2ull)
        if (sieve[r])
            for (unsigned long long mr2 = r2; mr2 < limit; mr2 += r2)
                sieve[mr2] = false;
    // Числа 2 и 3 — заведомо простые
    if (limit > 2u)
        sieve[2u] = true;
    if (limit > 3u)
        sieve[3u] = true;
    return sieve;
}

Пример использования:

// atkin_sieve_test.cpp
#include <iostream>

#include "atkin_sieve.h"


int main () {
    unsigned limit;
    std::cin >> limit;
    std::vector<bool> primes = AtkinSieve::getPrimesUpTo(limit);
    // Вывод последовательности
    for (unsigned number = 2u; number < limit; ++number)
        if (primes[number])
            std::cout << number << ' ';
    std::cout << '\n';
    return 0;
}

C#

// AtkinSieve.cs
using System;
using System.Collections;


public static class AtkinSieve
{
    public static BitArray GetPrimesUpTo (int limit)
    {
        var sieve = new BitArray(limit + 1);
        // Предварительное просеивание
        for (long x2 = 1L, dx2 = 3L; x2 < limit; x2 += dx2, dx2 += 2L)
            for (long y2 = 1L, dy2 = 3L, n; y2 < limit; y2 += dy2, dy2 += 2L)
            {
                // n = 4x² + y²
                n = (x2 << 2) + y2;
                if (n <= limit && (n % 12L == 1L || n % 12L == 5L))
                    sieve[(int)n] ^= true;
                // n = 3x² + y²
                n -= x2;
                if (n <= limit && n % 12L == 7L)
                    sieve[(int)n] ^= true;
                // n = 3x² - y² (при x > y)
                if (x2 > y2)
                {
                    n -= y2 << 1;
                    if (n <= limit && n % 12L == 11L)
                        sieve[(int)n] ^= true;
                }
            }
        // Все числа, кратные квадратам, помечаются как составные
        int r = 5;
        for (long r2 = r * r, dr2 = (r << 1) + 1L; r2 < limit; ++r, r2 += dr2, dr2 += 2L)
            if (sieve[r])
                for (long mr2 = r2; mr2 < limit; mr2 += r2)
                    sieve[(int)mr2] = false;
        // Числа 2 и 3 — заведомо простые
        if (limit > 2)
            sieve[2] = true;
        if (limit > 3)
            sieve[3] = true;
        return sieve;
    }
}

Пример использования:

// AtkinSieveTest.cs
using System;
using System.Collections;
using System.IO;


public static class AtkinSieveTest
{ 
    public static void Main ()
    {
        int limit = int.Parse(Console.ReadLine());
        var primes = AtkinSieve.GetPrimesUpTo(limit);
        // Вывод последовательности
        for (int number = 2; number < limit; ++number)
            if (primes[number])
                Console.Write("{0} ", number);
        Console.WriteLine();
    }
}

Java

// AtkinSieve.java
import java.util.BitSet;


public class AtkinSieve {
    public static BitSet getPrimesUpTo (int limit) {
        BitSet sieve = new BitSet();
        // Предварительное просеивание
        for (long x2 = 1L, dx2 = 3L; x2 < limit; x2 += dx2, dx2 += 2L)
            for (long y2 = 1L, dy2 = 3L, n; y2 < limit; y2 += dy2, dy2 += 2L) {
                // n = 4x² + y²
                n = (x2 << 2L) + y2;
                if (n <= limit && (n % 12L == 1L || n % 12L == 5L))
                    sieve.flip((int)n);
                // n = 3x² + y²
                n -= x2;
                if (n <= limit && n % 12L == 7L)
                    sieve.flip((int)n);
                // n = 3x² - y² (при x > y)
                if (x2 > y2) {
                    n -= y2 << 1L;
                    if (n <= limit && n % 12L == 11L)
                        sieve.flip((int)n);
                }
            }
        // Все числа, кратные квадратам, помечаются как составные
        int r = 5;
        for (long r2 = r * r, dr2 = (r << 1L) + 1L; r2 < limit; ++r, r2 += dr2, dr2 += 2L)
            if (sieve.get(r))
                for (long mr2 = r2; mr2 < limit; mr2 += r2)
                    sieve.set((int)mr2, false);
        // Числа 2 и 3 — заведомо простые
        if (limit > 2)
            sieve.set(2, true);
        if (limit > 3)
            sieve.set(3, true);
        return sieve;
    }
}

Пример использования:

// AtkinSieveTest.java
import java.util.BitSet;
import java.util.Scanner;


public class AtkinSieveTest {
    public static void main (String[] args) {
        Scanner scanner = new Scanner(System.in);
        int limit = scanner.nextInt();
        BitSet primes = AtkinSieve.getPrimesUpTo(limit);
        // Вывод последовательности
        for (int number = 2; number < limit; ++number)
            if (primes.get(number))
                System.out.format("%d ", number);
        System.out.println();
    }
}

Rust

Поскольку BitVec не входит в стандартную библиотеку Rust, его следует подключить через crates.io, указав bit-vec = "0.4.3" в разделе [dependencies] файла Cargo.toml.

// atkin_sieve.rs
extern crate bit_vec;

use self::bit_vec::BitVec;


pub fn get_primes_up_to (limit: u32) -> BitVec {
    let mut sieve = BitVec::from_elem(limit as usize, false);
    {
        let mut x2 = 1u64;
        let mut dx2 = 3u64;
        while x2 < limit as u64 {
            let mut y2 = 1u64;
            let mut dy2 = 3u64;
            let mut n: u64;
            let mut item: bool;
            while y2 < limit as u64 {
                // Просеивание:
                // n = 4x² + y²
                n = (x2 << 2) + y2;
                if n <= limit as u64 && (n % 12 == 1 || n % 12 == 5) {
                    item = sieve[n as usize];
                    sieve.set(n as usize, !item);
                }
                // n = 3x² + y²
                n -= x2;
                if n <= limit as u64 && n % 12 == 7 {
                    item = sieve[n as usize];
                    sieve.set(n as usize, !item);
                }
                // n = 3x² - y² (при x > y)
                if x2 > y2 {
                    n -= y2 << 1;
                    if n <= limit as u64 && n % 12 == 11 {
                        item = sieve[n as usize];
                        sieve.set(n as usize, !item);
                    }
                }
                y2 += dy2;
                dy2 += 2;
            }
            x2 += dx2;
            dx2 += 2;
        }
    }
    // Все числа, кратные квадратам, помечаются как составные
    {
        let mut r = 5u32;
        let mut r2 = (r * r) as u64;
        let mut dr2 = (r << 1) as u64 + 1;
        while r2 < limit as u64 {
            if sieve[r as usize] {
                let mut mr2 = r2;
                while mr2 < limit as u64 {
                    sieve.set(mr2 as usize, false);
                    mr2 += r2;
                }
            }
            r += 1;
            r2 += dr2;
            dr2 += 2;
        }
    }
    // 2 и 3 — заведомо простые
    if limit > 2 {
        sieve.set(2, true);
    }
    if limit > 3 {
        sieve.set(3, true);
    }
    sieve
}

Пример использования:

// atkin_sieve_test.rs
mod atkin_sieve;


// Функция для ввода данных
fn read_var <Type: std::str::FromStr> (var: &mut Type) -> bool {
    let mut input_text = String::new();
    std::io::stdin()
        .read_line(&mut input_text)
        .expect("failed to read from stdin");
    match input_text.trim().parse::<Type>() {
        Ok(value) => { *var = value; true },
        Err(..) => { false }
    }
}

fn main () {
    let mut limit = 0u32;
    if read_var(&mut limit) {
        let primes = atkin_sieve::get_primes_up_to(limit);
        for i in 2usize..limit as usize {
            if primes[i] {
                print!("{} ", i);
            }
        }
    }
}

В Cargo.toml следует добавить:

[[bin]]
name = "atkin_sieve_test"
path = "src/atkin_sieve_test.rs"

Go

package atkinsieve

import "math"

func AtkinSieve(limit uint) []bool {
	// Инициализация чиназеса
	sqr_lim := uint(math.Sqrt(float64(limit)))
	is_prime := make([]bool, limit+1)
	is_prime[2] = true
	is_prime[3] = true

	// Предположительно простые — это целые с нечётным числом
	// представлений в данных квадратных формах.
	// x2 и y2 — это квадраты i и j (оптимизация).
	var n uint
	x2 := uint(0)
	for i := uint(1); i <= sqr_lim; i++ {
		x2 += 2*i - 1
		y2 := uint(0)
		for j := uint(1); j <= sqr_lim; j++ {
			y2 += 2*j - 1
			n = 4*x2 + y2
			if (n <= limit) && (n%12 == 1 || n%12 == 5) {
				is_prime[n] = !is_prime[n]
			}
			// n = 3 * x2 + y2;
			n -= x2 // Оптимизация
			if (n <= limit) && (n%12 == 7) {
				is_prime[n] = !is_prime[n]
			}
			// n = 3 * x2 - y2;
			n -= 2 * y2 // Оптимизация
			if (i > j) && (n <= limit) && (n%12 == 11) {
				is_prime[n] = !is_prime[n]
			}

		}
	}
	// Отсеиваем кратные квадратам простых чисел в интервале [5, sqrt(limit)].
	// (основной этап не может их отсеять)
	for i := uint(5); i <= sqr_lim; i++ {
		if is_prime[i] {
			n = i * i
			for j := n; j <= limit; j += n {
				is_prime[j] = false
			}
		}
	}
	return is_prime
}

См. также

Решето Сундарама

Решето Эратосфена

Ссылки

Взятая за основу реализация алгоритма на C++