Problem 241

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

More Reading
Newer// Problem 243