Index Calculus - 조금은 특이한 이산로그 계산 방식

이산로그 (Dicrete Logarithm) 와 관련된 문제를 푸는 도중 매우 난감한 문제를 만났다.

 

https://www.acmicpc.net/problem/21094

21094번: Discrete Logarithm is a Joke

 

위 문제의 확장형인 아래의 문제이다.

 

https://level.goorm.io/exam/117941/discrete-logarithm-is-not-a-joke/quiz/1

구름LEVEL

https://www.acmicpc.net/problem/21864

21864번: 이산로그가 장난이냐?

(같은 문제이지만 BOJ 버전이다.)

 

요약하자면, $M = 10^{18}+31$, 원시근 $g=42$, $a_0 = 960002411612632915$ 으로 고정된 값이고, 

$f(x) = log_g{x}\ mod\ M$ , $a_{i+1}=f(a_i)$ 라고 정의된 점화식에서 $a_n$ 값을 구하는 문제이다.

 

BOJ 에 올라온 문제는 이산로그 문제인척 하는 단순한 유머성 문제이지만

sait2000 님이 만든 문제는 진짜로 이산로그를 풀어야지만 풀 수 있는 문제이다.ㅠㅠ

 

덤 : 21094번 문제 초고속 풀이.

$f^{-1}(x) = g^x$ 이기 때문에 $a_i = f^{-1}(a_{i+1}) = g^{a_{i+1}}$ 을 계산하면 되므로,

$a_0$ 부터 시작해서 $i$ 값을 증가해가면서가 아니라,

적당히 큰 $a_{1000000}$ 부터 시작해서 $i$ 값을 감소시켜가면서 값을 계산하면 되는 문제이다.

그리고 예제에서 $a_{1000000} = 300$ 이라고 알려주었다. 풀이 끝.

 

sait2000 님이 만든 문제를 풀기 위해서는 수단과 방법을 가리지 않고 $a_{2000000}$ 값을 구해낸 뒤, 위 풀이와 똑같은 과정을 거치면 된다.

 

 

일반적으로 $mod\ M$ 에 대한 이산로그 알고리즘들은 이 블로그 등에서 설명되어있듯이,

1. $\mathcal{O}(\sqrt{M})$ 에 계산

2. $M-1$ 을 이루는 소인수들이 적당히 작은 값들로 이루어짐

인 상황에서 답을 잘 구한다.

하지만 위 문제는 $M=10^{18} + 31$ 로, $\sqrt{M}$ 값도 매우 크고 $M-1=2*5*(10^{17}+3)$ 인데 $10^{17}+3$ 도 소수라서 $M-1$ 의 소인수도 매우 크다. 게다가 이런 연산을 적어도 $10^6$ 번 반복해서 계산해야한다. 꿈도 희망도 없다.

 

---

이 문제를 푼 사람들로부터 풀이를 들었는데, Index Calculus 라는 방법을 사용해서 해결하면 합리적인 시간 안에 답이 나온다고 알려주었다. 그래도 몇시간 단위이긴 하지만...

 

https://en.wikipedia.org/wiki/Index_calculus_algorithm 

Index calculus algorithm - Wikipedia

관련 위키피디아 링크이다. 사실은 위키링크 하나 보고 스스로 깨우치고 공부한 내용이나 다름없다. ㅠㅠ

 

Index Calculus 가 무슨 짓을 하는지 간략하게 로직을 나눠보자.

원시근 $g$ 가 알려진 상태일 때, $log_g{h}$ 를 계산하기 위해 다음과 같은 작업을 한다.

 

  step 1. $P$ 번째 소수까지 $log_g{p_i}$ 값을 계산한다.

  step 2. $g^s*h\ mod\ M$ 결과값을 소인수분해해서 $P$번째 소수보다 큰 소인수가 없기를 기대한다. 그 결과값을 $g^s*h = p_0^{c_{0}} * p_1^{c_{1}} * ... * p_{P-1}^{c_{P-1}}$ 라고 표현해보자.

  step 3. 양변에 로그를 취한다. $log_g{g^s*h} = s + log_g{h} = \sum{c_i * log_g{p_i}}$ 우리는 $s$ 와 $c_i$, $log_g{p_i}$ 을 모두 알고있으므로, $log_g{h}$ 값을 구할 수 있다.

 

 

