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

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

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

/*
 * Filename: divideEach.groovy
 *
 *  Purpose: Determine whether the given integer is a prime or not.
 *
 *  Execute: groovy divideEach.groovy [integer]
 *
 *     Date: 2009/03/23
 *   Author: PH Kim   [ pkim ((AT)) scripts.pe.kr ]
 */

/*
  Execution Examples (with Groovy 1.6.0):

      Prompt> groovy divideEach.groovy 1234567812343
      1234567812343 = 1 * 1234567812343
      1234567812343 is a prime
      Elapsed time: 4.953 sec
 
      Prompt> groovy divideEach.groovy 9999994200000841
      9999994200000841 = 99999971 * 99999971
      9999994200000841 is a not prime
      Elapsed time: 468.688 sec

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

      Prompt> groovy divideEach.groovy 10023859281455311421
      10023859281455311421 = 1308520867 * 7660450463
      10023859281455311421 is a not prime
      Elapsed time: 3521.562000 sec
*/

BigInteger n = 10006099720301
if (args.length > 0) {
    n = new BigInteger(args[0])
}

BigInteger z = n / 2G    // z = n / 2
if (n == 2G*z) {
    print "${n} = 2 * ${z}\n"
}

def time1 = new GregorianCalendar().getTime().getTime()

BigInteger d = 1
BigInteger k = 3
while (k*k <= n) {
    z = n / k    // z = n / k
    if (n == k*z) {
        d = k
        break
    }
    k = k + 2G
}

def time2 = new GregorianCalendar().getTime().getTime()

print "$n = $d * ${n/d}\n"
if (d == 1)
    print "$n is a prime\n"
else
    print "$n is a not prime\n"

print "Elapsed time: ${(time2 - time1)/1000.0} sec\n"

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

/*
 * Filename: pollardRho.groovy
 *
 *  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: groovy pollardRho.groovy [integer]
 *
 *     Date: 2009/03/24
 *   Author: PH Kim   [ pkim ((AT)) scripts.pe.kr ]
 */

/*
  Execution Examples (with Groovy 1.6.0):

      Prompt> groovy pollardRho.groovy 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: 33.937 sec

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

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

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

BigInteger f(BigInteger x, BigInteger c, BigInteger n) {
    return (x*x + c) % n
}

BigInteger g(BigInteger x, BigInteger c, BigInteger n) {
    return f(f(x, c, n), c, n)
}

BigInteger gcd(BigInteger x, BigInteger y) {
    BigInteger a, b, t
    a = x.abs()
    b = y.abs()
    if (b == 0G) {
        return a
    }
    while (b != 0G) {
        t = a % b
        a = b
        b = t
    }
    return a
}

BigInteger pollardRho(BigInteger n) {
    BigInteger c = 2G
    BigInteger x = 1G
    BigInteger y = 1G
    BigInteger d = 1G
    BigInteger savedX = x
    BigInteger count = 0
    println("Try first the Pollard rho algorithm with c = $c")
    while (d == 1G && count*count <= n) {
        x = f(x, c, n)
        if (x == savedX) {
         println("It is cyclic.  x = $x")
         break
        }
        y = g(y, c, n)
        d = gcd((x - y).abs(), n)
        count = count + 1
        // if (count % 5000 == 0) {
        //  println("  count = $count")
        // }
    }

    println("d = $d, count = $count")
    if (d > 1 && d < n) {
        return d
    }

    c = 3
    x = 1
    y = 1
    d = 1
    savedX = x
    count = 0
    println("Try second the Pollard rho algorithm with c = $c")
    while (d == 1 && count*count <= n) {
        x = f(x, c, n)
        if (x == savedX) {
         println("It is cyclic.  x = $x")
         break
        }
        y = g(y, c, n)
        d = gcd(Math.abs(x - y), n)
        count = count + 1
        // if (count % 5000 == 0) {
        //  println("  count = $count")
        // }
    }

    println("d = $d, count = $count")
    if (d > 1 && d < n) {
        return d
    }

    c = 1
    x = 1
    y = 1
    d = 1
    savedX = x
    count = 0
    println("Try third the Pollard rho algorithm with c = $c")
    while (d == 1 && count*count <= n) {
        x = f(x, c, n)
        if (x == savedX) {
         println("It is cyclic.  x = $x")
         break
        }
        y = g(y, c, n)
        d = gcd(Math.abs(x - y), n)
        count = count + 1
        // if (count % 5000 == 0) {
        //  println("  count = $count")
        // }
    }

    println("d = $d, count = $count")

    return d
}

BigInteger n, k, z
n = 9991
if (args.length > 0) {
    n = new BigInteger(args[0])
}

def time1 = new GregorianCalendar().getTime().getTime()

k = pollardRho(n)
z = n/k
// if (n == k*z) {
    println("$n = $k * $z")
// }

def time2 = new GregorianCalendar().getTime().getTime()
print "Elapsed time: ${(time2 - time1)/1000.0} sec\n"

 

크리에이티브 커먼즈 라이선스

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-동일조건변경허락 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.