30分プログラム 最大素因数分解の探索(Project Euler)

fromみずぴー日記

13195 の素因数は 5、7、13、29 である。

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

ポラード・ロー素因数分解法が面白そうだったのでチョイス。Gaucheで。

(use srfi-1)
(use srfi-42)
(use math.mt-random)

;;underの大きさは問題に応じて適当に調整してください。問題の数の素因数が全てunder以下なら解けます。
(define under 10000)
;;under以下の数をO(1)で素数判定するためのベクタ
(define under-vec (let ((v (make-vector under #f)))
		    (for-each (cut vector-set! v <> #t)
			      ;;under以下の素数のリスト
			      (list-ec (: n 2 under) (if (not (any?-ec (: j 2 (+ (floor (sqrt n)) 1)) (zero? (modulo n j))))) n))
		    v))

(define (prime-under? n)
  (and (< n under) (vector-ref under-vec n)))

(define mt (make <mersenne-twister> :seed (sys-time)))

;;ポラード・ロー素因数分解法
(define (p n)
  (let loop ((d 1))
    (cond 
     ((= d 1) (let ((x (mt-random-integer mt n))
		    (y (mt-random-integer mt n)))
		(loop (gcd (abs (- x y)) n))))
     ((= d n) #f)
     (else d))))

(define (find-largest-prime-factor n)  
  (let loop ((n n)
	     (l '())
	     (failed 0))
    (if (prime-under? n)
      ;;成功。最大の素因数を返す。
      (apply max (cons n l))
      (if (> failed 5)
	;;あきらめる。
	'giveup
	(let1 a (p n)
	  (if a
	    (begin
	      (format #t "~a = ~a * ~a.\n" n a (quotient n a))
	      (loop (quotient n a) (cons a l) 0))
	    (loop n l (+ failed 1))))))))

(find-largest-prime-factor 600851475143)
600851475143 = 71 * 8462696833.
8462696833 = 839 * 10086647.
10086647 = 1471 * 6857.
=> 6857
time は0.2秒@ThinkpadX60 でした。