Progetto Eulero 23: approfondimento su questo programma stackoverflow necessario

Question

Ciao haskell utenti. Attualmente sto lavorando allo 23rd problem di Project Euler. Dove sono atm è che il mio codice mi sembra giusto - non nel significato di "buon algoritmo", ma nel significato "dovrebbe funzionare" - ma produce un sovraccarico di memoria dello stack.

So che il mio algoritmo non è perfetto (in particolare potrei sicuramente evitare di calcolare un risultato intermedio così grande ad ogni passo di ricorsione nella mia funzione worker).

Tuttavia, essendo in procinto di apprendere Haskell, mi piacerebbe capire perché questo codice fallisce così miseramente, al fine di evitare questo tipo di errori la prossima volta.

Qualsiasi giudizio sul motivo per cui questo programma è sbagliato sarà apprezzato.

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 

Modifica: l'errore esatto è Stack space overflow: current size 8388608 bytes.. L'aumento del limite di memoria dello stack tramite +RTS -K... non risolve il problema.

Edit2: sulla cosa sqrt, ho appena copiato incollato dal collegamento nei commenti. Per evitare di dover lanciare Integer a Doubles e affrontare i problemi di arrotondamento ecc ...

Answer 1

In futuro, è educato a tentare un po 'di minimalizzazione da soli. Ad esempio, con un po 'di gioco, sono stato in grado di scoprire che il seguente programma stack overflow anche (con uno stack 8M):

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

... che in realtà inchioda proprio quello che la funzione è avvitamento voi sopra . Diamo uno sguardo a worker:

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

Anche durante la mia prima lettura, questa funzione è stata interrotta dalla bandiera rossa nella mia mente, perché è ricorsiva in coda. La ricorsione della coda in Haskell non è generalmente una grande idea come in altre lingue; la ricorsione custodita (dove si produce almeno un costruttore prima di ricorsero, o si recita un po 'di volte prima di produrre un costruttore) è generalmente migliore per una valutazione lenta. E infatti, qui, ciò che sta accadendo è che ogni chiamata ricorsiva a worker sta creando un thunk nidificato più profondo e più profondo nell'argomento prev. Quando arriva il momento di restituire finalmente prev, dobbiamo andare molto in profondità in una lunga catena di chiamate Set.\\ per capire esattamente che cosa abbiamo finalmente.

Questo problema è leggermente offuscato dal fatto che l'ovvia rigida annotazione non aiuta. Facciamo il massaggio worker finché non funziona. La prima osservazione è che la prima clausola è completamente sostituita dalla seconda. Questo è stilistico; non dovrebbe influire sul comportamento (tranne che sulle liste vuote).

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

Ora, l'ovvio severità della nota:

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

Sono stato sorpreso di scoprire che questo ancora pila overflow! La cosa subdola è che lo seq sugli elenchi valuta solo quanto basta per sapere se l'elenco corrisponde a [] o _:_.Quanto segue non Stack Overflow:

import Control.DeepSeq 

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

non ho collegare questa versione finale di nuovo nel codice originale, ma almeno funziona con la minimizzato main sopra. A proposito, come si potrebbe la seguente idea di attuazione, che pila anche overflow:

import Control.Monad 

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

ma che può essere risolto utilizzando Data.Set invece di Data.List, e nessun annotazioni severità:

import Control.Monad 
import Data.Set as Set 

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

Answer 2

Ok, L'ho caricato e gli ho dato un colpo. Il consiglio di Daniel Wagner è abbastanza buono, probabilmente migliore del mio. Il problema è in effetti con la funzione worker, ma stavo suggerendo di usare Data.MemoCombinators per memorizzare la tua funzione.

Inoltre, l'algoritmo dei divisori è un po 'sciocco. C'è un modo molto migliore per farlo. È un po 'mistico e richiederebbe molto TeX, quindi ecco un link a una pagina di math.stackexchange su come farlo. Quello di cui stavo parlando è stata la risposta accettata, anche se qualcun altro fornisce una soluzione ricorsiva che penso possa correre più velocemente. (Non richiede fattorizzazione in numeri primi.)

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

Answer 3

Come Daniel Wagner correctly said, il problema è che

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

costruisce un tonfo mal nidificato. È possibile evitarlo e ottenere prestazioni leggermente migliori rispetto a deepseq sfruttando il fatto che entrambi gli argomenti su worker sono ordinati in questa applicazione. Così è possibile ottenere uscita incrementale notando che ad ogni passo tutto in prev piccolo di 2*a non può essere la somma di due numeri abbondanti, così

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

fa meglio. Tuttavia, è ancora negativo perché (\\) non può utilizzare l'ordinamento dei due elenchi. Se si sostituisce con

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 

(o utilizza la versione del pacchetto di data-ordlist), calcolando il set-differenza è O (lunghezza) invece di O (lunghezza^2).