Pollard's rho method 소개: 정수의 인수분해(factorizing integers) with Common Lisp

정의 (소수와 합성수)
    1보다 큰 양의 정수 n에 대하여
    (i) n = a * b 를 만족하고 1보다 큰 두 양의 정수 a와 b가 존재하면,
        n을 합성수(合成數, composite number)라고 한다. 참고로 이 경우,
        a, b 중에 적어도 하나는 sqrt(n) 보다 작거나 같다.
        합성수의 예로는 4, 6, 9, 24, 143 등이 있다.
    (ii) n = a * b 를 만족하고 1보다 큰 두 양의 정수 a와 b가 존재하지 않으면,
         즉 n을 두 양의 정수의 곱으로 표현하는 방법이 1*n과 n*1 두 가지 뿐이면,
         n을 소수(素數, prime number)라고 한다.  소수의 예로는 2, 3, 5, 7, 11 등이 있다.
         n이 소수인지 아닌지 확인하려면, 
         n을 2 보다 크거나 같고 sqrt(n) 보다 작거나 같은 모든 정수로 나누어 본다.
         이 경우 언제나 나누어 떨어지지  않으면 n은 소수이고, 그렇지 않으면 n은 합성수이다.
    

우선 다음의 Common Lisp 소스 코드는 명령행 인자로 전달 받은 양의 정수 n을
2 및 3, 5, 7, ... , (sqrt n) 이하의 홀수들로 나누어 보아 n이 합성수인지 아닌지
확인하는 애플리케이션 소스이다. 확인하는데 걸린 경과 시간도 알려준다.

;;  Filename: divideEach.lsp
;;
;;  Purpose:  Determine whether the given integer is a prime or not.
;;
;;  Execute: clisp divideEach.lsp [integer]
;;
;;     Date:  2013. 9. 6.

#|
  Execution Examples:
      Prompt> clisp divideEach.lsp 1234567812343
      1234567812343 = 1 * 1234567812343
      1234567812343 is a prime.
      Elapsed time: 3 secs
 
      Prompt> clisp divideEach.lsp 9999994200000841
      9999994200000841 = 99999971 * 99999971
      9999994200000841 is a not prime
      Elapsed time: 51.437000 sec
 
      Prompt> clisp divideEach.lsp 18446744073709551617
      18446744073709551617 = 274177 * 67280421310721
      18446744073709551617 is a not prime
      Elapsed time: 0.141000 sec

      Prompt> clisp divideEach.lsp 10023859281455311421
      10023859281455311421 = 1308520867 * 7660450463
      10023859281455311421 is a not prime
      Elapsed time: 705.016000 sec
|#

(setf n 10006099720301)
(if (> (length ext:*args*) 0)
    (setf n (parse-integer (nth 0 ext:*args*)))
)

(setf z (floor (/ n 2)))
(if (= n (* 2 z))
    (progn
        (format t "~A = ~A * ~A~d" n 2 z)
        (quit)
    )
)

(setf time1 (get-universal-time))

(setf d 1)
(setf k 3)

(loop while (<= (* k k) n) do
    (setf z (floor (/ n k)))
    (if (= n (* k z))
        (progn
            (set d k)
            (return)
        )
    )
    (setf k (+ k 2))
)

(setf z (floor (/ n d)))

(setf time2 (get-universal-time))

(format t "~A = ~A * ~A~%" n d (floor (/ n d)))

(if (= d 1)
    (format t "~A is a prime.~%" n)
    (format t "~A is a not prime.~%" n)
)

(format t "Elapsed time: ~A secs" (- time2 time1))

 

이제 다음은 정수의 인수분해 능력이 뛰어난  Pollard의 rho 방법(일명 거북이와 토끼 알고리즘, tortoise-hair algorithm)을 구현한  Common Lisp 소스 코드이다. 이 알고리즘은 소수에 대해서는 시간이 많이 걸리지만, 합성수에 대해서는 시간이 매우 적게 걸린다는 것이 특징이다.

 

;;  Filename: pollardRho.lsp
;;
;;  Purpose:  By using the pollard rho method,
;;            determine whether the given integer is a prime or not.
;;
;;      See:  http://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
;;            http://en.wikipedia.org/wiki/Floyd%27s_cycle-finding_algorithm#Tortoise_and_hare#
;;
;;  Execute: clisp pollardRho.lsp [integer]
;;
;;     Date:  2013. 9. 6.

