Сито объяснения Аткина

Question

Сейчас я занимаюсь проектом, и мне нужен эффективный метод для вычисления простых чисел. Я использовал сито Эратосфена , но я искал вокруг и обнаружил, что сито Аткина является более эффективным методом. Мне было трудно найти объяснение (которое я смог понять!) Этого метода. Как это работает? Пример кода (желательно на C или Python) будет блестящим.

Редактировать: спасибо за вашу помощь, единственное, что я до сих пор не понимаю, это то, на что ссылаются переменные x и y в псевдокоде. Может ли кто-нибудь пролить свет на это для меня?

Answer 1

Вики-страница - это всегда хорошее место для начала, поскольку она полностью объясняет алгоритм и предоставляет закомментированный псевдокод. (Внимание: здесь много подробностей, и поскольку сайт вики надежно работает, я не буду здесь его цитировать.)

Для ссылок на конкретные языки, которые вы упомянули:

• C реализация (оптимизировано)
• Реализация Python

Надеюсь, это поможет.

Answer 2

Статья в Википедии имеет объяснение:

• "Алгоритм полностью игнорирует любые числа, делимые на два, три или пять. Все числа с четным остатком по модулю шестьдесят делятся на два и не просты. Все числа с остатком по модулю шестьдесят делятся на три также делится на три и не простое число. Все числа, имеющие деление на шестьдесят по модулю, делимое на пять, делятся на пять и не простое. Все эти остатки игнорируются. "

Начнем со знаменитого

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 & mdash; ни один номер в 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-element test), просто предполагает , что его аргумент list является упорядоченным, возрастающим списком чисел, так что он может обнаружить как присутствие, так и присутствие:

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.

см. Также: Какой самый быстрый алгоритм поиска простых чисел?

---

Временная сложность сита Эратосфена при получении простых чисел до N равна O (N log log N) , а у сита Аткина - O (N) (откладывая в сторону дополнительные log log N оптимизации, которые плохо масштабируются). Тем не менее, общепринятая мудрость заключается в том, что постоянные факторы в сите Аткина намного выше, и поэтому он может окупиться только выше 32-битных чисел ( N> 2 ³² ) если вообще .

Answer 3

Редактировать: единственное, что я до сих пор не понимаю, это то, на что ссылаются переменные x и y в псевдокоде. Может ли кто-нибудь пролить свет на это для меня?

Для базового объяснения общего использования переменных 'x' и 'y' в псевдокоде страницы Wikipedia (или других лучших реализациях Sieve of Atkin) вы можете найти мой ответ другому связанный вопрос полезно.

Answer 4

Вот реализация сита Аткинса на 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;
 }

Answer 5

// 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;
}