Liczby pierwsze można wyszukiwać poprzez eliminację ze zbioru liczbowego wszystkich wielokrotności wcześniejszych liczb.

Przykład:

Mamy następujący zbiór liczb naturalnych:

Ze zbioru wyrzucamy wszystkie wielokrotności pierwszej liczby 2. Otrzymujemy następujący zbiór:

W zbiorze pozostały liczby nieparzyste – z wyjątkiem pierwszej liczby 2. Liczby podzielne przez dwa zostały wyeliminowane. Teraz eliminujemy wielokrotności kolejnej liczby 3. Otrzymamy następujący zbiór:

Teraz w zbiorze pozostały liczby niepodzielne przez 2 i 3 – z wyjątkiem pierwszych 2 i 3. Zwróć uwagę, iż kolejna liczba 4 została już ze zbioru wyeliminowana. Skoro tak, to ze zbioru zniknęły również wszystkie wielokrotności liczby 4, ponieważ są one również wielokrotnościami liczby 2, a te wyeliminowaliśmy przecież na samym początku. Przechodzimy zatem do liczby 5 i eliminujemy jej wielokrotności otrzymując zbiór wynikowy:

Oprócz 2, 3 i 5 pozostałe w zbiorze liczby nie dzielą się już przez 2, 3 i 5. Liczba 6 jest wyeliminowana (wielokrotność 2 ), zatem przechodzimy do 7. Po wyeliminowaniu wielokrotności liczby 7 zbiór przyjmuje postać:

W zbiorze pozostały same liczby pierwsze.

Przy eliminacji wystarczy usunąć ze zbioru wielokrotności liczb leżących w przedziale od 2 do pierwiastka z n. Wyjaśnienie tego faktu jest identyczne jak w algorytmie szukania liczb pierwszych przez testowanie podzielności. Jeśli liczba p  jest złożona, to musi posiadać czynniki pierwsze w przedziale od 2 do pierwiastka z p.

