今日は Probrem 21 をやってみた。
cf. http://projecteuler.net/index.php?section=problems&id=21
数 n に対して約数（nを含まない）すべてを足しあわせる関数 d があるとして，d(a) = b かつ d(b) = a となるような a と b （ただしa≠b）を amicable pair，それぞれを amicable number という（んだそうだ）。で，10000未満の amicable number をすべて足しあわせろ，という問題。
内包表記大活躍……なのはいいけど，すっげー時間がかかる。10000未満なんてやったら待ってられないので，下のコードでは1000未満にしてある。いい方法はないだろうか。
module Main (main) where
d :: Int -> Int
d n = sum [x | x <- [1..(n `div` 2)], n `mod` x == 0]
euler021 :: Int -> Int
euler021 n = sum $ map (\(a,b) -> a + b) pairs
where
pairs = [(a,b) | a <- [1..n], b <- [1..n], a < b, d a == b, d b == a]
main :: IO ()
main = print $ euler021 1000
^o^ >euler021.exe 504
1000未満では 220 と 284 のペアしかないみたいだ。
追記：
いい方法はないだろうか，と書いたらコメントやらトラックバックやらでいろいろ教えてくれた。ありがとうございます。
cf. http://haskell.g.hatena.ne.jp/nobsun/20080314/p1
おおっ，速い！
追記2：
ううっ，うっかりとエントリに上書きして消してしまった。できるだけ復元したけど元のエントリとは少し違ってるかも。
取り敢えず
d n = (1+) $ sum $ do
x <- takeWhile (￥x -> x*x <= n) [2..]
let (y,r) = n `divMod` x
guard $ r == 0
if (y == x) then [x] else [x,y]
で少し速くなります。
後は
・素因数分解してオイラー関数を使う
http://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%81%AE%CF%86%E9%96%A2%E6%95%B0
・メモ化
かな。
ああ。pairs が同じ計算を何度もやるからおそいですよ。
eular21 n = sum $ map (uncurry (+)) pairs
where pairs = [(a,b) | a <- [1..n], let b = d a, a /= b && a == d b ]
てな感じでいかが?
おっと間違えた
eular021 n = sum [ x+y | x <- [1..n], let y = d x, x < y && x == d y ]
かな
ついでに d も
d n = 1 + sum [x+y | x<- [2..floor (sqrt (fromInteger (toInteger n)))], let (y,r) = divMod n x, r == 0]
ああ。またまちがえてる。orz
n が平方数のとき d n が変になる。。。ゴメンナサイ。m(_._)m