[백준 4149번] 큰 수 소인수분해 - 폴라드로(PollardRho) + 밀러라빈(MillerRabin)

문제 링크: https://www.acmicpc.net/problem/4149

4149번: 큰 수 소인수분해

<참조>

Pollrad's Rho Algorithm

https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm#The_results

https://www.geeksforgeeks.org/pollards-rho-algorithm-prime-factorization

Miller-Rabin

https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Miller%E2%80%93Rabin_test

 

<문제 풀이> 수학, 정수론, 폴라드로(Pollard-Rho), 밀러라빈(Miler-Rabin)

O(N^(1/4))에 소인수 분해를 할 수 있는 폴라드로 알고리즘을 사용하면 되는데 특정 소수에 대해서는 무한루프가 생긴다. 따라서 소수임을 O(k * log^3n) 다항 시간에 판별할 수 있는 밀러라빈 알고리즘을 함께 써야 한다.

밀러라빈을 구현할 때 빠른 거듭제곱 알고리즘 사용

 

<소스코드>

import sys
import random
sys.setrecursionlimit(10 ** 6)
n = int(input())
 
 
def power(x, y, p):
    res = 1
    x = x % p
    while y > 0:
        if y & 1:
            res = (res * x) % p
        y = y >> 1
        x = (x*x) % p
    return res
 
def gcd(a, b):
    if a < b:
        a, b = b, a
    while b != 0:
        r = a % b
        a = b
        b = r
    return a
 
def miller_rabin(n, a):
        r = 0
        d = n-1
        while d % 2 == 0:
            r += 1
            d = d//2
 
        x = power(a, d, n)
        if x == 1 or x == n-1:
            return True
 
        for i in range(r-1):
            x = power(x, 2, n)
            if x == n - 1:
                return True
        return False
 
def is_prime(n):
    alist = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    if n == 1:
        return False
    if n == 2 or n == 3:
        return True
    if n % 2 == 0:
        return False
    for a in alist:
        if n == a:
            return True
        if not miller_rabin(n, a):
            return False
    return True
 
 
def pollardRho(n):
    if is_prime(n):
        return n
    if n == 1:
        return 1
    if n % 2 == 0:
        return 2
    x = random.randrange(2, n)
    y = x
    c = random.randrange(1, n)
    d = 1
    while d == 1:
        x = ((x ** 2 % n) + c + n) % n
        y = ((y ** 2 % n) + c + n) % n
        y = ((y ** 2 % n) + c + n) % n
        d = gcd(abs(x - y), n)
        if d == n:
            return pollardRho(n)
    if is_prime(d):
        return d
    else:
        return pollardRho(d)
 
list = []
while n > 1:
    divisor = pollardRho(n)
    list.append(divisor)
    n = n // divisor
 
list.sort()
 
for i in list:
    print(i)