Project Euler Problem 3 - 思いついたら書く日記 Ex

　問題3: 600851475143 の素因数のうち最大のものを求める。[=>原文和訳]

Emacsのversion

今回は微妙に環境依存の話なので、最初にうちの開発実行環境を書いておきます。

(version)
>"GNU Emacs 22.3.1 (i386-msvc-nt5.1.2600)
 of 2009-08-10 on KOICHIRO-PC"

(Meadow-version)
>"Meadow-3.02-dev (RINDOU)"

32bit版 WindowsXP上でMeadowという環境です。

28bitの壁

　32bit版 Emacs Lisp の整数は符号部を除いて 28bit。扱える最大値は 268435456 で今回の問題を解くには足りません。

(/ 600851475143 2)
> 0

　多倍長ライブラリも用意されていないのでライブラリ自作か？と思いきや抜け道がありました。ヒントはこれでした。

;; 28bit 最大値に1を加算
(1+ (string-to-number "fffffff" 16))
> -268435456  ;オーバーフロー

; もう１桁増やしたら当然だめだよね？
(1+ (string-to-number "10000000" 16)) 
> 268435457.0   ;おや？

; もう一度 28bit 最大値を取得して
(string-to-number "fffffff" 16)
> 268435455

(1+ 268435455.0) ;「.0」をつけると
> 268435456.0       ; 計算できた。

　何をいまさらと言われそうですが、浮動小数を使えば整数限界を超えた計算ができます。浮動小数にも精度限界はあり、大きすぎる数は下の桁が切り捨てられて正確な整数計算ができませんが、ある程度の大きさまでなら整数型同様に整数計算に使えます。

素数列挙

　ようやく本題。「素因数を求めよ」ですから、素数が必要になります。今回の問題は上限値が決まってますので「ふるい」でいいですよね。

(require 'cl)
(defun primes-sieve (N)
  "エラトステネスのふるい"
  (let* ((vec (make-vector N nil))
         (vec-len (length vec))
         (n 1) (i 1) (lis '(2)))
    (while (< (incf i 2) vec-len)
      (incf n 2)
      (if (null (aref vec i))
          (setq lis (cons n lis)))
      (let ((j i))
        (while (> vec-len (setq j (+ j n)))
          (aset vec j t))))
    (reverse lis)))

　"2"は初期値で確定済にして、3から奇数だけをチェックしています。これやるだけで4割くらい高速化されました。他にもインクリメントの位置や、マクロ（dotimes, incf）使用/非使用などの速度比較をして、一応最適と思われる形がこれでした。
　今回は「ふるい」には vector を使っています。Emacs Lisp には、bool-vector という、エラトステネスのふるいにジャストフィットしたデータ型が用意されていますが、速度の面では vector の方が若干速いようです。数十億の素数をふるいにかけるようなメモリ節約に迫られる状況なら bool-vector のほうが良いと思います。

問題３の解答

(require 'cl)
(defun problem003 (N)
  (let ((primes-list
	 (primes-sieve
	  (1+ (floor (sqrt N))))))
    (defun get-max-factor (n plist)
      (let ((p 0))
        (dolist (p plist n)
          (if (> (* p p) n)
              (return n))
          (when (zerop (mod n p))
            (while (zerop (mod n p))
              (setq n (/ n p)))
            (return
             (let ((m (get-max-factor n (cdr plist))))
               (if (> p m) p m)))))))
    (get-max-factor N primes-list)))

(problem003 600851475143.0)

　整数型では % で剰余が求められますが、浮動小数の場合 mod じゃないとだめなようです。promes-sieve は上で作った「ふるい」関数です。これはライブラリとして保存しておくことにしましょう。素因数分解ではNの平方根までの素数があれば計算できます。ただし求める素因数は Nの平方根より大きい場合があります。

実は素数はいらなかった（別解）

　他の人の解を見てああそうかと納得。最大の素因数を得るだけなら、Nより小さい値で割れなくなるまで割っていけばいいんですね。

(require 'cl)
(defun problem003-2 (N)
  (let ((n N) (i 2.0))
    (while (>= n (* i i))
      (while (zerop (mod n i))
	(setq n (/ n i)))
      (incf i))
    (if (= n 1) (1- i) n)))

(problem003-2  600851475143.0)

除数 i が被除数 n の平方根を超えたところで終了。ものすごくスッキリしたし、素数リストを作る手間も無いので速いです。

おまけ（素因数分解）

　今回の problem003をちょっと変形すると素因数分解プログラムになります。こっちの方が使い道がありそうです。

(defun get-factors (N)
  "素因数分解"
  (let ((primes-list (primes-sieve
                      (1+ (floor (sqrt N))))))
    (defun get-factors-sub (n plist)
      (cond
       ((= n 1) nil)
       ((null plist) (list n))
       (t (let ((p 0) (dest nil))
	    (dolist (p plist dest)
	      (if (> (* (float p) p) n)
		  (return (cons n dest)))
	      (when (zerop (mod n p))
		(while (zerop (mod n p))
		  (setq n (/ n p))
		  (setq dest (cons p dest)))
		(return
		 (append
		  dest (get-factors-sub
			n (cdr plist))))))))))
    (get-factors-sub N primes-list)))

　dolist の次行の終了判定は、pの二乗が n を超えないという判定ですが、pの二乗は簡単に28bit整数を超えてしまうので、floatに変換しています。

　素因数分解も、別解のやり方を使ったほうが数倍早いみたいです。

(defun get-factors2 (N)
  (let ((n N) (i 2.0) (dest nil))
    (while (>= n (* i i))
      (while (zerop (mod n i))
	(setq n (/ n i))
	(setq dest (cons i dest)))
      (incf i))
    (if (= n 1)
	dest
      (cons n dest))))

整数で収まる範囲でも素因数に浮動小数が入ってしまうので、それがいやな場合は、こんなフィルタを作るといいかも。

(defvar *int-max* (string-to-number "fffffff" 16))
(defun truncate-if-int (x)
  (if (> x *int-max*) x (truncate x)))

(mapcar #'truncate-if-int (get-factors2 6008514751431.0))