Project Euler 23: требуется понимание этой программы переполнения стека

Question

Привет haskell молодцы. В настоящее время я работаю над 23-й проблемой из Project Euler . Я думаю, что мой код мне подходит - не в смысле «хороший алгоритм», а в смысле «должен работать», но вызывает переполнение стековой памяти.

Я знаю, что мой алгоритм не идеален (в частности, я бы, конечно, мог избежать вычисления такого большого промежуточного результата на каждом шаге рекурсии в моей worker функции).

Хотя, находясь в процессе изучения Хаскелла, я хотел бы понять, почему этот код так неудачно завершается, чтобы в следующий раз избежать подобных ошибок.

Любое понимание того, почему эта программа неправильна, будет оценено.

import qualified Data.List as Set ((\\))

main = print $ sum $ worker abundants [1..28123]

-- Limited list of abundant numbers
abundants :: [Int]
abundants = filter (\x -> (sum (divisors x)) - x > x) [1..28123]

-- Given a positive number, returns its divisors unordered.
divisors :: Int -> [Int]
divisors x  | x > 0     =   [1..squareRoot x] >>=
                            (\y ->  if      mod x y == 0
                                    then    let d = div x y in
                                            if y == d
                                            then [y]
                                            else [y, d]
                                    else    [])
            | otherwise = []


worker :: [Int] -> [Int] -> [Int]
worker (a:[]) prev = prev Set.\\ [a + a]
worker (a:as) prev = worker as (prev Set.\\ (map ((+) a) (a:as)))


-- http://www.haskell.org/haskellwiki/Generic_number_type#squareRoot
(^!) :: Num a => a -> Int -> a
(^!) x n = x^n

squareRoot :: Int -> Int
squareRoot 0 = 0
squareRoot 1 = 1
squareRoot n =
   let twopows = iterate (^!2) 2
       (lowerRoot, lowerN) =
          last $ takeWhile ((n>=) . snd) $ zip (1:twopows) twopows
       newtonStep x = div (x + div n x) 2
       iters = iterate newtonStep (squareRoot (div n lowerN) * lowerRoot)
       isRoot r  =  r^!2 <= n && n < (r+1)^!2
   in  head $ dropWhile (not . isRoot) iters

Редактировать: точная ошибка Stack space overflow: current size 8388608 bytes.. Увеличение лимита стековой памяти через +RTS -K... не решает проблему.

Edit2: о sqrt, я просто скопировал и вставил его по ссылке в комментариях. Чтобы избежать приведения Integer к Doubles и столкнуться с проблемами округления и т.д ...

Answer 1

В будущем будет вежливым попытаться немного минимизировать самостоятельно.Например, немного поиграв, я смог обнаружить, что следующая программа также переполняется стеком (со стеком 8M):

main = print (worker [1..1000] [1..1000])

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

worker (a:[]) prev = prev Set.\\ [a + a]
worker (a:as) prev = worker as (prev Set.\\ (map ((+) a) (a:as)))

Даже при первом прочтении эта функция была помечена красным, потому что она рекурсивна.Хвостовая рекурсия в Хаскеле, как правило, не такая хорошая идея, как в других языках;Охраняемая рекурсия (когда вы создаете хотя бы один конструктор перед рекурсией или повторяете несколько раз перед созданием конструктора), как правило, лучше для ленивой оценки.И на самом деле, здесь происходит то, что каждый рекурсивный вызов worker создает все более и более глубокий вложенный поток в аргументе prev.Когда приходит время, наконец, вернуть prev, нам нужно очень глубоко погрузиться в длинную цепочку Set.\\ вызовов, чтобы выяснить, что же у нас было в итоге.

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

worker []     prev = prev
worker (a:as) prev = worker as (prev Set.\\ map (a+) (a:as))

Теперь, очевидная аннотация строгости:

worker []     prev = prev
worker (a:as) prev = prev `seq` worker as (prev Set.\\ map (a+) (a:as))

Я был удивлен, обнаружив, что это все еще переполняет стек!Хитрость в том, что seq в списках оценивается достаточно далеко, чтобы узнать, соответствует ли список [] или _:_.Следующее не приводит к переполнению стека:

import Control.DeepSeq

worker []     prev = prev
worker (a:as) prev = prev `deepseq` worker as (prev Set.\\ map (a+) (a:as))

Я не включил эту финальную версию обратно в исходный код, но она, по крайней мере, работает с свернутым main выше.Кстати, вам может понравиться следующая идея реализации, которая также переполняет стек:

import Control.Monad

worker as bs = bs Set.\\ liftM2 (+) as as

, но которая может быть исправлена ​​с помощью Data.Set вместо Data.List и без аннотаций строгости:

import Control.Monad
import Data.Set as Set

worker as bs = toList (fromList bs Set.\\ fromList (liftM2 (+) as as))

Answer 2

Как Даниэль Вагнер правильно сказал , проблема в том, что

worker (a:as) prev = worker as (prev Set.\\ (map ((+) a) (a:as)))

создает плохо вложенный thunk. Вы можете избежать этого и получить несколько лучшую производительность, чем с deepseq, используя тот факт, что оба аргумента worker отсортированы в этом приложении. Таким образом, вы можете получить инкрементальный вывод, заметив, что на любом шаге все в prev меньше 2*a не может быть суммой двух чисел, поэтому

worker (a:as) prev = small ++ worker as (large Set.\\ map (+ a) (a:as))
  where
    (small,large) = span (< a+a) prev

делает лучше. Тем не менее, это все еще плохо, потому что (\\) не может использовать сортировку двух списков. Если вы замените его на

minus xxs@(x:xs) yys@(y:ys)
    = case compare x y of
        LT -> x : minus xs yys
        EQ -> minus xs ys
        GT -> minus xxs ys
minus xs _ = xs             -- originally forgot the case for one empty list

(или используйте версию пакета data-ordlist ), вычисляя разницу между наборами O (длина) вместо O (длина ^ 2).

Answer 3

Хорошо, я загрузил его и дал ему шанс. Совет Даниэля Вагнера довольно хорош, возможно, лучше, чем мой. Проблема действительно в рабочей функции, но я собирался предложить использовать Data.MemoCombinators для запоминания вашей функции.

Кроме того, ваш алгоритм делителей немного глупый. Есть намного лучший способ сделать это. Это немного математично и требует много TeX, так что вот ссылка на страницу math.stackexchange о том, как это сделать. Тот, о котором я говорил, был принятым ответом, хотя кто-то другой дает рекурсивное решение, которое, я думаю, будет работать быстрее. (Не требует первичной факторизации.)

https://math.stackexchange.com/questions/22721/is-there-a-formula-to-calculate-the-sum-of-all-proper-divisors-of-a-number