연분수(continued fraction)를 이용한 선형 디오판토스 방정식의 해를 구하는 방법

1. 연분수 표현에서 a,b의 최대공약수를 바로 구하는 방법

 

$\frac{A}{B}$의 연분수 표현을 구하고자 한다면, A/B가 기약분수 형태로 들어가줘야한다.

 

기약분수로 바꿀려면 A,B의 최대공약수가 g = gcd(A,B)일때, $p_{k} = A/g, q_{k} = B/g$이고 

 

그래서 $\frac{A}{B} = \frac{p_{k}}{q_{k}}$이다.

 

따라서, A/B가 기약분수가 아닐때, 연분수 표현 p,q를 구한 다음 A를 p[-1]로 나눠주면 A,B의 최대공약수가 된다.

 

#A/B가 기약분수가 아닐때, A,B의 최대공약수 구하기
def continued_fraction(p,q):
    
    f = []

    while q != 0:
        
        f.append(p//q)
        p,q = q,p%q
    
    return f

def convergent(f):
    
    p = [0,1]
    q = [1,0]

    for a in f:
        
        p.append(a*p[-1]+p[-2])
        q.append(a*q[-1]+q[-2])
    
    return p,q

#A/B = p_k/q_k이므로, A = p_k*g, B = q_k*g이다.
#g = A//p[-1] = B//q[-1]

a = 12
b = 8

p,q = convergent(continued_fraction(a,b))

print(a//p[-1])
print(b//q[-1])

 

 

2. ax+by = c의 해

 

정수 a,b,c에 대하여, ax+by = c의 정수해를 구하는 방법으로 가장 유명한 것은 확장 유클리드 알고리즘이다.

 

https://deepdata.tistory.com/576

확장된 유클리드 알고리즘(extended euclidean algorithm) 구현해보면서 익히기

 

하지만 연분수(continued fraction)를 이용해서 특수해를 쉽고 직관적으로 구할 수 있다.

 

$\frac{a}{b} = \frac{p_{k}}{q_{k}} = [a_{0}; a_{1}, a_{2}, ... , a_{k}]$로 나타낸다고 하자.

 

$$\frac{p_{k}}{q_{k}} = \frac{a_{k}p_{k-1} + p_{k-2}}{a_{k}q_{k-1} + q_{k-2}}$$가 성립함을 알고 있다.

 

https://deepdata.tistory.com/1003

컴퓨터가 실수를 연분수 표현(continued fraction)으로 나타내는 방법

 

$$p_{k} = a_{k}p_{k-1} + p_{k-2} , q_{k} = a_{k}q_{k-1} + q_{k-2}$$에서 각각 $q_{k-1}, p_{k-1}$을 곱하면..

 

$$p_{k}q_{k-1} = a_{k}q_{k-1}p_{k-1} + q_{k-1}p_{k-2}, p_{k-1}q_{k} = a_{k}q_{k-1}p_{k-1} + q_{k-2}p_{k-1}$$

 

두 식을 빼주면..

 

$$p_{k}q_{k-1} - p_{k-1}q_{k} = q_{k-1}p_{k-2} - q_{k-2}p_{k-1}$$

 

식을 자세히 관찰하면 $$A_{k} = -A_{k-1}$$임을 알 수 있다.

 

여기서 $A_{k} = p_{k}q_{k-1} - p_{k-1}q_{k}$이다.

 

점화식을 정리해보면...

 

$$A_{k} = -A_{k-1} = (-1)^{2}A_{k-2} = ... = (-1)^{k}A_{0} = (-1)^{k} * (q_{-1}p_{-2} - q_{-2}p_{-1}) = (-1)^{k+1}$$

 

여기서 $p_{-1} = 1, p_{-2} = 0, q_{-1} = 0,q_{-2} = 1$ 이기 때문이다.

 

A,B의 최대공약수가 g = gcd(A,B)일때, 

 

$p_{k}g = A, q_{k}g = B$이다. $\frac{A}{B} = \frac{p_{k}}{q_{k}}$로 A/B의 기약분수 형태이기 때문이다.

 

따라서, 양변에 g를 곱하면,

$$Aq_{k-1} - Bp_{k-1} = (-1)^{k+1}g$$

 

양변을 $(-1)^{k+1}g$로 나누고 C를 곱하면...

 

$$A*(-1)^{k+1}Cq_{k-1}/g - B*(-1)^{k+1}Cp_{k-1}/g = C$$이다.

 

따라서, Ax+By = C의 특수해는 C가 g로 나누어 떨어지면 존재하고..

 

$$x = (-1)^{k+1}Cq_{k-1}/g, y = (-1) *  (-1)^{k+1}Cp_{k-1}/g$$이다.

 

def continued_fraction(p,q):
    
    f = []

    while q != 0:
        
        f.append(p//q)
        p,q = q,p%q
    
    return f

def convergent(f):
    
    p = [0,1]
    q = [1,0]

    for a in f:
        
        p.append(a*p[-1]+p[-2])
        q.append(a*q[-1]+q[-2])
    
    return p,q

#ax+by = c의 특수해를 구하는 함수
def diophantine(a,b,c):
    
    #a/b의 연분수 표현 p,q를 구한다
    p,q = convergent(continued_fraction(a,b))
    
    #a = pk * g, b = qk * g이므로, g = a//pk
    g = a // p[-1]
    
    #c가 g = gcd(a,b)로 나누어 떨어지지 않으면 정수해가 존재하지 않는다
    if c % g == 0:
        
        c //= g
    
    else:
        
        return -1,-1,-1
    
    #(-1)^{k+1}의 부호를 정한다
    #k+1 = len(p)이므로, len(p)가 홀수이면 -1
    #len(p)가 짝수이면 1
    t = -1 if len(p) % 2 else 1
    
    #x = (-1)^(k+1)*C*q_(k-1)/g, y = -1*(-1)^(k+1)*C*p_(k-1)/g    
    #여기서 p[-1] = p_k이고 p[-2] = p_k-1
    x = t*c*q[-2]
    y = -t*c*p[-2]
    
    return g,x,y #a,b의 최대공약수, 특수해 x,y

 

연습문제는 다음에..