Powyższe operacje wyrzucania wielokrotności prowadzą do przesiania zbioru wejściowego. W zbiorze pozostają tylko liczby pierwsze, liczby będące wielokrotnościami poprzedzających je liczb zostają ze zbioru odsiane. Algorytm dokonujący tej eliminacji nosi nazwę sita Eratostenesa (ang. Eratosthenes' sieve )  i jest jednym z najszybszych sposobów generowania liczb pierwszych.

Sito Eratostenesa jest algorytmem dwuprzebiegowym. Najpierw dokonuje on eliminacji ze zbioru liczb złożonych oznaczając je w określony sposób, a w drugim obiegu przegląda zbiór ponownie, wyprowadzając na wyjście liczby, które nie zostały oznaczone. Podstawowe znaczenie w tym algorytmie ma wybór odpowiedniej struktury danych do reprezentacji  zbioru liczb. Na tym polu można dokonywać różnych optymalizacji. W pierwszym podejściu zastosujemy tablicę wartości logicznych S. Element S [ i  ] będzie odpowiadał liczbie o wartości i. Zawartość S [ i  ] będzie z kolei informowała o tym, czy liczba i  pozostała w zbiorze (S [ i  ] = true ) lub została usunięta (S [ i  ] = false ). 

Najpierw przygotowujemy tablicę reprezentującą zbiór liczbowy wypełniając ją wartościami logicznymi true. Odpowiada to umieszczeniu w zbiorze wszystkich liczb wchodzących w zakres od 2 do n. Następnie z tablicy będziemy usuwali kolejne wielokrotności początkowych liczb od 2 do pierwiastka całkowitego z n w pisując do odpowiednich elementów wartość logiczną false. Na koniec przeglądniemy zbiór i wypiszemy indeksy elementów zawierających wartość logiczną true – odpowiadają one liczbom, które w zbiorze pozostały.

Za pierwszą wielokrotność do wyrzucenia ze zbioru przyjmiemy kwadrat liczby i. Przyjrzyj się naszemu przykładowi. Gdy wyrzucamy wielokrotności liczby 2, to pierwszą z nich jest 4 = 2 ². Następnie dla wielokrotności liczby 3 pierwszą do wyrzucenia jest 9 = 3 ², gdyż 6 zostało wyrzucone wcześniej jako wielokrotność 2. Dla 5 będzie to 25 = 5 ², gdyż 10 i 20 to wielokrotności 2, a 15 jest wielokrotnością 3, itd. Pozwoli to wyeliminować zbędne obiegi pętli usuwającej wielokrotności.

Algorytm sita Eratostenesa

// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

program prg;

type
  Tbarray = array of boolean;

var i, n, w : longword;
    S   : Tbarray;

begin
  readln ( n );
  setlength ( S, n + 1 );
  for i := 2 to n do S [ i ] := true;
  for i := 2 to round ( sqrt ( n ) ) do
    if S [ i ] then
    begin
      w := i * i;
      while w <= n do
      begin
        S [ w ] := false; inc ( w, i );
      end;
    end;
  for i := 2 to n do if S [ i ] then write ( i, ' ' );
  writeln;
end.

// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>
#include <cmath>

using namespace std;

int main( )
{
  unsigned int g, i, n, w;
  bool * S;

  cin >> n;
  S = new bool [ n + 1 ];
  for( i = 2; i <= n; i++ ) S [ i ] = true;
  g = ( unsigned int )sqrt ( n );
  for( i = 2; i <= g; i++ )
    if( S [ i ] )
    {
     w = i * i;
      while( w <= n )
      {
        S [ w ] = false; w += i;
      }
    } 
  for( i = 2; i <= n; i++ ) if( S [ i ] ) cout << i << " ";
  cout << endl;
  delete [ ] S;
  return 0;
}

' Sito Eratostenesa
' Data   : 17.03.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

Dim As Uinteger g, i, n, w

Input n
  
Dim As Byte S ( n += 1 )

For i = 2 To n: S ( i ) = 1: Next
g = Cuint ( Sqr ( n ) )
For i = 2 To g
  If S ( i ) > 0 Then
    w = i * i
    While w <= n
      S ( w ) = 0: w += i
    Wend
  End If
Next
For i = 2 To n
  If S ( i ) > 0 Then Print i;" ";
Next
Print
End

Jedyną parzystą liczbą pierwszą jest 2. Zatem wszystkie pozostałe liczby parzyste ( 4, 6, 8, ... ) już nie mogą być pierwsze. Utwórzmy tablicę S, której elementy będą reprezentowały kolejne liczby nieparzyste, począwszy od 3. Poniżej przedstawiamy początkowe przypisania indeksów liczbom nieparzystym.

Komórka o indeksie i  oznacza liczbę o wartości 2i + 1. Np. komórka 7 reprezentuje liczbę 2 × 7 + 1 = 15.

Liczba o wartości k  jest reprezentowana przez komórkę ( k - 1 ) / 2. Np. liczbę 45 reprezentuje komórka ( 45 - 1 ) / 2 = 22.

Z poprzedniego algorytmu pamiętamy, że pierwszą wyrzucaną wielokrotnością jest zawsze kwadrat liczby podstawowej. Nasza tablica rozpoczyna się od liczby 3, zatem pierwszą wielokrotnością, którą usuniemy, będzie liczba 9 na pozycji 4:

Zauważ, że kolejne wielokrotności liczby 3, które są reprezentowane w tablicy, znajdują się w niej w odstępach co 3:

Oznaczmy przez p  odstępy pomiędzy wielokrotnościami, a przez q  pozycję kwadratu liczby, której wielokrotności usuwamy z tablicy. Na początku mamy:

Po usunięciu tych wielokrotności tablica S  nie będzie zawierała dalszych liczb podzielnych przez 3. Przejdźmy do kolejnej liczby, czyli do 5. Kwadrat 5 znajduje się na pozycji 12:

Kolejne wielokrotności 5 znajdują się w naszej tablicy w odstępach co 5 (niektóre z nich są usunięte, ponieważ są również wielokrotnościami liczby 3) :

Zapiszmy:

Następna liczba to 7 i jej pierwsza wielokrotność 49 na pozycji 24. Kolejne wielokrotności są w odstępach co 7:

Znów zapiszmy:

Następna liczba to 9 i jej pierwsza wielokrotność 81 na pozycji 40. Kolejne wielokrotności są w odstępach co 9:

Znów zapiszmy:

Wyrzucanie wielokrotności liczby 9 możemy pominąć, ponieważ sama liczba 9 jest już wyrzucona z tablicy przez 3. Ostatnią rozważaną tu liczbą jest 11 o kwadracie równym 121. Liczba 121 znajduje się na pozycji ( 121-1 ) / 2 = 60:

W tablicy ta pozycja już nie występuje (indeksy kończą się na wartości 40 ), zatem kończymy. W S  pozostały jedynie liczby pierwsze. Na koniec wyprowadźmy wzory rekurencyjne dla kolejnych wartości p  oraz q. W tym celu wystarczy zauważyć prostą zależność (dla dociekliwych - wynika ona z rozkładu kwadratów liczb nieparzystych na sumę różnic - patrz: całkowity pierwiastek kwadratowy ):

Zatem możemy zapisać wzór rekurencyjny:

Po tych rozważaniach przystępujemy do zapisu algorytmu.

Algorytm ulepszonego sita Eratostenesa

// Program: 0014
// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

program prg;

type
  Tbarray = array of boolean;

var i, k, p, q, n, m : longword;
    S           : Tbarray;

begin
  readln ( n );
  if n and 1 = 1 then inc ( n );
  m  := n shr 1;
  setlength ( S, m+1 );
  for i := 1 to m - 1 do S [ i ] := true;
  i := 1; p := 3; q := 4;
  repeat
    if S [ i ] then
    begin
      k := q;
      while k < m do
      begin
        S [ k ] := false; inc ( k, p );
      end;
    end;
    inc ( i ); inc ( p, 2 );
    inc ( q, ( p shl 1 ) - 2 );
  until q >= m;
  write ( 2, ' ' );
  for i := 1 to m - 1 do if S [ i ] then write ( ( i  shl 1 ) + 1, ' ' );
  writeln;
end.

// Program: 0014
// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>

using namespace std;

int main( )
{
  unsigned int i, k, p, q, n, m;
  bool * S;

  cin >> n;
  if( n & 1 ) n++;
  m = n >> 1;
  S = new bool [ m + 1 ];
  for( i = 1; i < m; i++ ) S [ i ] = true;
  i = 1; p = 3; q = 4;
  do
  {
    if( S [ i ] )
    {
      k = q;
      while( k < m )
      {
        S [ k ] = false; k += p;
      }
    }
    i++; p += 2;
    q += ( p << 1 ) - 2;
  } while( q < m );
  cout << 2 << " ";
  for( i = 1; i < m; i++ )
    if( S [ i ] ) cout << ( i << 1 ) + 1 << " ";
  cout << endl;
  delete [ ] S;
  return 0;
}

' Program: 0014
' Sito Eratostenesa
' Data   : 17.03.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

Dim As Uinteger i, k, p, q, n, m

Input n
If n And 1 = 1 Then n += 1
m = n Shr 1

Dim As Byte S ( m + 1 )

For i = 1 To m - 1:
  S ( i ) = 1:
Next
i = 1: p = 3: q = 4
Do
  If S ( i ) = 1 Then
    k = q
    While k < m
      S ( k ) = 0: k += p
    Wend
  End If
  i += 1: p += 2:
  q += ( p Shl 1 ) - 2
Loop Until q >= m
Print 2;" ";
For i = 1 To m - 1
  If S ( i ) = 1 Then
    Print ( i Shl 1 ) + 1;" ";
  End If
Next
Print
End

Autorem prezentowanego algorytmu jest chiński informatyk Xuedong Luo. W rozwiązaniu 2 zbiór liczbowy ograniczyliśmy do liczb nieparzystych. Teraz pójdziemy dalej i wrzucimy z tego zbioru dodatkowo wszystkie liczby podzielne przez 3. W efekcie nasz zbiór S przybierze postać:

Element S [ i ] odwzorowuje liczbę niepodzielną ani przez 2, ani przez 3. Wartość liczby określamy w zależności od tego, czy indeks i  jest parzysty, czy nieparzysty:

Dodatkowe założenie obejmuje n, które powinno spełniać równanie:

Algorytm wyrzuca ze zbioru S wszystkie wielokrotności początkowych liczb. Najpierw wyznaczana jest pozycja kwadratu liczby na podstawie poprzedniej pozycji. W dalszej kolejności algorytm wyznacza wszystkie pozycje wielokrotności liczby począwszy od jej kwadratu i pomijając wielokrotności podzielne przez 2 lub przez 3. Elementy zbioru S znajdujące się na tych pozycjach są zaznaczane jako wyrzucone.

Algorytm jest trochę trudny do zrozumienia, lecz działa znakomicie. Zaletą jest zmniejszone 3 krotnie zapotrzebowanie na pamięć w stosunku do podstawowego algorytmu Euklidesa. Przy określaniu n  pamiętaj, iż musi spełniać zależność n  = 3m  + 2, gdzie m  jest liczbą nieparzystą. W przeciwnym razie ostatnie liczby pierwsze mogą nie zostać wyliczone.

Algorytm ulepszonego sita Eratostenesa

// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

program prg;

type
  Tbarray = array of boolean;

var
  c, k, t, q, m, n, i, j, ij : longword;
  S : Tbarray;

begin
  readln ( n );
  m := n div 3;
  if( m and 1 ) = 0 then inc ( m ); 
  setlength ( S, m+1 );
  c := 0; k := 1; t := 2;
  q := round ( sqrt ( n ) ) div 3;
  for i := 1 to m do S [ i ] := true;
  for i := 1 to q do
  begin
    k  := 3 - k;
    inc ( c, ( k shl 2 ) * i );
    j  := c;
    ij := ( i shl 1 ) * ( 3 - k ) + 1;
    inc ( t, k shl 2 );
    while j <= m do
    begin
      S [ j ] := false;
      inc ( j, ij );
      ij   := t - ij;
    end;
  end;
  write ( 2, ' ', 3, ' ' );
  for i := 1 to m do
    if S [ i ] then
    begin
      if( i and 1 ) = 1 then write ( 3 * i + 2 )
      else                   write ( 3 * i + 1 );
      write ( ' ' );
    end;
  writeln;
end.

// Sito Eratostenesa
// Data   : 17.03.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>
#include <cmath>

using namespace std;

int main( )
{
  unsigned int c, k, t, q, m, n, i, j, ij;
  bool * S;

  cin >> n;
  m = n / 3;
  if( ! ( m & 1 ) ) m++; 
  
  S = new bool [ m + 1 ];
  
  c = 0; k = 1; t = 2;
  q = ( ( unsigned int )sqrt ( n ) ) / 3;
  for( i = 1; i <= m; i++ ) S [ i ] = true;
  for( i = 1; i <= q; i++ )
  {
    k  = 3 - k;
    c  += ( k << 2 ) * i;
    j  = c;
    ij = ( i << 1 ) * ( 3 - k ) + 1;
    t  += k << 2;
    while( j <= m )
    {
      S [ j ] =  false;
      j    += ij;
      ij   =  t - ij;
    }
  }
  cout << 2 << " " << 3 << " ";
  for( i = 1; i <= m; i++ )
    if( S [ i ] )
    {
      if( i & 1 ) cout << 3 * i + 2;
      else        cout << 3 * i = 1;
      cout << " ";
    }
  cout << endl;
  delete [ ] S;
  return 0;
}

' Sito Eratostenesa
' Data   : 17.03.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

Dim As Uinteger c, k, t, q, m, n, i, j, ij

Input n
m = n \ 3
if( m And 1 ) = 0 Then m += 1 
  
Dim As Byte S ( m + 1 )

c = 0: k = 1: t = 2
q = Cuint ( Sqr ( n ) ) \ 3
For i = 1 To m: S ( i ) = 1: Next
For i = 1 To q
  k  =  3 - k
  c  += ( k Shl 2 ) * i
  j  =  c
  ij = ( i Shl 1 ) * ( 3 - k ) + 1
  t  += k Shl 2
  While j <= m
    S ( j ) = 0
    j    += ij
    ij   = t - ij
  Wend
Next
Print 2;" ";3;" ";
For i = 1 To m
  If S ( i ) = 1 Then
    if( i And 1 ) = 1 Then
      Print 3 * i + 2;
    Else
      Print 3 * i + 1;
    End If
    Print " ";
  End If
Next
Print
End