Probrem 21

今日は 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:

ううっ,うっかりとエントリに上書きして消してしまった。できるだけ復元したけど元のエントリとは少し違ってるかも。

カテゴリー: Haskell パーマリンク

5 Responses to Probrem 21

  1. _ のコメント:

    取り敢えず
    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
    ・メモ化
    かな。

  2. nobsun のコメント:

    ああ。pairs が同じ計算を何度もやるからおそいですよ。
    eular21 n = sum $ map (uncurry (+)) pairs
    where pairs = [(a,b) | a <- [1..n], let b = d a, a /= b && a == d b ]
    てな感じでいかが?

  3. nobsun のコメント:

    おっと間違えた
    eular021 n = sum [ x+y | x <- [1..n], let y = d x, x < y && x == d y ]
    かな

  4. nobsun のコメント:

    ついでに d も
    d n = 1 + sum [x+y | x<- [2..floor (sqrt (fromInteger (toInteger n)))], let (y,r) = divMod n x, r == 0]

  5. nobsun のコメント:

    ああ。またまちがえてる。orz
    n が平方数のとき d n が変になる。。。ゴメンナサイ。m(_._)m

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください