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

정의 (소수와 합성수)
    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은 합성수이다.
    

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

(*
 *  Filename: divideEach.ml
 *
 *  Purpose:  Determine whether the given integer is a prime or not.
 *
 *  Execute: ocaml nums.cma unix.cma divideEach.ml 1234567
 *
 *   Or
 *
 *  Compile: ocamlc -o divideEach.exe nums.cma unix.cma divideEach.ml 
 *  Execute: divideEach [integer]
 *
 *     Date:  2013. 2. 3.
 *   Author:  pkim __AT__ scripts.pe.kr
 *)

(*
  Execution Examples with OCaml:
      Prompt> divideEach 10006099720301
      10006099720301 = 2673067 * 3743303
      10006099720301 is a not prime.
      Elapsed time: 13.000000 sec.

      Prompt> divideEach 1234567812343
      1234567812343 = 1 * 1234567812343
      1234567812343 is a prime.
      Elapsed time: 6.000000 sec.

      Prompt> divideEach 9999994200000841
      9999994200000841 = 99999971 * 99999971
      9999994200000841 is a not prime.
      Elapsed time: 504.000000 sec.

      Prompt> divideEach 18446744073709551617
      18446744073709551617 = 274177 * 67280421310721
      18446744073709551617 is a not prime.
      Elapsed time: 1.000000 sec.

      Prompt> divideEach 10023859281455311421
      10023859281455311421 = 1308520867 * 7660450463
      10023859281455311421 is a not prime.
      Elapsed time: 5328.000000 sec.
*)

open Num ;;
open Unix ;;
open Printf ;;

let bignum n = num_of_int n ;;

let bignum_s n = num_of_string n ;;

let n = ref (bignum_s "10006099720301") in

(* Begin here *)
let cmdArgs = Sys.argv in
if ((Array.length cmdArgs) > 1) then begin
    n := bignum_s cmdArgs.(1) ;
end ;

let z = ref (quo_num !n (bignum 2)) in
if !n =/ (bignum 2) */ !z then
    printf "%s = %s * %s\n" (string_of_num !n) (string_of_num (bignum 2)) (string_of_num !z) ;

let time1 = fst ( mktime (localtime (time())) ) in

