День добрый, всем!
> Прошу прощения за некропостинг, но на письмо отвечу.
> Я остановился на 40 000, поскольку решето Эратосфена на Рефале
> оказалось сильно медленным, и его время растёт нелинейно от длины
> исходного списка. Причём, мне показалось, что даже квадратично.
> Поэтому результата в 176 100 я просто не дождусь.
> Время я замерял для переславской реализации Рефала-5. Рефал-5λ
> сильно медленнее (в разы), поэтому мне стыдно было показывать эти числа .
> В последующей переписке Антон Орлов продолжил оффтопик, мне уже лень
> пытаться переписывать его примеры на Рефал.
Вся переписка очень интересная и содержательная.
То, что нарыл Антон я ввел в свой курс по FP (Хаскелль) под заголовком
"Эратосфен это про вычеркивание, а не про делимость".
Действительно, Эратосфен говорил: оставь P как очередное простое и иди
с шагом P дальше и вычеркивай. То есть, он не про делимость, а про
удаление из изходного множества арифметической прогрессии:
[2P, 3P ..]
С этой точки зрения, чистый Эротосфен выглядит так (это все можно
описать на Рефале, но у меня нет Рефала):
primesE1 :: [Integer]
primesE1 = sieve [2..]
where sieve (p:xs) = let c1 = 2*p
c2 = c1 + p
in p : sieve (minusOS xs [c1, c2 ..])
> primesE1!!2000
> 17389
> (14831010 reductions, 27298448 cells, 2 garbage collections)
здесь minusOS -- вычитание "упорядоченных множеств", представленных
отсортированными ленивыми бесконечными списками:
minusOS :: [Integer] -> [Integer] -> [Integer]
minusOS (x:xs) (y:ys)
| x < y = x : minusOS xs (y:ys)
| x == y = minusOS xs ys
| x > y = x : minusOS xs ys
То, что нашел Антон, это следующих шаг: давайте вычеркивать будем не
как сказал Эратосфен, а вычеркнем ОБъЕДИНЕНИЕ всего, что советовал
вычеркивать Эратосфен. Появляется функция объединения бесконечного
множетсва "упорядоченных множеств", причем они перечислены в порядке
роста первого элемента:
primesE :: [Integer]
primesE = 2 : minusOS [3 ..] cs
where
cs = unionOS
[ [c1, c2 ..] | p <- primesE,
let c1 = 2*p,
let c2 = c1+p
]
unionOS :: [[Integer]] -> [Integer]
unionOS ((x:xs) : xss) = x : unionOS (ins xs xss)
where ins (x:xs) ((y:ys) : yss)
| x < y = (x:xs) : (y:ys) : yss
| x > y = (y:ys) : ins (x:xs) yss
| x == y = (y:ys) : ins xs yss
и этот вариант не лучше, а хуже (в 3+ раз) обычного Эратосфена:
> primesE!!2000
> 17393
> (47292483 reductions, 83601402 cells, 8 garbage collections)
А вот дальше можно сообразить, что первое число, которое удастся
простому P действительно вычеркнуть своей арифметической прогрессией
[2P, 3P..] будет P*P -- все меньшие члены прогрессии уже вычеркнуты до
него. Поэтому прогрессию можно заменить на [P*P, P*P+P ..]
В коде выше это замена всех c1 = 2*p на c1 = p*P
И это ускоряет обычный Эратосфен совсем немного, а Эратосфен с
объединением unionOS -- радикально и он на порядок (раз в 6) обгоняет
обычного!
> primesE1!!2000
> 17387
> (14692824 reductions, 27114235 cells, 2 garbage collections)
> primesE!!2000
> 17393
> (2731729 reductions, 4790869 cells)
Тупое прозрачное опредленеие "N -- простое, если не делется ни на одно
простое P, которое P*P<=K" дает лишь втрое хуже результат:
primes1 = 2 : filter (isPrime1) [3..]
where isPrime1 n = null(filter (\ p -> n`mod`p == 0)
(takeWhile (\ k -> (k*k <= n)) primes1))
> primes1!!2000
> 17393
> (7394774 reductions, 10355414 cells, 1 garbage collection)
И, наконец, мное описанное решето Абрамова-Дейкстры(?), которое не про
вычеркивание, не про делимость, а про взаимную простоту с
произведением простых (gcd pp n)==1 дает такой результат:
primesA = sieve 1 1 2 [2..]
where sieve pp n' lim ns@(n:ns1)
| n <= lim = n : sieve (pp*n) n lim ns1
| otherwise = sieve 1 1 (n'*n') (filter ((== 1).(gcd pp)) ns)
> primesA!!2000
> 17393
> (4254611 reductions, 7465819 cells)
что хуже ("всего" в 1.55 раз) Эратосфена вычеркивания с объединением
unionOS.
Это все.
Есть еще одна идея ускорения, но она применимы ко всем. И если
исследовать ее эффект, то применяя ее ко всем. В коде, который Антон
Орлов нашек на Haskell.org эта идея применена для No=1.
Идея состоит в том, что можно "в уме" вычеркнуть из [2..] все числа,
которые делятся на 2, 3, 5, ... p (взято No первых простых чисел).
Пусть ns результат вычеркивания. Тогда он определяется вот так:
Пусть pp произведение 2, 3, 5, ... p. Тогда:
ns = [ n+r | n <- [pp, 2*pp ...], r <- rs ]
where rs = filter ((==1).(gcd pp)) [p+1-pp .. p]
Если взять No простых, то после фильтрации в rs останется не pp
элементов, a меньше -- length rs. И от исходного множества у нас
останется вот такая доля: (length rs) / pp. Как-то так:
No length rs pp Доля
1 1 2 50.0%
2 2 6 33.3%
3 8 30 26.7%
4 48 210 22.9%
5 480 2310 20.8%
6 5760 30030 19.2%
7 92160 510510 18.1%
На мой вкус, разумно внедрять случай No=4 или No=5. Ускорение должно
быть вполне себе.
И, последнее, жалко, что не имею рефала. Можно все упомянутое
написать на рефале. Даже в условиях отсутствия ленивости, ее
эммуляция в этих алгоритмах на сложна.
С уважением,
Абрамов С.М.
ab...@botik.ru
мобильный: +7(903)2928308
