(Источник рисунка )

Общеизвестно, что Решето Эратосфена (РЭ) один из древнейших алгоритмов, появившийся задолго до изобретения компьютеров. Поэтому можно подумать, что за века этот алгоритм изучен вдоль и поперек и добавить к нему ничего невозможно. Если посмотреть Википедию – там море ссылок на авторитетные источники, в которых запросто утонуть. Поэтому удивился, когда на днях случайно обнаружил, что вариант, который в Википедии преподносится как оптимальный, можно заметно оптимизировать.

Дело было так. В обсуждении статьи по функциональному программированию (ФП) задал вопрос: как в этой парадигме написать РЭ. Обладая более чем скудными знаниями по ФП, не берусь судить об ответах, но другие участники обсуждения забраковали некоторые из предложенных сходу решений, указав, что вместо теоретической сложности

$O(n \log \log n)$ (1)

предложенная реализация будет обладать вычислительной сложностью

$O(n^2/\log n) $ (2)

и что с такой сложностью не дождаться, когда, например, просеются 10 миллионов чисел. Мне стало интересно и я попробовал реализовать оптимальную согласно Википедии версию, используя привычное мне процедурное программирование. В Delphi-7 у меня получился следующий код:

Листинг 1
program EratosthenesSieve;
// Sieve of Eratosthenes
{$APPTYPE CONSOLE}

uses
  SysUtils, DateUtils,Math;
const
 version ='1.0.1d1';
 N = 1000000000; // number of primes
var
 sieve : array [2..N] of boolean; // false if prime
 t0, t1,dt : TDateTime;
 O,C : Extended;

procedure init;
  var
   i : integer;
  begin
   for i:=2 to n do
    sieve [i] := false;
  end; //init

procedure calc (start : integer);
  var
   prime, i : integer;
   breakLoop, exitProc : Boolean;
  begin
   prime := start;
   exitProc := false;
   repeat
// find next prime
      prime := prime+1;
      while (prime<N) and sieve[prime] do
       inc (prime);
      i:= sqr(prime);
// delete
      if i<= N then
        begin
          breakLoop := false;
          repeat
            if i<= N then
             begin
              sieve [i] := true;
              inc (i,prime);
             end
            else  // if i<= N
             breakLoop := true;
          until breakLoop;
        end  
      else // if prime+prime<= N
       exitProc := true;
   until exitProc;
  end; //calc

procedure print;
  var
    i :integer;
    found : integer;

  begin
    found := 0;
    for i:=2 to N do
     if not sieve [i] then
      begin
//       write (i,', ');
        inc(found);
      end;
    writeln;
    writeln ('Found ',found,' primes.');
  end; //

begin // program body
  writeln ('Sieve of Eratosthenes version ', version);
  writeln('N= ',N);

  init;
  t0 := now;
  writeln('Program started ',DateTimeToStr(t0));
  t0 := now;
  calc (1);
  t1 := now;
  writeln('Program finished ',DateTimeToStr(t1));
  dt := SecondSpan(t1,t0);
  writeln ('Time is ',FloatToStr(dt),' sec.');

  O := N* ln(ln(N));
  C := dt/O;
  writeln ('O(N ln ln N)= ',O,' C=',C);

  O := sqr(N)/ln(N);
  C := dt/O;
  writeln ('O(N*N/ln N)= ',O,' C=',C);

  print;
  
  writeln ('Press Enter to stop...');
  readln;
end.


РЭ представлено булевым массивом sieve с инверсными значениями — если число простое, оно обозначается как false, что позволяет сократить количество операций отрицания (not) при просеивании. В программе 3 процедуры: инициализации РЭ — init, вычислений (просеивание и зачеркивание чисел в РЭ) — calc и вывода найденных в результате простых чисел — print, при этом подсчитывается количество найденных чисел. Особо обращу внимание на закомментированный вывод простых чисел в процедуре print: для тестирования при N=1000 комментарий снимается.

Здесь в процедуре calc использую рекомендацию Википедии: для очередного простого числа i вычеркивать из РЭ числа

$i^2, i^2+i, i^2+2i, …$



Эта программа просеяла миллиард чисел за 17.6 сек. на моем ПК (CPU Intel Core i7 3.4 ГГц).
Сделав эту программу, я вдруг вспомнил о широко известных свойствах четных и нечетных чисел.

