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:

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

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

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

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:

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

Teraz w zbiorze pozostały liczby niepodzielne przez 2 i 3 – z wyjątkiem pierwszych liczb 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:

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

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

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

W zbiorze pozostały same liczby pierwsze.

{2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50}

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 wpisują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

K01: Dla i = 2, 3, …, n: ; zbiór początkowo zawiera
     wykonuj S[i] ← true ; wszystkie liczby
K02: g ← [sqr(n)] ; obliczamy granicę eliminowania
                  ; wielokrotności
K03: Dla i = 2, 3, …, g: ; w pętli wyrzucamy ze
     wykonuj kroki K04…K08 ; zbioru wielokrotności i
K04:   Jeśli S[i] = false, ; sprawdzamy, czy liczba
       to następny obieg pętli K03 ; i jest w zbiorze
K05:   w ← i2 ; jeśli tak, wyrzucamy jej
              ; wielokrotności ze zbioru
K06:   Dopóki w ≤ n:
       wykonuj kroki K07…K08
K07:   S[w] ← false
K08:   w ← w+i ; następna wielokrotność
K09: Dla i = 2, 3, …, n:
     wykonuj krok K10 ; przeglądamy zbiór wynikowy
K10:   Jeśli S[i] = true, ; wyprowadzając pozostałe
       to pisz i          ; w nim liczby
K11: Zakończ

Przykładowe programy

// 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
 
n += 1
Dim As Byte S(n)

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

# Sito Eratostenesa
# Data: 18.06.2024
# (C)2024 mgr Jerzy Wałaszek
#---------------------------

import math

n = int(input())
s = [True]*(n+1)
g = int(math.sqrt(n))

for i in range(2, g+1):
    if s[i]:
        w = i*i
        while w <= n:
            s[w] = False
            w += i
for i in range(2, n+1):
    if s[i]:
        print(i, end=" ")
print()

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

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.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

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:

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

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:

p0 = 3
q0 = 4

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 liczby 5 znajduje się na pozycji 12:

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

Zapiszmy:

p1 = p0+2 = 3+2 = 5
q1 = q0+8 = 4+8 = 12

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

Znów zapiszmy:

p2 = p1+ 2 =  5+ 2 = 7
q2 = q1+12 = 12+12 = 24

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

3

5

7

9

11

13

15

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

47

49

51

53

55

57

59

61

63

65

67

69

71

73

75

77

79

81

Znów zapiszmy:

p3 = p2+ 2 =  7+ 2 = 9
q3 = q2+16 = 24+16 = 40

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:

p4 = p3+ 2 =  9+ 2 = 11
q4 = q2+20 = 40+20 = 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):

Kolejne p powstaje z poprzedniego przez dodanie 2. Kolejne q powstaje z poprzedniego przez dodanie 2(p-1):

0

3

4

1

3+2 =  5

2×(5-1)  =  8

 4+ 8 = 12

2

5+2 =  7

2×(7-1)  = 12

12+12 = 24

3

7+2 =  9

2×(9-1)  = 16

24+16 = 40

4

9+2 = 11

2×(11-1) = 20

40+20 = 60

Zatem możemy zapisać wzór rekurencyjny:

p0 = 3, q0 = 4 : wartości startowe
Dla i > 0:
pi = pi-1+2 qi = qi-1+2(pi-1)

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

Algorytm ulepszonego sita Eratostenesa

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

K01: Jeśli n jest nieparzyste, ; n sprowadzamy
     to n ← n+1                ; do parzystego
K02: m ← n shr 1 ; m to połowa z n
K03: Dla i = 1, 2, …, m-1: ; w S są początkowo wszystkie
     wykonuj S[i] ← true ; liczby nieparzyste
K04: i ← 1 ; indeks kolejnych liczb liczb
           ; pierwszych w S
K05: p ← 3 ; odstęp 3, pierwsza wielokrotność 4 (9)
     q ← 4
K06: Jeśli S[i] = false, ; przeskakujemy liczby
     to idź do kroku K11 ; wyrzucone z S
K07: k ← q ; w k indeks pierwszej wielokrotności
           ; liczby 2i = 1
K08: Dopóki k < m, ; wyrzucamy z S wielokrotności
     wykonuj kroki K09…K10
