Хабрахабр

Отсеиваем простые из миллиарда чисел быстрее, чем в Википедии


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

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

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

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

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

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

Листинг 1

program EratosthenesSieve;
// Sieve of Eratosthenes
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.

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

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

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

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

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

ЧТД.
2)$ (2n+1)+(2m)=2n+2m+1$ не делится на 2 без остатка. Доказательство.
1) $(2n+1)+(2m+1)=2n+2m+2 $делится на 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$ – теоретическая оценка.

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

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

Показать больше

Похожие публикации

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»