Лемма 1. 1) нечетное+нечетное=четное; 2) нечетное+четное=нечетное; 3) четное+четное= четное.

Доказательство

1) $(2n+1)+(2m+1)=2n+2m+2 $делится на 2. ЧТД.
2)$ (2n+1)+(2m)=2n+2m+1$ не делится на 2 без остатка. ЧТД.
3)$ (2n)+(2m)=2n+2m$ делится на 2. ЧТД.

Лемма 2. Квадрат нечетного числа есть нечетное число.
Доказательство. $(2n+1)^{2}=4n^{2}+4n+1$ не делится на 2 без остатка. ЧТД.

Замечание. Простое число, большее 2, нечетное.

Поэтому можно вычеркивать только нечетные числа:

$i^2, i^2+2i, i^2+4i, … $ (3)

Но прежде надо вычеркнуть все четные числа. Это делает измененная процедура инициализации init.

Листинг 2
program EratosthenesSieve;
// Sieve of Eratosthenes
{$APPTYPE CONSOLE}

uses
  SysUtils, DateUtils,Math;
const
 version ='1.0.1d1';
 N = 1000000000; // number of primes
var
 sieve : array [2..N] of boolean; // false if prime
 t0, t1,dt : TDateTime;
 O,C : Extended;

procedure init;
  var
   i : integer;
  begin
   for i:=2 to n do
    sieve [i] := not odd(i); 
  end; //init

procedure calc (start : integer);
  var
   prime,prime2, i : integer;
   breakLoop, exitProc : Boolean;
  begin
   prime := start;
   exitProc := false;
   repeat
// find next prime
      prime := prime+1;
      while (prime<N) and sieve[prime] do
       inc (prime);
//      i:= prime*prime;
      i:= sqr(prime);
      prime2 := prime+prime;   
// delete
      if i<= N then
        begin
          breakLoop := false;
          repeat
            if i<= N then
             begin
              sieve [i] := true;
              inc (i,prime2);
             end
            else  // if i<= N
             breakLoop := true;
          until breakLoop;
        end  
      else // if prime+prime<= N
       exitProc := true;
   until exitProc;
   sieve [2] := false;
  end; //calc

procedure print;
  var
    i :integer;
    found : integer;

  begin
    found := 0;
    for i:=2 to N do
     if not sieve [i] then
      begin
//      write (i,', ');
        inc(found);
      end;
    writeln;
    writeln ('Found ',found,' primes.');
  end; //

begin // program body
  writeln ('Sieve of Eratosthenes version ', version);
  writeln('N= ',N);

  init;
  t0 := now;
  writeln('Program started ',DateTimeToStr(t0));
  t0 := now;
  calc (2);
  t1 := now;
  writeln('Program finished ',DateTimeToStr(t1));
  dt := SecondSpan(t1,t0);
  writeln ('Time is ',FloatToStr(dt),' sec.');

  O := N* ln(ln(N));
  C := dt/O;
  writeln ('O(N ln ln N)= ',O,' C=',C);

  O := sqr(N)/ln(N);
  C := dt/O;
  writeln ('O(N*N/ln N)= ',O,' C=',C);

  print;
  
  writeln ('Press Enter to stop...');
  readln;
end.


Эта программа сработала за 9.9 сек. – почти вдвое быстрее.

Для оценки соответствия реального времени работы программы теоретическому я предположил, что

$dt=C*O,$



где $dt$ – измеренное время работы;
$C $– константа с размерностью времени;
$O$ – теоретическая оценка.

Вычислил отсюда $C $ для оценки (1) и (2). Для $N= 10^6$ и меньше $dt$ близко нулю. Поэтому привожу данные по первой программе для больших $N.$

$N $ (1) (2)
$10^7 $ $1.69\cdot 10^{-9} $ $2.74\cdot 10^{-9}$
$10^8 $ $5.14 \cdot 10^{-9}$ $1.47 \cdot 10^{-8}$
$10^9$ $5.80 \cdot 10^{-9}$ $1.29 \cdot 10^{-7}$

Как видим, оценка (1) гораздо ближе к реальным результатам. Для второй программы наблюдается похожая картина. Сильно сомневаюсь, что открыл Америку с применением последовательности (3) и буду очень благодарен за ссылку на работу, где применялся этот подход. Скорее всего, авторы Википедии сами утонули в море информации по РЭ и пропустили эту работу.

PS О приведенном в Википедии алгоритме с «линейным временем работы» см.