Решето Аткина, объяснение
https://stackoverflow.com/questions/1023768
-
06-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
В данный момент я занимаюсь проектом, и мне нужен эффективный метод вычисления простых чисел.Я использовал решето Эратосфена но я искал вокруг и обнаружил, что сито Аткина это более эффективный метод.Мне было трудно найти объяснение (которое я смог понять!) этого метода.Как это работает?Пример кода (желательно на C или Python) был бы великолепен.
Редактировать: спасибо за вашу помощь, единственное, чего я до сих пор не понимаю, это то, на что ссылаются переменные x и y в псевдокоде.Может ли кто-нибудь пролить свет на это для меня?
Решение
А вики-страница всегда хорошее место для начала, поскольку оно полностью объясняет алгоритм и предоставляет псевдокод с комментариями.(Н.Б.Там много подробностей, и поскольку вики-сайт надежно работает, я не буду цитировать их здесь.)
Для ссылок на конкретных языках, которые вы упомянули:
- C-реализация (оптимизировано)
- Реализация Python
Надеюсь, это поможет.
Другие советы
Статья в Википедии есть объяснение:
- «Алгоритм полностью игнорирует любые числа, делящиеся на два, три или пять.Все числа с четным остатком по модулю шестьдесят делятся на два и не являются простыми.Все числа, у которых остаток по модулю шестьдесят делится на три, также делятся на три и не являются простыми.Все числа, у которых остаток по модулю шестьдесят делится на пять, делятся на пять и не являются простыми.Все эти остатки игнорируются».
Начнем со знаменитого
primes = sieve [2..] = sieve (2:[3..])
-- primes are sieve of list of 2,3,4... , i.e. 2 prepended to 3,4,5...
sieve (x:xs) = x : sieve [y | y <- xs, rem y x /= 0] -- set notation
-- sieve of list of (x prepended to xs) is x prepended to the sieve of
-- list of `y`s where y is drawn from xs and y % x /= 0
Давайте посмотрим, как это происходит на нескольких первых шагах:
primes = sieve [2..] = sieve (2:[3..])
= 2 : sieve p2 -- list starting w/ 2, the rest is (sieve p2)
p2 = [y | y <- [3..], rem y 2 /= 0] -- for y from 3 step 1: if y%2 /= 0: yield y
p2 не должно содержать кратных 2.IOW, он будет содержать только числа, взаимно простые с 2 — нет номера в p2 имеет 2 как его фактор.Найти p2 нам на самом деле не нужно проверять деление на 2 каждое число в [3..] (это 3 и вверх, 3,4,5,6,7,...), потому что мы можем перечислить все кратные 2 заранее:
rem y 2 /= 0 === not (ordElem y [2,4..]) -- "y is not one of 2,4,6,8,10,..."
ordElem как elem (т.е. is-элемент тест), это просто предполагает что его аргумент списка представляет собой упорядоченный, увеличивающийся список чисел, поэтому он может безопасно обнаруживать как отсутствие, так и присутствие:
ordElem y xs = take 1 (dropWhile (< y) xs) == [y] -- = elem y (takeWhile (<= y) xs)
Обычный elem ничего не предполагает, поэтому должен проверять каждый элемент своего аргумента списка, поэтому не может обрабатывать бесконечные списки. ordElem может.Итак, тогда,
p2 = [y | y <- [3..], not (ordElem y [2,4..])] -- abstract this as a function, diff a b =
= diff [3..] [2,4..] -- = [y | y <- a, not (ordElem y b)]
-- 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
-- . 4 . 6 . 8 . 10 . 12 . 14 . 16 . 18 . 20 . 22 .
= diff [3..] (map (2*) [2..] ) -- y > 2, so [4,6..] is enough
= diff [2*k+x | k <- [0..], x <- [3,4]] -- "for k from 0 step 1: for x in [3,4]:
[2*k+x | k <- [0..], x <- [ 4]] -- yield (2*k+x)"
= [ 2*k+x | k <- [0..], x <- [3 ]] -- 2 = 1*2 = 2*1
= [3,5..] -- that's 3,5,7,9,11,...
p2 это просто список всех шансов выше 2.Ну, мы знали что.А что насчет следующего шага?
sieve p2 = sieve [3,5..] = sieve (3:[5,7..])
= 3 : sieve p3
p3 = [y | y <- [5,7..], rem y 3 /= 0]
= [y | y <- [5,7..], not (ordElem y [3,6..])] -- 3,6,9,12,...
= diff [5,7..] [6,9..] -- but, we've already removed the multiples of 2, (!)
-- 5 . 7 . 9 . 11 . 13 . 15 . 17 . 19 . 21 . 23 . 25 . 27 .
-- . 6 . . 9 . . 12 . . 15 . . 18 . . 21 . . 24 . . 27 .
= diff [5,7..] (map (3*) [3,5..]) -- so, [9,15..] is enough
= diff [2*k+x | k <- [0..], x <- [5]] (map (3*)
[2*k+x | k <- [0..], x <- [ 3]] )
= diff [6*k+x | k <- [0..], x <- [5,7,9]] -- 6 = 2*3 = 3*2
[6*k+x | k <- [0..], x <- [ 9]]
= [ 6*k+x | k <- [0..], x <- [5,7 ]] -- 5,7,11,13,17,19,...
И следующий,
sieve p3 = sieve (5 : [6*k+x | k <- [0..], x <- [7,11]])
= 5 : sieve p5
p5 = [y | y <- [6*k+x | k <- [0..], x <- [7,11]], rem y 5 /= 0]
= diff [ 6*k+x | k <- [0..], x <- [7,11]] (map (5*)
[ 6*k+x | k <- [0..], x <- [ 5, 7]]) -- no mults of 2 or 3!
= diff [30*k+x | k <- [0..], x <- [7,11,13,17,19,23,25,29,31,35]] -- 30 = 6*5 = 5*6
[30*k+x | k <- [0..], x <- [ 25, 35]]
= [ 30*k+x | k <- [0..], x <- [7,11,13,17,19,23, 29,31 ]]
Именно в такой последовательности и работает решето Аткина.Нет кратных 2, 3 или 5 присутствуют в нем.Аткин и Бернштейн работают по модулю 60, т.е.они удваивают диапазон:
p60 = [ 60*k+x | k <- [0..], x <- [1, 7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59]]
Затем они используют некоторые сложные теоремы, чтобы узнать некоторые свойства каждого из этих чисел и обращаться с каждым соответствующим образом.Все повторяется (а-ля «колесо») с периодом 60.
- «Все числа n с остатком по модулю шестьдесят 1, 13, 17, 29, 37, 41, 49 или 53 (...) являются простыми тогда и только тогда, когда число решений задачи
4x^2 + y^2 = nнечетно и число не содержит квадратов».
Что это значит?Если мы знаем, что число решений этого уравнения нечетно для такого числа, то оно простое. если это без квадратов.Это означает, что в нем нет повторяющихся факторов (например, 49, 121, и т. д).
Аткин и Бернштейн используют это, чтобы сократить общее количество операций:для каждого простого числа (от 7 и выше) мы перечисляем (и отмечаем для удаления) кратные его площадь, поэтому они находятся гораздо дальше друг от друга, чем в решете Эратосфена, поэтому в данном диапазоне их меньше.
Остальные правила такие:
«Все числа n с остатком по модулю шестьдесят 7, 19, 31 или 43 (...) являются простыми тогда и только тогда, когда число решений задачи
3x^2 + y^2 = nнечетно и число не содержит квадратов».«Все числа n с остатком по модулю шестьдесят 11, 23, 47 или 59 (...) являются простыми тогда и только тогда, когда число решений задачи
3x^2 − y^2 = nнечетно и число не содержит квадратов».
Это учитывает все 16 основных чисел в определении p60.
смотрите также: Какой алгоритм поиска простых чисел является самым быстрым?
Временная сложность решета Эратосфена при получении простых чисел до Н является O (N журнал журнал N), а у решета Аткина НА) (без учета дополнительных log log N оптимизации, которые плохо масштабируются).Однако общепринятая мудрость заключается в том, что постоянные коэффициенты в решете Аткина намного выше, и поэтому оно может окупиться только выше 32-битных чисел (Н > 232), если вообще.
Редактировать: единственное, что я до сих пор не понимаю, это то, на что ссылаются переменные x и y в псевдокоде. Может ли кто-нибудь пролить свет на это для меня?
Для базового объяснения общего использования переменных 'x' и 'y' в псевдокоде страницы Wikipedia (или других лучших реализациях Sieve of Atkin) вы можете найти мой ответ на другой связанный вопрос полезно.
Вот реализация сита atkins на c ++, которая может вам помочь ...
#include <iostream>
#include <cmath>
#include <fstream>
#include<stdio.h>
#include<conio.h>
using namespace std;
#define limit 1000000
int root = (int)ceil(sqrt(limit));
bool sieve[limit];
int primes[(limit/2)+1];
int main (int argc, char* argv[])
{
//Create the various different variables required
FILE *fp=fopen("primes.txt","w");
int insert = 2;
primes[0] = 2;
primes[1] = 3;
for (int z = 0; z < limit; z++) sieve[z] = false; //Not all compilers have false as the default boolean value
for (int x = 1; x <= root; x++)
{
for (int y = 1; y <= root; y++)
{
//Main part of Sieve of Atkin
int n = (4*x*x)+(y*y);
if (n <= limit && (n % 12 == 1 || n % 12 == 5)) sieve[n] ^= true;
n = (3*x*x)+(y*y);
if (n <= limit && n % 12 == 7) sieve[n] ^= true;
n = (3*x*x)-(y*y);
if (x > y && n <= limit && n % 12 == 11) sieve[n] ^= true;
}
}
//Mark all multiples of squares as non-prime
for (int r = 5; r <= root; r++) if (sieve[r]) for (int i = r*r; i < limit; i += r*r) sieve[i] = false;
//Add into prime array
for (int a = 5; a < limit; a++)
{
if (sieve[a])
{
primes[insert] = a;
insert++;
}
}
//The following code just writes the array to a file
for(int i=0;i<1000;i++){
fprintf(fp,"%d",primes[i]);
fprintf(fp,", ");
}
return 0;
}
// Title : Seive of Atkin ( Prime number Generator)
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
int main()
{
ios_base::sync_with_stdio(false);
long long int n;
cout<<"Enter the value of n : ";
cin>>n;
vector<bool> is_prime(n+1);
for(long long int i = 5; i <= n; i++)
{
is_prime[i] = false;
}
long long int lim = ceil(sqrt(n));
for(long long int x = 1; x <= lim; x++)
{
for(long long int y = 1; y <= lim; y++)
{
long long int num = (4*x*x+y*y);
if(num <= n && (num % 12 == 1 || num%12 == 5))
{
is_prime[num] = true;
}
num = (3*x*x + y*y);
if(num <= n && (num % 12 == 7))
{
is_prime[num] = true;
}
if(x > y)
{
num = (3*x*x - y*y);
if(num <= n && (num % 12 == 11))
{
is_prime[num] = true;
}
}
}
}
// Eliminating the composite by seiveing
for(long long int i = 5; i <= lim; i++)
{
if(is_prime[i])
for(long long int j = i*i; j <= n; j += i)
{
is_prime[j] = false;
}
}
// Here we will start printing of prime number
if(n > 2)
{
cout<<"2\t"<<"3\t";
}
for(long long int i = 5; i <= n; i++)
{
if(is_prime[i])
{
cout<<i<<"\t";
}
}
cout<<"\n";
return 0;
}
