Algorytmy i Struktury Danych - Sito Atkina Bernsteina

Problem

W przedziale <2,n> znaleźć wszystkie liczby pierwsze (ang. primes).

Sito Atkina-Bernsteina (ang. Atkin-Bernstein Sieve) należy do nowoczesnych algorytmów wyszukujących liczby pierwsze w przedziale od 2 do zadanej liczby naturalnej n. Jest to zoptymalizowany algorytm w stosunku do starożytnego sita Eratostenesa, który najpierw wykonuje wstępne przetwarzanie zbioru liczb, a następnie wyrzuca z niego wszystkie wielokrotności kwadratów początkowych liczb pierwszych zamiast samych liczb pierwszych. Twórcami algorytmu są dwaj profesorowie University of Illinois w Chicago: Arthur Oliver Lonsdale Atkin – emerytowany profesor matematyki oraz Daniel Julius Bernstein – profesor matematyki, kryptolog i programista.

Przedstawiony poniżej algorytm i jego realizacja nie są optymalne. Umieściłem je tutaj ze względów dydaktycznych. Prawdziwą implementację sita Atkina można znaleźć na stronie domowej prof. Bernsteina, czyli u źródeł i tam odsyłam wszystkich zainteresowanych tym tematem.

W przeciwieństwie do sita Eratostenesa, sito Atkina-Bernsteina rozpoczyna pracę ze zbiorem S, w którym wszystkie liczby są zaznaczone jako złożone (czyli nie pierwsze). Ma to sens, ponieważ algorytm ignoruje kompletnie liczby podzielne przez 2, 3 lub 5. Zaoszczędzamy czas na wyrzucaniu ich ze zbioru.

Algorytm opiera swoje działanie na następujących faktach matematycznych:

Wszystkie liczby dające resztę z dzielenia przez 60 równą 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56 lub 58 są podzielne przez 2, zatem nie są pierwsze – algorytm je ignoruje.
Wszystkie liczby dające resztę z dzielenia przez 60 równą 3, 9, 15, 21, 27, 33, 39, 45, 51 lub 57 są z kolei podzielne przez 3 i również nie są pierwsze – algorytm je ignoruje.
Wszystkie liczby dające resztę z dzielenia przez 60 równą 5, 25, 35 lub 55 są podzielne przez 5 i także nie są pierwsze – algorytm je ignoruje.
Wszystkie liczby dające resztę z dzielenia przez 60 równą 1, 13, 17, 29, 37, 41, 49 lub 53 posiadają resztę z dzielenia przez 12 równą 1 lub 5. Liczby te są pierwsze wtedy i tylko wtedy, gdy liczba rozwiązań równania 4x² + y² = n jest nieparzysta dla x,y N, a liczba n jest niepodzielna przez kwadraty liczb naturalnych.
Wszystkie liczby dające resztę z dzielenia przez 60 równą 7, 19, 31 lub 43 posiadają resztę z dzielenia przez 12 równą 7. Są one liczbami pierwszymi wtedy i tylko wtedy, gdy liczba rozwiązań równania 3x² + y² = n jest nieparzysta dla x,y N, a liczba n jest niepodzielna przez kwadraty liczb naturalnych.
Wszystkie liczby dające resztę z dzielenia przez 60 równą 11, 23, 47 lub 59 posiadają resztę z dzielenia przez 12 równą 11. Są one liczbami pierwszymi wtedy i tylko wtedy, gdy liczba rozwiązań równania 3x² − y² = n jest nieparzysta dla x,y N, a liczba n jest niepodzielna przez kwadraty liczb naturalnych.

Algorytm sita Atkina-Bernsteina

Wejście

–

liczba określająca górny kraniec przedziału poszukiwania liczb pierwszych, n

N, n > 4

Wyjście:

Kolejne liczby pierwsze w przedziale od 2 do n.

Zmienne pomocnicze

S	–	zbiór liczbowy. S ⊂ {5, 6, ..., n}
g	–	granica przetwarzania liczb w S[]. g N
x,y	–	używane w równaniach kwadratowych. x,y N
xx,yy	–	kwadraty x i y, xx,yy N
z	–	rozwiązanie równania, z N
i	–	indeks, i N

Lista kroków:

