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)  :: 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)  :: 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..]