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

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

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

(*
 *  Filename: DivideEach.fs
 *
 *  Purpose:  Determine whether the given integer is a prime or not.
 *
 *  Compile: fsc DivideEach.fs
 *  Execute: DivideEach [integer]
 *
 *     Date:  2010/07/17
 *   Author:  PH Kim   [ pkim ((AT)) scripts.pe.kr ]
 *)

(*
  Execution Examples:
      Prompt> DivideEach 1234567812343
      1234567812343 = 1 * 1234567812343
      1234567812343 is a prime
      Elapsed time: 0.750000 sec

      Prompt> DivideEach 9999994200000841
      9999994200000841 = 99999971 * 99999971
      9999994200000841 is a not prime
      Elapsed time: 121.140000 sec

      Prompt> DivideEach 18446744073709551617
      18446744073709551617 = 274177 * 67280421310721
      18446744073709551617 is a not prime
      Elapsed time: 0.218000 sec

      Prompt> DivideEach 10023859281455311421
      10023859281455311421 = 1308520867 * 7660450463
      10023859281455311421 is a not prime
      Elapsed time: 1471.500000 sec
/*)

#light

open System

let mutable n = 10006099720301I

// Begin here
let cmdArgs = System.Environment.GetCommandLineArgs()
if (Array.length cmdArgs > 1) then
    // n <- new Numerics.BigInteger(cmdArgs.[1])
    n <- (cmdArgs.[1]) |> Numerics.BigInteger.Parse

let mutable z = n / 2I
if n = 2I*z then
    printfn "%O = %O * %O" n  2I z

let time1 = DateTime.UtcNow

let mutable d = 1I
let mutable k = 3I
let mutable cont_flag = true
while (k*k <= n && cont_flag) do
    z <- n / k
    if n = k*z then
        d <- k
        cont_flag <- false
    else
        k <- k + 2I

let time2 = DateTime.UtcNow
let elapsed = float (int64 (time2 - time1).TotalMilliseconds) / 1000.0

printfn "%O = %O * %O" n d (n/d)
if d = 1I then
    printfn "%O is a prime" n
else
    printfn "%O is a not prime" n

printfn "Elapsed time: %f sec" elapsed

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

(*
 *  Filename: PollardRho.fs
 *
 *  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#
 *
 *  Compile: fsc PollardRho.fs
 *  Execute: PollardRho [integer]
 *
 *     Date:  20010/07/18
 *   Author:  PH Kim   [ pkim ((AT)) scripts.pe.kr ]
 *)

(*
  Execution Examples:
      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: 63.640000 sec

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

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

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

#light

open System

let f (x : bigint) (c : bigint) (n : bigint) =
    (x*x + c) % n

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

let gcd x  y =
    let mutable a = abs x
    let mutable b = abs y
    if b = 0I then
        a
    else
        while b <> 0I do
            let t = a % b
            a <- b
            b <- t
        a

let pollardRho (n : bigint) =
    let mutable c = 2I
    let mutable x = 1I
    let mutable y = 1I
    let mutable d = 1I
    let mutable savedX = x
    let mutable count = 0I
    printfn "Try first the Pollard rho algorithm with c = %O" c
    let mutable cont_flag = true
    while d = 1I && count*count <= n && cont_flag do
        x <- f x c n
        if x = savedX then
            printfn "It is cyclic.  x = %O" x
            cont_flag <- false
        else
            y <- g y c n
            d <- gcd (abs (x - y)) n
        count <- count + 1I

    printfn "d = %O, count = %O" d count
    if d > 1I && d < n then
        d
    else
        c <- 3I
        x <- 1I
        y <- 1I
        d <- 1I
        savedX <- x
        count <- 0I
        printfn "Try second the Pollard rho algorithm with c = %O" c
        cont_flag <- true
        while d = 1I && count*count <= n do
            x <- f x c n
            if x = savedX then
                 printfn "It is cyclic.  x = %O" x
                 cont_flag <- false
            else
                 y <- g y c n
                 d <- gcd (abs(x - y)) n
                 count <- count + 1I

        printfn "d = %O, count = %O" d count
        if d > 1I && d < n then
            d
        else
            c <- 1I
            x <- 1I
            y <- 1I
            d <- 1I
            savedX <- x
            count <- 0I
            printfn "Try third the Pollard rho algorithm with c = %O" c
            let mutable cont_flag2 = true
            while d = 1I && count*count <= n do
                x <- f x c n
                if x = savedX then
                    printfn "It is cyclic.  x = %O" x
                    cont_flag2 <- false
                else
                    y <- g y c n
                    d <- gcd (abs(x - y)) n
                    count <- count + 1I

            printfn "d = %O, count = %O" d count
            d

// Begin here
let mutable n = 9991I

let cmdArgs = System.Environment.GetCommandLineArgs()
if (Array.length cmdArgs > 1) then
    n <- (cmdArgs.[1]) |> Numerics.BigInteger.Parse

let time1 = DateTime.UtcNow
let k = pollardRho n
let time2 = DateTime.UtcNow
let elapsed = float (int64 (time2 - time1).TotalMilliseconds) / 1000.0

let z = n/k
if n = k*z then
    printfn "%O = %O * %O" n k z

printfn "Elapsed time: %f sec" elapsed