Отсеиваем простые из миллиарда чисел быстрее, чем в Википедии
(Источник рисунка )
Поэтому можно подумать, что за века этот алгоритм изучен вдоль и поперек и добавить к нему ничего невозможно. Общеизвестно, что Решето Эратосфена (РЭ) один из древнейших алгоритмов, появившийся задолго до изобретения компьютеров. Поэтому удивился, когда на днях случайно обнаружил, что вариант, который в Википедии преподносится как оптимальный, можно заметно оптимизировать.
Дело было так. Если посмотреть Википедию – там море ссылок на авторитетные источники, в которых запросто утонуть. Обладая более чем скудными знаниями по ФП, не берусь судить об ответах, но другие участники обсуждения забраковали некоторые из предложенных сходу решений, указав, что вместо теоретической сложности В обсуждении статьи по функциональному программированию (ФП) задал вопрос: как в этой парадигме написать РЭ.
(1)
предложенная реализация будет обладать вычислительной сложностью
(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 вычеркивать из РЭ числа
6 сек. Эта программа просеяла миллиард чисел за 17. 4 ГГц).
Сделав эту программу, я вдруг вспомнил о широко известных свойствах четных и нечетных чисел. на моем ПК (CPU Intel Core i7 3.
Лемма 1. 1) нечетное+нечетное=четное; 2) нечетное+четное=нечетное; 3) четное+четное= четное.
ЧТД.
2) не делится на 2 без остатка. Доказательство.
1) делится на 2. ЧТД. ЧТД.
3) делится на 2.
ЧТД. Лемма 2. Квадрат нечетного числа есть нечетное число.
Доказательство. не делится на 2 без остатка.
Замечание. Простое число, большее 2, нечетное.
Поэтому можно вычеркивать только нечетные числа:
(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. – почти вдвое быстрее.
Для оценки соответствия реального времени работы программы теоретическому я предположил, что
где – измеренное время работы;
– константа с размерностью времени;
– теоретическая оценка.
Для и меньше близко нулю. Вычислил отсюда для оценки (1) и (2). Поэтому привожу данные по первой программе для больших
Для второй программы наблюдается похожая картина. Как видим, оценка (1) гораздо ближе к реальным результатам. Скорее всего, авторы Википедии сами утонули в море информации по РЭ и пропустили эту работу. Сильно сомневаюсь, что открыл Америку с применением последовательности (3) и буду очень благодарен за ссылку на работу, где применялся этот подход.