Problem 3

Project Euler

一週間，追記できなかったので，2008-03-23 - dachs_hippoの日記を仕切りなおし．

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

訳

13195の素因数は5，7，13，29である．

600851475143の最大の素因数は何か？

単純にはできないのは明らか．

[Haskell]「単純な方法」での失敗

「単純に」というのは，素数を列挙して下から順番にチェックするということです．これだとチェックする対象の数が大きいと膨大な無駄が生じます．

以下のコードは600851475143で実行すると，かなり悲惨です．

{-
e3.hs Project Euler Problem 3
大きい数字になると途端に遅くなるぞバージョン
-}
import System

primes::[Integer]
primes=2:(eratos [3,5..])
    where eratos (p:ps) = p: (eratos [ x | x<-ps, x `mod` p /= 0 ])


min_fact::Integer->Integer
min_fact n = head $ filter (\x -> (n `mod` x ==0)) primes


largest_factor::Integer->Integer
largest_factor n = last $ filter (nmod) $ filter (>= min) $ takeWhile (<= max) primes
    where min = min_fact n
          max = n `div` min
          nmod p = (n `mod` p == 0)

main = getArgs >>= putStr.show.largest_factor.read.head

primesはエラトステネスの篩をそのまましてるだけです．min_factorは最小の素因数を求めています．largest_factor（max_facorは避けました（笑））は min <= p <= n/min の素数を一個一個，割り切れるかどうかチェックして割り切れるものだけをリストにして，最後のそれの最後の値を取り出しています．

実は最初は

largest_factor::Integer->Integer
largest_factor n = last $ filter (nmod) $ takeWhile (<= n) primes
    where  nmod p = (n `mod` p == 0)

だったのですが，当然ながらあまりに遅いし，無駄が多すぎたので少し細工したのがe3.hsです．13195程度が相手なら多少の改善はありますが，600851475143が相手では無駄です．

[Haskell]中学一年生に戻る

中学一年生では素因数分解をやります．この問題はその計算に戻ればいいんです．
\begin{array} 5 & ) & 13195 \\\hline 7 & ) & 2639 \\\hline13 & ) & 377 \\\hline & & 17\end{array}
（mimeTeXのパーサ甘いなあ。。。もうちょっときちんとTeXで書けるようにならないかな）
最小の素因数を求めてそれで割って，その商に対して同じことをして，最後に素因数が出てきたら（最小の素因数が対象そのものになったら）終了するわけです．これなら対象がどんどん（等比数列的に）小さくなるので効率がよいでしょう．

{-
e3-2.hs Project Euler Problem 3
素直に素因数分解してみたぞバージョン
-}
import System

primes::[Integer]
primes=2:(eratos [3,5..])
    where eratos (p:ps) = p: (eratos [ x | x<-ps, x `mod` p /= 0 ])

min_fact::Integer->Integer
min_fact n = head $ filter (\x -> (n `mod` x ==0)) primes

factorization::Integer->[Integer]
factorization n 
    | n == p = [n]
    | otherwise = p :  factorization(n `div` p)
          where p = min_fact n

main = getArgs >>= putStr.show.last.factorization.read.head

却ってすっきりしました．実行結果は

   ___         ___ _
  / _ \ /\  /\/ __(_)
 / /_\// /_/ / /  | |      GHC Interactive, version 6.6.1, for Haskell 98.
/ /_\\/ __  / /___| |      http://www.haskell.org/ghc/
\____/\/ /_/\____/|_|      Type :? for help.

Loading package base ... linking ... done.
Prelude> :l e3-2.hs
[1 of 1] Compiling Main             ( e3-2.hs, interpreted )
Ok, modules loaded: Main.
*Main> (last.factorization) 600851475143
Loading package haskell98 ... linking ... done.
6857
*Main>

時間は一瞬とはいいがたいですが，耐えられないほどではないです（私の環境では体感で一秒くらい）．

しかし，ちょっと数字をずらして``600851475145''を相手にするともうだめです．``600851475145=5*120170295029''なので，結局下から順番に探すのと大差はなくなってしまいます．人間が計算するのと同様，素因数に極めて巨大な値が現われる場合は，当然ながら駄目です．

まじめにやるなら，既存の「素因数分解のアルゴリズム」を勉強しろとなるわけですね．

[Haskell] \rho法（rho法）を実装してみる

Wikipediaを参照して簡単に翻訳してみました．

{-
e3-3.hs Project Euler Problem 3
$\rho$法を単純に実装したぞバージョン
-}
import System

r:: Integer -> Integer -> Integer
r k n  = ((k*k) `mod` n) + 1

g::Integer -> Integer -> Integer -> Integer
g n a b = gcd (abs((r a n)- (r (r b n) n ))) n 

rho:: Integer -> Integer -> Integer -> Integer
rho n a b  | g n a b /= 1 = g n a b
           | otherwise    = rho n (r a n) (r (r b n) n)


factorization'::Integer->[Integer]
factorization' n 
        | n == p = [n]
        | otherwise = p :  factorization'(n `div` p)
    where p = rho n 1 1

main = getArgs >>= putStr.show.last.factorization'.read.head

結果は，かなり高速です．モンテカルロな方法（答えが明確にでるものに乱数的な（確率的な）要素を使って解を出す）なので，失敗することもあるはずですが，そういうところは目をつぶります．

*Main> factorization' 10023859281455311421
[7660450463,1308520867]
*Main> factorization' 600851475143
[71,839,1471,6857]
*Main> factorization' 600851475145
[5,120170295029]
*Main> factorization' 18446744073709551617
[274177,67280421310721]

最後の二つはそこそこの時間がかかりますが，それでも数分のレベルです．

しっかしmaximaとかはどうやってるんでしょうね．最適化されてるのは当然ですけども，すごく差が出ます（コードを見ろというのはごもっとも）．