let d = ref (bignum 1) in
let k = ref (bignum 3) in
let cont_flag = ref true in
while (!k */ !k <=/ !n && !cont_flag) do
    (* z := !n // !k ; *)
    z := quo_num !n !k ;
    if !n =/ !k */ !z then (
        d := !k ;
        cont_flag := false
    )
    else
        k := !k +/ (bignum 2)
done;

let time2 = fst ( mktime (localtime (time())) ) in
let elapsed = (time2 -. time1)  in

printf "%s = %s * %s\n" (string_of_num !n) (string_of_num !d) (string_of_num (!n // !d)) ;

if !d =/ (bignum 1) then
    printf "%s is a prime.\n" (string_of_num !n)
else
    printf "%s is a not prime.\n" (string_of_num !n) ;

printf "Elapsed time: %f sec.\n" elapsed

 

 

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

(*
 *  Filename: pollardRho.ml
 *
 *  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: ocaml nums.cma unix.cma pollardRho.ml [integer]
 *
 *  Compile: ocamlc -o pollardRho.exe nums.cma unix.cma pollardRho.ml
 *  Execute: pollardRho [integer]
 *
 *     Date:  20013. 2. 3.
 *   Author:  PH Kim   [ pkim ((AT)) scripts.pe.kr ]
 *)

(*
  Execution Examples with OCaml:
      Prompt> pollardRho 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 = 1234567812343 * 1
      Elapsed time: 124.000000 sec

      Prompt> pollardRho 9999994200000841
      Try first the Pollard rho algorithm with c = 2
      d = 99999971, count = 3593
      9999994200000841 = 99999971 * 99999971
      Elapsed time: 0.000000 sec

      Prompt> pollardRho 18446744073709551617
      Try first the Pollard rho algorithm with c = 2
      d = 274177, count = 1028
      18446744073709551617 = 274177 * 67280421310721
      Elapsed time: 0.000000 sec

      Prompt> pollardRho 10023859281455311421
      Try first the Pollard rho algorithm with c = 2
      d = 1308520867, count = 20350
      10023859281455311421 = 1308520867 * 7660450463
      Elapsed time: 2.000000 sec
*)

open Unix ;;
open Printf ;;
open Num ;;

let bignum n = num_of_int n

let bignum_s n = num_of_string n

let f x c n =
    mod_num (x */ x +/ c)  n

let g x c n =
    f (f x c n) c n

let gcd x y =
    let a = ref (abs_num x) in
    let b = ref (abs_num y) in
    if !b =/ (bignum 0) then
        !a
    else begin
        while !b <>/ (bignum 0) do
            let t = mod_num !a !b in
            a := !b ;
            b := t
        done ;
        !a
    end ;;

let pollardRho n =
    let c = ref (bignum 2) in
    let x = ref (bignum 1) in
    let y = ref (bignum 1) in
    let d = ref (bignum 1) in
    let savedX = ref (!x +/ (bignum 0)) in
    let count = ref (bignum 0) in
    printf "Try first the Pollard rho algorithm with c = %s\n" (string_of_num !c) ;
    let cont_flag = ref true in
    while !d =/ (bignum 1) && !count */ !count <=/ n && !cont_flag do
        x := f !x !c n ;
        if !x =/ !savedX then begin
            printf "It is cyclic.  x = %s\n" (string_of_num !x) ;
            cont_flag := false
        end
        else begin
            y := g !y !c n ;
            d := gcd (abs_num (!x -/ !y)) n ;
        end ;
        count := !count +/ (bignum 1) ;
    done;

    printf "d = %s, count = %s\n" (string_of_num !d) (string_of_num !count) ;
    if !d >/ (bignum 1) && !d </ n then begin
        !d
    end
    else begin
        c := (bignum 3) ;
        x := (bignum 1) ;
        y := (bignum 1) ;
        d := (bignum 1) ;
        savedX := !x ;
        (* printf " (0.2) x = %s\n" (string_of_num !x) ; *)
        count := (bignum 0) ;
        printf "Try second the Pollard rho algorithm with c = %s\n" (string_of_num !c) ;
        cont_flag := true ;
        while !d =/ (bignum 1) && !count */ !count <=/ n do
            x := f !x !c n ;
            if !x =/ !savedX then begin
                 printf "It is cyclic.  x = %s\n" (string_of_num !x) ;
                 cont_flag := false ;
            end
            else begin
                 y := g !y !c n ;
                 d := gcd (abs_num (!x -/ !y)) n ;
                 count := !count +/ (bignum 1) ;
            end
        done ;

        printf "d = %s, count = %s\n" (string_of_num !d) (string_of_num !count) ;
        if !d >/ (bignum 1) && !d </ n then
            !d
        else begin
            c := (bignum 1) ;
            x := (bignum 1) ;
            y := (bignum 1) ;
            d := (bignum 1) ;
            savedX := !x ;
            count := (bignum 0) ;
            printf "Try third the Pollard rho algorithm with c = %s\n" (string_of_num !c) ;
            let cont_flag2 = ref true in
            while !d =/ (bignum 1) && !count */ !count <=/ n do
                x := f !x !c n ;
                if !x =/ !savedX then begin
                    printf "It is cyclic.  x = %s\n" (string_of_num !x) ;
                    cont_flag2 := false ;
                end
                else begin
                    y := g !y !c n ;
                    d := gcd (abs_num (!x -/ !y)) n ;
                    count := !count +/ (bignum 1) ;
                end
            done ;

            printf "d = %s, count = %s\n" (string_of_num !d) (string_of_num !count) ;
            !d
        end
    end ;;

(* Begin here *)
let n = ref (bignum 9991) in

let cmdArgs = Sys.argv in
if (Array.length cmdArgs > 1) then
    n := bignum_s cmdArgs.(1) ;

let time1 = fst ( mktime (localtime (time())) ) in
let k = pollardRho !n in
let time2 = fst ( mktime (localtime (time())) ) in
let elapsed = time2 -. time1 in

let z = !n // k in
if !n =/ k */ z then
    printf "%s = %s * %s\n" (string_of_num !n) (string_of_num k) (string_of_num z) ;

printf "Elapsed time: %f sec\n" elapsed ;