Algorytmy i Struktury Danych - Sito Atkina Bernsteina

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 z 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:

n : 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.

Elementy 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.
sqr(x) : pierwiastek kwadratowy z x.

Lista kroków:

K01: Dla i = 5, 6, …, n: ; inicjujemy zbiór liczbowy – wszystkie
     wykonuj S[i] ← false ; liczby zaznaczamy jako liczby złożone
K02: g ← [sqr(n)] ; granica wyznaczania początkowych liczb pierwszych
K03: Dla każdego (x, y) z <1; g>: ; szukamy rozwiązań równań kwadratowych
     wykonuj kroki K04…K12
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), ; jeśli spełnione
       to S[z] ← ¬S[z] ; są warunki, 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, ; ostatnie równanie sprawdzamy tylko dla x > y
       to następny obieg pętli K03
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: ; eliminujemy sitem liczby złożone
     wykonuj kroki K14…K18
K14:   Jeśli S[i] = false, 
       to następny obieg pętli K13
K15:   xx = i×i ; z S eliminujemy wielokrotności
       z ← xx ; 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: ; przeglądamy zbiór wypisując znalezione
     wykonuj krok K21 ; liczby pierwsze
K21:   Jeśli S[i], 
       to pisz i
K22: Zakończ

Przykładowe programy

Uwaga:

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.

Pascal

// Sito Atkina-Bernsteina
// Data: 21.03.2008
// (C)2020 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.

C++

// Sito Atkina-Bernsteina
// Data: 21.03.2008
// (C)2020 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;
}

Basic

' Sito Atkina-Bernsteina
' Data: 21.03.2008
' (C)2020 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

Python (dodatek)

# Sito Atkina-Bernsteina
# Data: 11.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ---------------------------

import math

n = int(input())
s = [False]*(n+1)
g = int(math.sqrt(n))
for x in range(1, g+1):
    xx = x*x
    for y in range(1, g+1):
        yy = y*y
        z = (xx << 2)+yy
        if (z <= n) and \
           ((z%12 == 1) or (z%12 == 5)):
            s[z] = not s[z]
        z -= xx
        if (z <= n) and (z%12 == 7):
            s[z] = not s[z]
        if x > y:
            z -= yy << 1
            if (z <= n) and (z%12 == 11):
                s[z] = not s[z]
for i in range(5, g+1):
    if s[i]:
        xx = i*i
        z = xx
        while z <= n:
            s[z] = False
            z += xx
print(2, 3, end=" ")
for i in range(5, n+1):
    if s[i]:
        print(i, end=" ")
print()

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

JavaScript

<html>
  <head>
    <title>
      Sito Atkina-Bernsteina
    </title>
  </head>
  <body>

<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm"
        style="text-align: center;
               background-color:
               #E7E7DA">
          <b>
          Sito Atkina-Bernsteina
          </b><br/>
          &nbsp;&nbsp;(C)2026
          mgr Jerzy Wałaszek
          &nbsp;&nbsp;
          <hr>
          n = <input
              type="text"
              name="inp_n"
              size="16"
              value="100"
              style="text-align:
                     right">
          <hr>
          <input
          type="button"
          value="Wykonaj"
          name="B1"
          onclick="main()">
          <hr>
          <b>Wynik:</b>
          <div id="out">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script type=text/javascript
        language=javascript>

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

function main()
{
  var c,n,g,x,y,xx,yy,z,i,t,S;

  n = parseInt(document.frm
      .inp_n.value);
  c = 0;
  if(!isNaN(n))
  {
    S = new Array(n + 1);
    for(i = 5; i <= n; i++)
      S[i] = false;
    g = Math.floor(Math.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;
        }
      }
    t = "2 3 ";
    for(i = 5; i <= n; i++)
      if(S[i])
      {
        t += i+" ";
        c++;
        if(c == 10)
        {
          c = 0;
          t += "<br/>";
        }
      }
  }
  else t = "Złe dane wejściowe!";
  document.getElementById ("out")
  .innerHTML = t;
}

</script>

  </body>
</html>

do podrozdziału do strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone pod warunkiem podania źródła oraz niepobierania za to pieniędzy.
Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Liczby pierwsze – sito Atkina Bernsteina