Problem 241

Project Euler Haskell

241 Perfection Quotients
Problem 241 - Project Euler
ひさしぶりのproject euler．
最近は忙しくて，やっていなかった（というか，頭の中になかった）．
が，偶然思いだし，とりかかる．

この問題は難しいと思われる．

  • 探索範囲が広い(10^18)
  • 計算量の算出が難しそう

いつもは，アルゴリズム考えて，計算量考えて，実装って流れ．しかし，今回はどうも，計算量が決定できない．
つまり，なんとなく，解けそうなんだけど，理論的保証がない．って，ことで，なかなか実装に移らなかった．
でも，実際やってみたら，まぁ，なんとか．


アルゴリズムのアイデアは，小さい範囲での探索からの観察．

あとは，約数の和と素因数分解の関係を考えた．素因数を組合せて解を構成するイメージをもつと，
\frac{\sigma(n)}{n}の分母，分子で余分な素数は約分される必要がある．

そこから紆余曲折をへて，できたhaskellのコード．
ひさしぶりに書いたから，どこか変かもしれない（が，そこは御愛嬌）

import Data.Function (on)
import Number (primes, factors, isPrime)
import Data.List (group, inits, findIndex)

getSieve :: Integer -> Integer -> [Integer] -> Maybe (Integer, Integer, Integer)
getSieve 1 2 [] = Just (2, 1, 3)
getSieve 1 _ _  = Nothing
getSieve 2 2 [] = Just (2, 1, 3)
getSieve a p ps
    | any ((< p).fst) $ qs = Nothing
    | elem pn ps           = Nothing
    | otherwise            = Just (pn, toInteger k, pm')
    where qs = map (\xs -> (head xs, length xs)).group.factors $ a
          (pn, k) = last qs
          pm' | p < pn    = p
              | otherwise = head [ q | q <- [pn+2..], isPrime q, not.elem q $ ps]

sigmaEq :: Integer -> Integer -> Integer -> Integer -> [Integer] -> [Integer]
sigmaEq 1 1 u p ps = [1]
sigmaEq a b u p ps
    | getSieve a p ps == Nothing = []
    | otherwise                  =
        concat $ do
          let Just (pn, i, p') = getSieve a p ps
          j <- takeWhile (\j -> pn^j <= u) [i..]
          let a' = a * ( (pn^(j+1) - 1) `div` (pn - 1))
              d  = gcd a' ( b * pn^j )
              a''= a' `div` d
              b' = (b * pn^j ) `div` d
              u' = u `div` ( pn^j )
          return.map (* pn^j ).sigmaEq a'' b' u' p' $ ps++[pn]

bound :: Integer -> Integer
bound u = floor.(2*).product $ [on (/) fromIntegral p (p-1) | p <- take k primes]
    where Just k = findIndex (>u).scanl1 (*) $ primes

p241 :: Integer -> Integer
p241 u = sum.(2:).concatMap (\k -> sigmaEq 2 k u 2 []) $ [5,7..bound u]

main :: IO ()
main = print.p241 $ 10^18

結局は\sigma(n)/n=b/aを再帰的に解いている．
デスクトップPC(Core 2 Duo E8500)で10秒ぐらい．マシンパワーを非常に頼っている気もするが，利用できるものは利用する方針で．
まぁ，どうせ，C++とかで書けば10倍ぐらい速くなるでs