step 3. 는 단순 계산이고, step 2. 는 천지신명이 지정해주는 영역이기 때문에 평범한 인간따위가 컨트롤할 수 있는 영역이 아니다. (수학 잘하는 사람들은 컨트롤 할 수 있으려나...?)

문제를 미리 풀고 온 자로써 할 수 있는 말은.. step 2. 를 만족하는 상황이 의외로 잘 발생한다는 것 정도..?

 

 

따라서 step 1. 이 가장 중요하다는 것을 알 수 있다.

심지어 step 1. 을 한번 계산해두면 $h$ 값이 여러 번 들어오는 상황애도 재활용할 수 있다는 것을 알 수 있다!

이산로그를 $10^6$ 번 계산해야하는 우리에게 아주 알맞은 알고리즘이다.

 

 

그럼 이제 $log_g{p_i}$ 를 구해보자.

시작은 step 2. 의 방식과 비슷하다. 적당한 $k$ 에 대해, $g^k = p_0^{c_{0}} * p_1^{c_{1}} * ... * p_{P-1}^{c_{P-1}}$ 를 만족하기를 기대한다. (사실은 어떠한 $k$가 되어도 상관없기 때문에, 그냥 1부터 천천히 증가시켰다.)

이 결과를 길이 $P+1$ 인 벡터로 표현해보자.

 

$v = \left( c_0, c_1, .... ,c_{P-1}, k\right)$

 

이 벡터의 각각 요소들은 $\mathbb{Z}_{M-1}$ 위에서 움직인다. ($log_g$ 값이고, $M$ 의 배수가 아닌 어떤 값이든 $x^{M-1} = x^{0} = 1$ 이 되기 때문이다.)

이런 벡터들을 여러 개 모아서, 그 중 선형독립 적인 벡터 $P$ 개를 찾으면, 가우스 소거법 을 사용해서 실제 답을 구할 수 있다.

 

실제로 구현한다고 생각하면... 선형독립적인 $R$ 개의 벡터 집합을 미리 가지고있고, ($R$ 값을 보통 rank라고 부른다.)

새로 구한 벡터 $v$ 가 위 $R$ 개의 벡터간의 산술연산으로 만들어지는 벡터인지 테스트한다.

만약 $R$ 개의 벡터와 독립적인 새로운 벡터인게 확인되면, $v$ 를 벡터 집합에 추가하고 rank 값 $R$ 을 1 증가시킨다.

만약 테스트를 통과히자 않았다면, $v$ 를 버리고, 새로운 벡터가 들어올 때 까지 기다린다.

이 과정을 rank 가 $P$ 가 될 때까지 반복하면 된다.

 

실제 구현할 때 각 요소들이 $\mathbb{Z}_{M-1}$ 위에서 움직이기 때문에 나눗셈을 계산할 수 없는 경우가 발생할 수 있다. 이런 상황이 발생하지 않도록 주의하면서 코드를 작성해야한다.

이 문제에서는 $M-1$ 의 소인수의 가지수가 적어서 역원 계산시 곤란해지는 경우가 오히려 적어지는 소소한 특징이 있다.

 

여기서 덤으로 $log_g{-1} = \frac{M-1}{2}$ 임을 활용해서, $g^{k+\frac{M-1}{2}}$ 에 대해 한번 더 테스트 할 수도 있다.

---

 

아래부터는 실제로 python3 코드로 작성한 내용들이다.

수들이 $10^{18}$ 범위에서 이것저것 곱셈 연산을 많이 하므로, C++ 보다는 python 을 활용하기로 결정했다.

import time
start_time = time.time()
 
