241 Perfection Quotients

Problem 241 – Project Euler

ひさしぶりのproject euler．

最近は忙しくて，やっていなかった（というか，頭の中になかった）．

が，偶然思いだし，とりかかる．

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

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

いつもは，アルゴリズム考えて，計算量考えて，実装って流れ．しかし，今回はどうも，計算量が決定できない．

つまり，なんとなく，解けそうなんだけど，理論的保証がない．って，ことで，なかなか実装に移らなかった．

でも，実際やってみたら，まぁ，なんとか．


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

  • 解は少ない
  • kは小さい
  • 解を素因数分解すると比較的小さい素数ばかりが出現

あとは，約数の和と素因数分解の関係を考えた．素因数を組合せて解を構成するイメージをもつと，

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