K01:	Dla i = 5,6,...,n wykonuj S[i] ← false	; inicjujemy zbiór liczbowy – wszystkie liczby zaznaczamy jako liczby złożone
K02:	g ← [√n]	; granica wyznaczania początkowych liczb pierwszych
K03:	Dla każdego (x,y) z <1,g> x <1,g> wykonuj kroki K04...K12	; szukamy rozwiązań równań kwadratowych
K04:	xx = x × x	; xx = x²
K05:	yy = y × y	; yy = y²
K06:	z ← (xx shl 2) + yy	; z = 4x² + y²
K07:	Jeśli (z ≤ n) ∧ (z mod 12 = 1 ∨ z mod 12 = 5), to S[z] ← ¬ S[z]	; jeśli spełnione są warunki, to zmieniamy na przeciwny ; stan liczby w S
K08:	z ← z - xx	; z = 3x² + y²
K09:	Jeśli (z ≤ n) ∧ (z mod 12 = 7), to S[z] ← ¬ S[z]
K10:	Jeśli x ≤ y, to następny obieg pętli K03	; ostatnie równanie sprawdzamy tylko dla x > y
K11:	z ← z - (yy shl 1)	; z = 3x² - y²
K12:	Jeśli (z ≤ n) ∧ ( z mod 12 = 11), to S[z] ← ¬ S[z]
K13:	Dla i = 5,6,...,g wykonuj kroki K14...K18	; eliminujemy liczby złożone sitem
K14:	Jeśli S[i] = false, to następny obieg pętli K13
K15:	xx = i × i; z ← xx	; z S eliminujemy wielokrotności kwadratów liczb pierwszych
K16:	Dopóki z ≤ n wykonuj kroki K17...K18
K17:	S[z] ← false
K18:	z ← z + xx
K19:	Pisz 2 3	; dwie pierwsze liczby wyprowadzamy bezwarunkowo
K20:	Dla i = 5,6,...,n wykonuj krok K21	; przeglądamy zbiór wypisując znalezione liczby pierwsze
K21:	Jeśli S[i], to pisz i
K22:	Zakończ

Program

Ważne:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program w pierwszym wierszu czyta liczbę n i w następnych wierszach wypisuje kolejne liczby pierwsze zawarte w przedziale od 2 do n.

Lazarus

// Sito Atkina-Bernsteina
// Data   : 21.03.2008
// (C)2012 mgr Jerzy Wałaszek
//----------------------------

program prg;

type
  Tbarray = array of boolean;

var n,g,x,y,xx,yy,z,i : longword;
    S : Tbarray;

begin
  readln(n);
  setlength(S,n+1);
  for i := 5 to n do S[i] := false;
  g := round(sqrt(n));
  for x := 1 to g do
  begin
    xx := x * x;
    for y := 1 to g do
    begin
      yy := y * y;
      z := (xx shl 2) + yy;
      if (z <= n) and ((z mod 12 = 1) or (z mod 12 = 5)) then S[z] := not S[z];
      dec(z,xx);
      if (z <= n) and (z mod 12 = 7) then S[z] := not S[z];
      if (x > y) then
      begin
        dec(z,yy shl 1);
        if (z <= n) and (z mod 12 = 11) then S[z] := not S[z];
      end; 
    end;
  end;
  for i := 5 to g do
    if S[i] then
    begin
      xx := i * i;
      z := xx;
      while z <= n do
      begin
        S[z] := false;
        inc(z,xx);
      end;
    end;
  write(2,' ',3,' ');
  for i := 5 to n do
    if S[i] then write(i,' ');
  writeln;
end.

Code::Blocks

// Sito Atkina-Bernsteina
// Data   : 21.03.2008
// (C)2012 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>
#include <cmath>

using namespace std;

int main()
{
  unsigned int n,g,x,y,xx,yy,z,i;
  bool * S;

  cin >> n;
  S = new bool[n+1];
  for(i = 5; i <= n; i++) S[i] = false;
  g = (unsigned int)(sqrt(n));
  for(x = 1; x <= g; x++)
  {
    xx = x * x;
    for(y = 1; y <= g; y++)
    {
      yy = y * y;
      z = (xx << 2) + yy;
      if((z <= n) && ((z % 12 == 1) || (z % 12 == 5))) S[z] = !S[z];
      z -= xx;
      if((z <= n) && (z % 12 == 7)) S[z] = !S[z];
      if(x > y)
      {
        z -= yy << 1;
        if((z <= n) && (z % 12 == 11)) S[z] = !S[z];
      } 
    }
  }
  for(i = 5; i <= g; i++)
    if(S[i])
    {
      xx = i * i;
      z = xx;
      while(z <= n)
      {
        S[z] = false;
        z += xx;
      }
    }
  cout << 2 << " " << 3 << " ";
  for(i = 5; i <= n; i++)
    if(S[i]) cout << i << " ";
  cout << endl;
  delete [] S;
  return 0;
}

Free Basic

' Sito Atkina-Bernsteina
' Data   : 21.03.2008
' (C)2012 mgr Jerzy Wałaszek
'----------------------------

Dim As Uinteger n,g,x,y,xx,yy,z,i

Input n

Dim As Byte S(n+1)

For i = 5 To n: S(i) = 0: Next
g = Cuint(Sqr(n))
For x = 1 To g
  xx = x * x
  For y = 1 To g
    yy = y * y
    z = (xx Shl 2) + yy
    If (z <= n) And ((z Mod 12 = 1) Or (z Mod 12 = 5)) Then S(z) = 1-S(z)
    z -= xx
    If (z <= n) And (z Mod 12 = 7) Then S(z) = 1-S(z)
    If (x > y) Then
      z -= yy Shl 1
      If (z <= n) And (z Mod 12 = 11) Then S(z) = 1-S(z)
    End If
  Next
Next
For i = 5 To g
  If S(i) = 1 Then
    xx = i * i
    z  = xx
    While z <= n
      S(z) = 0
      z += xx
    Wend
  End If
Next
Print 2;" ";3;" ";
For i = 5 To n
  If S(i) Then Print i;" ";
Next
Print
End

Wynik

100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97