Project Euler 21

http://projecteuler.net/index.php?section=problems&id=21

約数を求めるために、エラトステネスのふるい的に素因数分解する。
その前に配列の更新方法。

import Data.Array
a = array (1, 5) [ (n, 1) | n <- [1..5] ]

update a = a // [(2, 2)]

main = print(update a)

//で更新できて、

array (1,5) [(1,1),(2,2),(3,1),(4,1),(5,1)]

と表示される。

import Data.Array
sieve f p | f!p == 1  = sieve (f // [ (m, p) | m <- a, f!m == 1 ]) (p + 1)

          | p * p > n = f

          | otherwise  = sieve f (p + 1)

                where a = if p == 2 then [p*2,p*3..n] else [p*3,p*5..n]
n = 10000

f = sieve (array (2, n) [ (k, 1) | k <- [2..n] ]) 2
div_pow n p | mod n p /= 0 = (n, 0)

            | otherwise    = (\(n, e) -> (n, e + 1)) (div_pow (div n p) p)
factorize 1 = []

factorize n | f!n == 1  = [ (n, 1) ]

            | otherwise = [(p, e)] ++ (factorize m)

                where

                    p = f!n

                    (m, e) = div_pow n (f!n)
divs [] = [1]

divs ((p,e):fs) = [ p^e1 * n | e1 <- [0..e], n <- (divs fs) ]
d n = sum (filter (/= n) (divs (factorize n)))
a = [ (p, d p) | p <- [2..n], (d p) < p && d (d p) == p ]

main = print(sum (map (\(p,q) -> p + q) a))