BOJ 13926 gcd(n, k) = 1

 

이 문제는 이전에 풀었던 BOJ 11689 GCD(n,k) = 1 의 강화버전인 문제이다. 

문제에서 구하고자 하는 바와 기본적인 코드 구조는 거의 동일하다고 볼 수 있다.

단지, 입력되는 숫자의 범위가 기존 10^12에서 10^18로 늘어났단 점만 제외하면 말이다.

 

이렇게 범위가 늘어남으로 인해 많은 것이 바뀌었다.

 

기존에 11689번 문제는 에라토스테네스의 체를 이용해서 소수의 리스트들을 뽑아낸 뒤 (10^6까지의 소수만 뽑으면 된다) 이를 직접 입력된 n에 대입해가면서 인수인지 아닌지를 구분해 낸 뒤, 이를 소인수 리스트에 저장해서 오일러 피 함수를 구하는 방식으로 진행됐다.

 

하지만 수의 범위가 10^18이 되면서, 에라토스테네스의 체를 사용하는 것이 불가능하게 되었다.

에라토스테네스의 체는 10^6을 넘어가는 순간 엄청나게 많은 시간이 소요되는 알고리즘이다.

따라서 최소 10^9까지의 소수를 구해야 하는 위 문제에는 어울리지 않았다.

 

그래서 아래 알고리즘 분류를 봤더니,,

 

예??? 폴라드 로는 뭐고 밀러-라빈 소수 판별법은 뭐죠???
생전 처음보는 알고리즘의 등장해 나는 나의 코딩메이트, 대 황 GPT형님을 부르게 됐다.

 

그렇게  GPT형님에게 물어본 결과, 폴라드 로는 결론적으로 인수분해를 O(n^(1/4)) 에 할 수 있게 해주는 인수정리 알고리즘이고 밀러-라빈 소수 판별법은 10^9 까지 빠른 속도로 소수를 판별할 수 있는 소수 판별법이라는 것을 알았다.

 

그렇게 폴라드 로 알고리즘에 대해 찾아본 결과 구현에 큰 어려움은 없을것으로 예상됐다. 그냥 함수 하나 세워두고 삐까뻔쩍 계산하면 나오는 거라서 딱히 어려워보이진 않았다.

 

 

하지만 진짜 문제는 밀러-라빈 소수 판별법에 있었다. 이 알고리즘은 사실 구현 자체는 매우 쉬웠다. 그냥 소수 리스트 만들어둔다음 공식대로 제곱시키고 모듈러 연산 하고 뭐시기 하고 어떻게 다 대입해서 소수인지 판별하면 되는거니까 크게 상관은 없었다. 

 

 

하지만, 진짜 큰 문제는 이 소수 판별법이 100% 보증된 소수를 판별할 수 없다는 것이다. (이딴게 무슨 소수판별법?)

그나마 밑으로 들어가는 a의 소수 리스트들을 많이 만들어서(많이 만든것이 2~37까지의 소수 리스트들이다. GPT형님 피셜 2^64-1 까지의 숫자까지 거의 정확히 소수를 분간할 수 있다고 한다.) 확률을 높이는 수밖에 없었다. 

 

골치가 아파졌다. 어떻게 하면 정확한 소수를 구할 수 있을까,,, 게다가 내가 짠 코드에서 자꾸 무한 루프가 일어나서 머리가 더 아파졌다. 

 

그러다 문득 생각해냈다. '소수 판별법으로 구한 소수도 다시 소수 판별기에 넣어버리면 되는거 아님?' 

그렇게 나는 기존 코드에서 소인수를 더욱 정확히 판별하기 위해 소인수가 될때까지 해당 숫자를 무한 재귀시켜서 소수를 100% 정확도로 구할 수 있도록 짜게 됐다. 

 

하지만 이 과정에서 자꾸 무한재귀가 일어나거나 엉뚱한 인수가 나오거나, 값이 거지같이 나오거나,,, 

참  많은 시련이 있었지만 꾸역꾸역 해낸 끝에 코드를 짜고 결국 정답을 맞추는 것에 성공했다.

 

 

from math import gcd
from random import randint

def IsPrime(n):
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    a = [2,3,5,7,11,13,17,19,23,29,31,37]
    s = 0
    d = n - 1
    while d % 2 == 0:
        d //= 2
        s += 1
    for ra in a:
        if ra >= n:
            continue
        x = pow(ra, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

def FindPrimeFactor(n):
    def F(x, c): return (x * x + c) % n

    if n % 2 == 0:
        return 2

    while True:
        c = randint(1, 10)
        x, y, d = 2, 2, 1
        while d == 1:
            x = F(x, c)
            y = F(F(y, c), c)
            d = gcd(abs(x - y), n)
        if d == n or d == 1:
            continue
        return d

def EulerRun(n):
    if n == 1:
        return 1
    x = n
    primelst = []

    def recursion(n):
        if n == 1:
            return
        if IsPrime(n):
            primelst.append(n)
            return
        d = FindPrimeFactor(n)
        recursion(d)
        recursion(n // d)

    recursion(n)

    primeset = set(primelst)
    for p in primeset:
        x = x // p * (p - 1)
    return x

N = int(input())
print(EulerRun(N))

 

위 코드에서 IsPrime이 밀러-라빈 소수 판별기 부분, FindPrimeFactor가 폴라드 로 알고리즘 부분, 그리고 EulerRun이 마지막 오릴러 피함수를 구하는 부분이다.