Решето́ А́ткина — метод нахождения всех простых чисел, не превышающих заданное натуральное 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

	// Предположительно простые — это целые с нечётным числом ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
	// представлений в данных квадратных формах. 卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐卐
	// 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
}