Problem 126

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

立方体を数える簡単なお仕事、とはいかなく。高速化は結構大変だった。
でいろいろやった結果が下。

import Data.Array.IO
import Data.List

upper = 50000

layer (h,w,d) n = 4*(h+w+d+n-2)*(n-1) + 2*h*w + 2*h*d + 2*w*d

cuboid = [(h,w,d)| h<-[1..boundH], w<-[1..boundW h],d<-[1..boundD h w]]
         where boundH = (upper -2 ) `div` 4
               boundW h' = min h' $ (upper-2*h')`div`(2+2*h')
               boundD h' w' = min w' $ (upper -2*h'*w') `div` (2*h'+2*w')

count :: IO[(Int,Int)]
count = do count <- newArray (1,upper) 0 :: IO (IOUArray Int Int)
           let succArr ix = readArray count ix>>=writeArray count ix .succ
               addLayer n cb | layer cb n <= upper = do succArr $layer cb n
                                                        addLayer (n+1) cb
                             | otherwise = return ()
           mapM_ (addLayer 1) cuboid
           getAssocs count

main = count>>=print.find((==1000).snd)

上のほうが下のほうより倍ぐらい早い、理由はよく分からない。

count' :: IO[(Int,Int)]
count' = do count <- newArray (1,upper) 0 :: IO (IOUArray Int Int)
            let succArr ix = readArray count ix>>=writeArray count ix .succ
            mapM_ succArr $concatMap (takeWhile (<=upper).layers) cuboid
            getAssocs count
    where layers cb = map (layer cb) [1..]