Project Euler - Problem 3

問題

• 原文

  What is the largest prime factor of the number 600851475143 ?

• 日本語訳

  600851475143 の素因数のうち最大のものを求めよ。

解答

いきなり難易度が上がりました。素因数分解の問題です。 結局のところ片っ端から割ってみるしかないのですが、探索範囲を狭めることはできます。

正の整数の範囲で考えると、l = m * nかつm >= nであれば、floor(sqrt(l)) >= nであることが直感的に分かります。ここでfloorは実数を0の方向に丸めて整数にする関数、sqrtは平方根を返す実数関数です。 nが分かれば、mは除算で求められます。つまりlが与えられた時、それを2つの因数に分解するのに探索する範囲は高々2 .. floor(sqrt(l))で十分ということです。 2つに分けられたならこっちのもの、得られたmとnを再帰的に分解して、その結果を合わせれば素因数分解の結果が得られます。 アルゴリズムを書き下すと次のようになります：

1. n <- floor(sqrt(l))とする。
2. n = 1なら、これ以上は分解できない（lは素数である）のでlをそのまま返す。
3. nがlの因数なら、m <- l / nとし、mとnに対してこのアルゴリズムを再帰的に適用し、その結果を連結して返す。
4. 因数でないなら、n <- n - 1とし、2.から繰り返す。

これで計算は可能です。しかしもう少し最適化を検討してみましょう。

偶数は一般に2mの形で表せます。同様に奇数は2m + 1と書けます。 偶数同士の積は2m * 2n = 4mnとなり、l = 2mnとおけば2lなのでこれも偶数です。 奇数同士の積は(2m + 1)(2n + 1) = 4mn + 2(m + n) + 1であり、l = 2mn + m + nとすると2l + 1なのでこれも奇数です。 偶数と奇数の積は2m(2n + 1) = 4mn + 2mなので、l = 2mn + mとおいて2lなので偶数です。 つまり、奇数の因数は必ず奇数のみから成るということが分かります。先ほどのアルゴリズムでは、4.のステップでn <- n - 1としていました。これではlとnが両方とも奇数の時、1回不要な繰り返しが生じることになります。つまりこの場合を特別扱いして、n <- n - 2とすれば繰り返しを減らすことができます。 幸い、問題で与えられた600851475143という数値は奇数なので、この最適化の恩恵を十分に受けることができます。

この最適化を施したところ、40%ほど高速化しました。上の考察におけるlはn、nはdivという名前になっています。

(define (divisor-of? n m) (zero? (modulo m n)))
(define (factorize n)
  (define div (inexact->exact (floor (sqrt n))))
  (define (factorize-1 n div)
    (cond
     ((= div 1) (list n))
     ((and (odd? n) (even? div)) (factorize-1 n (- div 1)))
     ((divisor-of? div n) (append (factorize (/ n div)) (factorize div)))
     ((odd? n) (factorize-1 n (- div 2)))
     (else (factorize-1 n (- div 1)))))
  (factorize-1 n div))
(define (solve)
  (apply max (factorize 600851475143)))
(define (main argv)
  (display (solve))
  (newline))

追記

コメントの匿名氏のコードを基にしたところ500倍高速になりました。よって上記の議論は忘れて手続き的に書いた方が良さげ。 多値を扱うのにSRFI-8のreceiveマクロを使っています。

(define (divisor-of? m n) (zero? (modulo n m)))
(define (factorize n)
  (define (factorize-1 n i result)
    (cond
     ((< n (* i i)) (reverse (if (= n 1) result (cons n result))))
     ((divisor-of? i n) (factorize-1 (/ n i) i (cons i result)))
     (else (factorize-1 n (+ i 2) result))))
  (receive (n result)
           (let prepare-loop ((n n) (result '()))
             (if (divisor-of? 2 n) (prepare-loop (/ n 2) (cons 2 result))
                 (values n result)))
           (factorize-1 n 3 result)))
(define (solve)
  (apply max (factorize 600851475143)))
(define (main argv)
  (display (solve))
  (newline))