Problem 245

Project Euler Haskell

245 Coresilience
Problem 245 - Project Euler
オイラーのφ関数がらみの問題．
c(n) = (n - φ(n)) / (n - 1) = 1 / k となるような，合成数の和を求める．
探索範囲が 2 <= n <= 10^11 と広いが…
結論から言うと，全探索．
重要な点はどのように，探索範囲を狭くするか．
次のことが分かる（これぐらいしか分からなかった）：

  • n は奇数，
  • n は同じ素因数を2以上含まない，
  • k は n の最小素因数以下で偶数．

やったことは，n = p * q であるか，ないかで場合分け．
n = p * q のときはpを固定して，ある値の約数を列挙して，qを探した．
そうでないときは,最大の素因数を k を loop して，探した．

import Number (primes, isPrime, divisors)
import Data.List (findIndex)
import Control.Monad (guard)

main :: IO ()
main = print.sum.p245 $ 10^11

search :: Integer -> [Integer] -> [Integer]
search u [p] =
    do d <- takeWhile (<= div u p + p - 1).dropWhile (< 2 * p + 1).divisors $ p^2 - p + 1
       let q = d - p + 1
       guard $ isPrime q
       return $ p * q
search u ps =
    do k <- filter even [kL..minimum [kU1, kU2, p1 - 1]]
       let (p, r) = divMod (k * f + 1) (m - k * (m - f))
       guard $ r == 0 && isPrime p
       return $ m * p
    where m = product ps
          f = product.map pred $ ps
          (p1, p2) = (head ps, last ps)
          kL = div (m * p2 - 1) ((m - f) * p2 + f) + 1
          kU1 = div m (m - f)
          kU2 = div (m * (u - 1)) (m * f + u * (m - f))

candidates :: Integer -> Integer -> Integer -> [[Integer]]
candidates 1 _ _ = [[]]
candidates a l u = concat [map (p:).candidates (a-1) p $ div u p | p <- takeWhile (<=u').dropWhile (<=l) $ primes]
    where u' = floor.(** (1 / (fromInteger a))).fromInteger $ u

solutions :: Integer -> Integer -> [Integer]
solutions u a = concatMap (search u).candidates a 2 $ u

p245 :: Integer -> [Integer]
p245 u = map (sum.solutions u) $ [2..toInteger $ aU + 1]
    where Just aU = findIndex (>u).scanl (*) 1.tail $ primes

実行時間は1分切っているが，どうもマシンパワーに頼っているという，感じがしてならない．
(Core 2 Duo 3.16GHz)
まいどのことだが，どうせ C++ で実装したら100倍ぐらい速くなるんだって．多分．