g = 42
m = 1000000000000000031
PAM = 10000
p = []
pa = [0 for _ in range(PAM)]
pl = 0
for i in range(2,PAM):
    if pa[i] == 0:
        p.append(i)
        pl += 1
        for j in range(i, (PAM-1)//i+1):
            pa[i*j] = 1
print("prime count : "+str(pl))
f = [0 for _ in range(pl)]
 
def egcd(a1, a2):
  x1, x2 = 1, 0
  y1, y2 = 0, 1
  while a2:
    q = a1 // a2
    a1, a2 = a2, a1 - q * a2
    x1, x2 = x2, x1 - q * x2
    y1, y2 = y2, y1 - q * y2
  return (x1, y1, a1)
 
mi = [i for i in range(pl+2)]
mat = [[0 for _ in range(pl+2)] for __ in range(pl+2)]
rank = 0
def increase_matrix(k):
    global rank, mat, mi
    for i in range(pl):
        mat[rank][i] = f[mi[i]]
    mat[rank][pl] = k
    for i in range(rank):
        if mat[rank][i] > 0:
            scale = mat[rank][i]
            for j in range(pl+1):
                mat[rank][j] = (mat[rank][j] - scale * mat[i][j]) % (m-1)
    ranki = -1
    key = -1
    rev = -1
    for j in range(rank, pl):
        if mat[rank][j] > 0:
            key = mat[rank][j]
            rev, y, g = egcd(key, m-1)
            if g == 1:
                ranki = j
                break
    if ranki < 0: return False
    rev %= (m-1)
    if ranki != rank:
        for i in range(rank+1):
            mat[i][rank], mat[i][ranki] = mat[i][ranki], mat[i][rank]
        mi[rank], mi[ranki] = mi[ranki], mi[rank]
    
    for j in range(rank, pl+1):
        mat[rank][j] = mat[rank][j] * rev % (m-1)
    for i in range(rank):
        if mat[i][rank] > 0:
            scale = mat[i][rank]
            for j in range(rank, pl+1):
                mat[i][j] = (mat[i][j] -  scale * mat[rank][j])%(m-1)
    rank += 1
    print("rank : "+str(rank)+" / "+str(k))
#    for j in range(pl+1):
#        print(mi[j], end=' ')
#    print()
#    for i in range(rank):
#        for j in range(pl+1):
#            print(mat[i][j], end=' ')
#        print()
    return rank == pl
 
tg = 1
for k in range(1, 1000000000):
    tg = (tg * g) % m
    tgg = tg
    for i in range(pl):
        while tgg % p[i] == 0:
            f[i] += 1
            tgg //= p[i]
    if tgg == 1:
#        for i in range(pl):
#            print(f[i], end=' ')
#        print(" : " + str(tg))
        if increase_matrix(k):
            break
    for i in range(pl):
        f[i] = 0
 
    tgg = m - tg
    for i in range(pl):
        while tgg % p[i] == 0:
            f[i] += 1
            tgg //= p[i]
    if tgg == 1:
        if increase_matrix(k + m//2):
            break
    for i in range(pl):
        f[i] = 0
 
print("time :" + str(time.time() - start_time) + " s") 
 
file = open("result.txt","a")
file.write('%d\n' % pl)
for i in range(pl):
    file.write('%d ' % p[mi[i]])
    file.write('%d\n' % mat[i][pl])
file.close()

 

위 코드는 $PAM\ (=10000)$ 이하인 모든 소수에 대해서 $log_g{p}$ 값을 계산하는 프로그램이다.

전체적인 시간복잡도는 $\mathcal{O}(P^2*(K+P))$ 정도 걸린다. ($P$ 는 소수의 개수, $K$ 는 검사한 수의 개수)

(확장 유클리드 함수 egcd 는 직접 짜기 귀찮아서 https://www.secmem.org/blog/2019/05/17/%EC%9D%B4%EC%82%B0-%EB%A1%9C%EA%B7%B8/ 여기서 복붙해왔다. 감사합니다.)

 

계산 과정에서 $\mathbb{Z}_{M-1}$ 의 역원이니 뭐니 귀찮은게 많아질 것 같아서...

애초에 행렬에 백터를 추가할 때 $mat[0\sim rank][0\sim rank]$ 는 전부 단위행렬이 되도록 열의 위치를 적절히 조절해버렸다.

적절히 조절당한 열의 위치를 찾아주는 변수는 코드상에서 $mi$ 라는 이름으로 있다.

 

나의 구데기 노트북 에서 pypy3 로 계산한 결과 $K=1700700$ 이 나왔고... 계산은 약 268초가 걸렸다.

 

이렇게 하면 result.txt 에 소수 $p$ 와, $log_g{p}$ 값이 한줄에 하나씩 총 $P$ 줄에 걸쳐서 계산이 되있을 것이다.

 

해당 정보를 활용해서 진짜로 이산로그값을 구해보자.

import time
start_time = time.time()
 
inp = open("result.txt","r")
g = 42
m = 1000000000000000031
pl = int(inp.readline())
 
p = [0 for _ in range(pl)]
q = [0 for _ in range(pl)]
for i in range(pl):
    p[i], q[i] = map(int, inp.readline().split())
    if p[i] != pow(g, q[i], m):
        print("wrong at line "+str(i)+" "+str(p[i])+" "+str(q[i]))
inp.close()
 
n = 1000000
an = 300
#n = 0
#an = 960002411612632915
f = [0 for _ in range(pl)]
def cal(h):
    tg = h
    k = 0
    while True:
        tgg = tg
        fl = pl
        for i in range(pl):
            while tgg % p[i] == 0:
                tgg //= p[i]
                f[i] += 1
            if tgg == 1:
                fl = i+1
                break
        if tgg == 1:
            res = 0
            for i in range(fl):
                res = (res + f[i] * q[i]) % (m-1)
                f[i] = 0
            return (res - k) % (m-1)
        for i in range(fl):
            f[i] = 0
        tg = tg * g % m
        k += 1
 
 
file = open("log.txt","a")
while n < 2000000:
    tn = cal(an)
    print(n)
    n += 1
    if n%10000 == 0:
        file.write('%d ' % n)
        file.write('%d\n' % tn)
        print(n, end=' ')
        print(tn)
    an = tn
file.close()
print("time :" + str(time.time() - start_time) + " s") 

 

우선, 이산로그를 구하고자 하는 값 $h$ 에다가 $g^s$ 를 곱한 결과값의 소인수가 앞서 계산한 $P$ 개의 소수들로만 이루어져있기를 기도한다.

그 값 $g^s*h$ 를 log 한 결과값은 앞서 계산한 $log_g{p}$ 테이블을 활용해서 값을 구할 수 있다.

즉, $log_g{(g^s * h)} = s + log_g{h}$ 를 $\mathcal{O}(P)$ 만에 계산할 수 있다.

따라서 $log_g{h}$  는 위에서 계산한 결과값에서 $s$ 만 빼주면 된다.

 

$log_g{-1} = \frac{M-1}{2}$ 인 점을 이용해서 테스트를 한번 더 할 수 있긴 하지만, 귀찮아서 따로 하진 않았다ㅋㅋㅋㅋ

 

전체 복잡도는 전처리에 시간 $\mathcal{O}(P^2(K+P))$, 공간 $\mathcal{O}(P^2)$. 실제 계산에는 공간 $\mathcal{O}(P)$, 시간은 쿼리당 $\mathcal{O}(sP)$ 만큼의 시간이 소요된다.

 

여기서 우리는 $P$ 값만 조절할 수 있고, $K$ 나 $s$ 값은 예측할 수 없는 값이라는 특징이 있다.

대신, $P$ 값이 너무 작으면 $g^K$ 값이 $P$개의 소인수로만 이루어져있을 확률이 점점 작아지므로, $P$ 가 어느정도 커야지만 그나마 합리적으로 동작하지 않을까? 라고 예측할 수 있는 정도이다.

 

나는 이 문제를 풀 때, 이론 학습에 약 1시간, 코딩에 1시간 (파이썬이 익숙하지 않아서 ㅠㅠ), 전처리에 5분, $10^6$개의 쿼리 계산에 10시간 이하가 소요되었다. (전날 2시에 돌려놓고 자고 일어나니 계산이 끝나있었다. 시간 측정은 못함 ㅠㅠ)