K09:   S[k] ← false
K10:   k ← k+p ; wyznaczamy pozycję kolejnej wielokrotności, 
               ; przesuniętą o odstęp p
K11: i ← i+1   ; indeks następnej liczby w S
K12: p ← p+2   ; zwiększamy odstęp wielokrotności
K13: q ← q+(p shl 1)-2 ; wyznaczamy pierwszą wielokrotność
K14: Jeśli q < m, 
     to idź do kroku K06
K15: Pisz 2 ; początkowa liczba pierwsza
            ; wyprowadzana bezwarunkowo
K16: Dla i = 1, 2, …, m-1: ; przeglądamy S wypisując pozostałe
     wykonuj krok K17 ; w niej liczby pierwsze
K17:   Jeśli S[i] = true, 
       to pisz 2i+1
K18: Zakończ

Przykładowe programy

// 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.

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

' 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

# Sito Eratostenesa
# Data: 6.02.2024
# (C)2024 mgr Jerzy Wałaszek
# --------------------------

n = int(input())
if n & 1: n += 1
m = n >> 1
s = [True]*m
i, p, q = 1, 3, 4
while True:
    if s[i]:
        k = q
        while k < m:
            s[k] = 0
            k += p
    i += 1
    p += 2
    q += (p << 1)-2
    if q >= m: break
print(2, end=" ")
for i in range(1, m):
    if s[i]:
        print((i << 1)+1, end=" ")
print()

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

…

in

ip

…

m

5

7

11

13

17

19

23

25

29

31

35

37

41

43

47

49

53

55

59

61

…

3in+2

3ip+1

n

gdzie:
in – indeks nieparzysty, 
ip – indeks parzysty.

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:

S[in] → 3in+2
S[ip] → 3ip+1

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

n = 3m+2, gdzie m jest liczbą nieparzystą

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

K01: m ← (n-2) div 3 ; obliczamy liczbę elementów w S
K02: Dla i = 1, 2, …, m: ; w zbiorze S są początkowo
     wykonuj S[i] ← true ; wszystkie liczby
K03: c ← 0 ; c będzie wskazywało pozycję
           ; kwadratu kolejnej liczby w zbiorze
K04: k ← 1; t ← 2 ; elementy pomocnicze:
K05: q ← [sqr(n)/3] ; granica pozycji liczb, których
                    ; wielokrotności eliminujemy
K06: Dla i = 1, 2, …, q: ; w pętli wyznaczamy pozycje
     wykonuj kroki K07…K15 ; liczb do usunięcia
K07:   k ← 3-k
K08:   c ← c+4k×i ; pozycja kwadratu liczby o indeksie i
K09:   j ← c      ; do wyrzucania liczb posłuży indeks j
K10:   ij ← 2i×(3-k)+1
K11:   t ← t+4k
K12:   Dopóki j ≤ m, ; pętla usuwająca liczby
       wykonuj kroki K13…K15
K13:     S[j] ← false ; usuwamy j-tą liczbę
K14:     j ← j+ij ; obliczamy pozycję następnej
                  ; do usunięcia liczby
K15:     ij ← t-ij
K16: Pisz 2 3 ; dwie początkowe liczby pierwsze
              ; wyprowadzamy bezwarunkowo
K17: Dla i = 1, 2, …, m: ; przeglądamy zbiór S
     wykonuj kroki K18…K19
K18:  Jeśli S[i] = false, ; pomijamy liczby usunięte
       to następny obieg pętli K17
K19:   Jeśli i nieparzyste,  ; inaczej wyprowadzamy
       to pisz 2i+2      ; wartości liczb, które w zbiorze
       inaczej pisz 2i+1 ; pozostały
K20: Zakończ

Przykładowe programy

// Ulepszone 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.

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

' Ulepszone 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

# Ulepszone Sito Eratostenesa
# Data: 6.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ---------------------------

import math

n = int(input())
m = n//3
if m & 1 == 0: m += 1
s = [True]*(m+1)
c, k, t = 0, 1, 2
q = int(math.sqrt(n))//3
for i in range(1, q+1):
    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
print(2, 3, end=" ")
for i in range(1, m+1):
    if s[i]:
        if i & 1:
            print(3*i+2, end="")
        else:
            print(3*i+1, end="")
        print(" ", end="")
print()

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 101