Решето Эратосфена

[egorius]

Ещё один пример ленивых вычислений и рекурсии, менее тривиальный. Есть такое решето Эратосфена — простой, но эффективный алгоритм нахождения простых чисел. Выписываем дцать первых целых чисел, начиная с двойки; дальше мы их будем вычёркивать. Двойку, как невычеркнутую, объявляем простым числом, и вычёркиваем все числа, кратные её (4, 6, 8 и так далее). Берём следующее ещё невычеркнутое число (это будет три), объявляем его простым и вычёркиваем кратные ему. И так повторяем до конца (достаточно дойти до корня из дцати — все оставшиеся невычеркнутые заведомо будут простыми).

Написать такой алгоритм на традиционном языке — дело нескольких минут, но опять-таки первый вопрос: какого размера заводить массив? А нельзя ли подойти к этой задаче функционально, в терминах множеств, а не алгоритмически? Попробуем.

Итак, множество простых чисел — это множество всех чисел, начиная с двойки, за вычетом множества составных чисел.
Множество составных чисел — это множество чисел, кратных простым числам.

И это не какая-нибудь тавтология, это — рекурсия. Достаточно явно назвать двойку простым числом, чтобы это определение завертелось, вычисляя всё новые и новые простые числа.

Дальше не для слабонервных.

В Хаскеле как таковом множеств нет, хотя они и есть в стандартной библиотеке. Только непонятно, как с ними дружат ленивые вычисления, поэтому сделаем себе сами множества на основе списка — одного из базовых типов языка. Нам нужны две операции: объединение и вычитание. Написать их просто, если считать, что множество хранится в виде упорядоченного по возрастанию списка:

union (a:as) (b:bs)
  | a < b   = a : (union as (b:bs))
  | a == b  = a : (union as bs)
  | a > b   = b : (union (a:as) bs)

minus (a:as) (b:bs)
  | a < b   = a : (minus as (b:bs))
  | a == b  = minus as bs
  | a > b   = minus (a:as) bs

Читается так: если голова первого списка (a) меньше головы второго списка (b), то в результате сначала пойдёт голова первого списка (a), а затем объединение остатка первого списка (as) со всем вторым (b:bs); ну и так далее. Единственный тонкий момент: это определение не говорит, что делать со списком из одного элемента или с пустым списком. Но нам это принципиально не нужно, мы собираемся работать только с бесконечными списками!

Теперь можно написать первую часть определения про простые числа:

primes = 2 : (minus [3..] composite)

Это дословно двойка (как мы говорили, это необходимая отправная точка для рекурсии) и разность множества чисел, начиная с тройки и дальше, и множества составных чисел.

С составными числами несколько сложнее:

composite = foldr union' [] (map multiple primes)
  where 
    union' (x:xt) ys = x : (union xt ys)
    multiple p = [ p*p, p*(p+1) .. ]

Здесь используется классический map-reduce: каждое простое число отображается (map) в множество кратных ему (multiple), а затем получившееся множество множеств сворачивается в одно множество (foldr) с помощью операции объединения (union'):

map multiple [2,3,5...] =
[
  [4,6..],
  [9,12..],
  [25,30..],
  ...
]

foldr union' [] (map multiple [2,3,5...]) =
union' [4,6..] ( union' [9,12..] ( union' [25,30..] ... ) ) =
[4,6,8,9,10,12,14,15,16,18,20,21,22,24,25,26,27,28,30...]

Тонкий момент состоит в том, что для свёртывания нельзя использовать операцию union в том виде, в каком мы её определили раньше, поскольку она будет пытаться вычислить оба своих аргумента, а реально известно только начало первого. Но поскольку мы знаем, что в нашем случае голова первого аргумента всегда меньше головы второго, то небольшая модификация (union') решает дело.

В итоге получается красивейшая, короткая, и притом довольно эффективная штука. Пользоваться ей просто:

> take 30 primes
[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,103,107,109,113]
(0.02 secs, 0 bytes)
> primes!!10000
104743
(1.91 secs, 259296496 bytes)

Проблема только в том, что на начальные потуги её написать, изучение мирового опыта, попытки понять предложенные в инете решения и написание окончательного варианта ушло дня три или четыре. Мозг заржавел, похоже. Кое-что я и до сих пор не понимаю, например, почему программа отказывается работать с foldr1 вместо foldr.

Кстати, ещё предстоит изучить The Genuine Sieve of Eratosthenes, чтобы понять, как эта задача решается ещё эффективнее на приоритетных очередях... Да и Кнута бы почитать надо...

Comments

[pigdeon.livejournal.com]

Мне эта тема как-то неблизка, но возможно, тебе будет интересна вот эта задача и ее обсуждение: http://petchik.livejournal.com/223706.html#cutid1