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 – 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 \({4\cdot x^2+y^2=n}\) jest nieparzysta dla \({x, y\in\mathbb{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 \({3\cdot x^2+y^2=n}\) jest nieparzysta dla \({x,y\in\mathbb{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 \({3\cdot x^2−y^2=n}\) jest nieparzysta dla \({x,y\in\mathbb{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\in\mathbb{N},n>4}\).

Wyjście:

Kolejne liczby pierwsze w przedziale od \(2\) do \(n\).

Elementy pomocnicze:

\(S\) : zbiór liczbowy; \({S \in \{5, 6, …, n\}}\).
\(g\) : granica przetwarzania liczb w \(S\); \({g\in\mathbb{N}}\).
\(x,y\) : używane w równaniach kwadratowych; \({x,y\in\mathbb{N}}\).
\(xx, yy\) : kwadraty \(x\) i \(y\); \({xx,yy\in\mathbb{N}}\).
\(z\) : rozwiązanie równania; \({z\in\mathbb{N}}\).
\(i\) : indeks; \({i\in\mathbb{N}}\).

Lista kroków:

K01:	Dla \(i=5,6,\dots,n:\) ; inicjujemy zbiór liczbowy – wszystkie wykonuj \(S[i]\leftarrow\text{false}\) ; liczby zaznaczamy jako liczby złożone
K02:	\(g\leftarrow\left[\sqrt{n}\right]\) ; 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\leftarrow x\cdot x\) ; \(xx = x^2\)
K05:	\(yy\leftarrow y\cdot y\) ; \(yy = y^2\)
K06:	\(z\leftarrow (xx\ll 2)+yy\) ; \(z = 4\cdot x^2+y^2\)
K07:	Jeśli \((z\le n)\wedge((z\bmod12=1)\lor(z\bmod12 = 5)),\) ; jeśli spełnione to \(S[z]\leftarrow\neg S[z]\) ; są warunki, zmieniamy na przeciwny stan liczby w \(S\)
K08:	\(z\leftarrow z-xx\) ; \(z=3\cdot x^2+y^2\)
K09:	Jeśli \((z\le n)\wedge(z\bmod12=7),\) to \(S[z]\leftarrow\neg S[z]\)
K10:	Jeśli \(x\le y,\) ; ostatnie równanie sprawdzamy tylko dla \(x>y\) to następny obieg pętli K03
K11:	\(z\leftarrow z-(yy\ll1)\) ; \(z = 3\cdot x^2-y^2\)
K12:	Jeśli \((z\le n)\wedge(z\bmod12 = 11),\) to \(S[z]\leftarrow\neg S[z]\)
K13:	Dla \(i = 5,6,\dots,g:\) ; eliminujemy sitem liczby złożone wykonuj kroki K14…K18
K14:	Jeśli \(S[i]=\text{false},\) to następny obieg pętli K13
K15:	\(xx\leftarrow i\cdot i\) ; z \(S\) eliminujemy wielokrotności \(z\leftarrow xx\) ; kwadratów liczb pierwszych
K16:	Dopóki \(z\le n,\) wykonuj kroki K17…K18
K17:	\(S[z]\leftarrow\text{false}\)
K18:	\(z\leftarrow z+xx\)
K19:	Pisz \(2\;3\) ; dwie początkowe liczby pierwsze ; wyprowadzamy bezwarunkowo
K20:	Dla \(i=5,6,\dots,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; overflow-y: hidden;"
     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

Liczby pierwsze – sito Atkina Bernsteina