#|
  Execution Examples:
      Prompt> clisp pollardRho.lsp 1234567812343
      Try first the Pollard rho algorithm with c = 2
      d = 1234567812343, count = 466951
      Try second the Pollard rho algorithm with c = 3
      d = 1, count = 1111112
      Try third the Pollard rho algorithm with c = 1
      d = 1234567812343, count = 799441
      1234567812343 = 1 * 1234567812343
      Elapsed time: 104 secs      Try first the Pollard rho algorithm with c = 2

      Prompt> clisp pollardRho.lsp 9999994200000841
      Try first the Pollard rho algorithm with c = 2
      d = 99999971, count = 3593
      9999994200000841 = 99999971 * 99999971
      Elapsed time: 0 secs

      Prompt> clisp pollardRho.lsp 18446744073709551617
      Try first the Pollard rho algorithm with c = 2
      d = 274177, count = 1028
      18446744073709551617 = 274177 * 67280421310721
      Elapsed time: 0 secs

      Prompt> clisp pollardRho.lsp 10023859281455311421
      Try first the Pollard rho algorithm with c = 2
      d = 1308520867, count = 20350
      10023859281455311421 = 1308520867 * 7660450463
      Elapsed time: 1 secs
|#

(defun f(x c n)
    (floor (mod (+ (* x x) c) n))
)

(defun g(x c n)
    (f (f x c n) c n)
)

(defun  my-gcd(x y):
 (let* ((a (abs x))
        (b (abs y))
        (tt 0))
        (if (zerop b) a)
        (loop while (not (zerop b)) do
            (setf tt (floor (mod a b)))
            (setf a b)
            (setf b tt)
        )
        a
    )
)

(defun pollardRho(n)
 (let* ((c 2)
        (x 1)
        (y 1)
        (d 1)
        (savedX x)
        (count 0)
        (tt 0))

     (setf c 2)
     (setf x 1)
        (setf y 1)
     (setf d 1)
     (setf savedX x)
     (setf count 0)

        (format t "Try first the Pollard rho algorithm with c = ~A~%" c)
        (loop while (and (= d 1) (<= (* count count) n)) do
            (setf x (f x c n))
            (if (= x savedX)
                (progn
                    (format t "It is cyclic.  x = ~A~%" x)
                    (return)   ;; break loop
                )
            )
            (setf y (g y c n))
            (setf d (my-gcd (abs (- x y)) n))
            (setf count (1+ count))
        )

        (format t "d = ~A, count = ~A~%" d count)
        (if (and (> d 1) (< d n))
                d
      (progn
     (setf c 3)
     (setf x 1)
        (setf y 1)
     (setf d 1)
     (setf savedX x)
     (setf count 0)

        (format t "Try second the Pollard rho algorithm with c = ~A~%" c)
        (loop while (and (= d 1) (<= (* count count) n)) do
            (setf x (f x c n))
            (if (= x savedX)
                (progn
                    (format t "It is cyclic.  x = ~A~%" x)
                    (return)   ;; break loop
                )
            )
            (setf y (g y c n))
            (setf d (my-gcd (abs (- x y)) n))
            (setf count (1+ count))
        )

        (format t "d = ~A, count = ~A~%" d count)
        (if (and (> d 1) (< d n))
                d

      (progn
     (setf c 1)
     (setf x 1)
        (setf y 1)
     (setf d 1)
     (setf savedX x)
     (setf count 0)

        (format t "Try third the Pollard rho algorithm with c = ~A~%" c)
        (loop while (and (= d 1) (<= (* count count) n)) do
            (setf x (f x c n))
            ; (format t "  x = ~A~%" x)
            (if (= x savedX)
                (progn
                    (format t "It is cyclic.  x = ~A~%" x)
                    (return)   ;; break loop
                )
            )
            (setf y (g y c n))
            (setf d (my-gcd (abs (- x y)) n))
            (setf count (1+ count))
        )

        (format t "d = ~A, count = ~A~%" d count)
        (if (and (> d 1) (< d n))
                d

         1  )
      ))
      ))
    )
)

(setf n 9991)
(if (> (length ext:*args*) 0)
    (setf n (parse-integer (nth 0 ext:*args*)))
)

(setf time1 (get-universal-time))

(setf k (pollardRho n))

(setf time2 (get-universal-time))

(setf z (floor (/ n k)))
(if (= n (* k z))
    (format t "~A = ~A * ~A~%" n k z)
)

(format t "Elapsed time: ~A secs" (- time2 time1))