[정수론](-)[펠 방정식의 해를 코드로 구하기]

펠 방정식 정리

• 항상 양의 정수해를 가짐
• 해 (x1, y1)이 가장 작은 정수해라면?
  모든 해 (xk, yk)를 $x_k+y_k\sqrt{D} \ = (x_1+y_1\sqrt{D})^k$ 로 구할 수 있다.

펠 방정식 정리에 의하면, 제일 작은 해만 구하면 나머지 해를 알 수 있다고 한다.
근데 손으로 구하기엔 너무 노가다임...
-> 컴공 식으로 해결해보자

 

$(x_1+y_1\sqrt{D})^k$에서 어떻게 계수만 뽑아낼 수 있을까?
바로 이항정리를 응용하면 된다.

$$(x_1+y_1\sqrt{D})^k = \sum_{k=0}^{n}{{n \choose k}\space x_1^{n-k}\space {(\sqrt{D}y_1)}^{k}}$$
로 표현할 수 있고, k에 따라 $\sqrt{D}$가 나오는지 여부가 갈린다
    k = 홀수: $\sqrt{D}$ 나옴
    k = 짝수: $\sqrt{D}$ 안나옴
즉,
$$x_k = \sum_{k=짝수}^{n}{{n \choose k}\space x_1^{n-k}\space {(\sqrt{D}y_1)}^{k}}$$
홀수일땐 $\sqrt{D}$ 를 처리하기가 애매하므로 걍 $\sqrt{D}$를 나눠주겠다.
$$y_k = \sum_{k=홀수}^{n}{{n \choose k}\space x_1^{n-k}\space {y_1}^{k}\space {\sqrt{D}^{k-1}}}$$

def factorial(n):
    if (n == 1 or n == 0):
        return 1
    else:
        return factorial(n-1)*n

def binomial_even(n, k):
    bin = factorial(n)/(factorial(n-k)*factorial(k))
    return bin*(X1**(n-k))*(Y1**k)*(D**(k/2))

def binomial_odd(n, k):
    bin = factorial(n)/(factorial(n-k)*factorial(k))
    return bin*(X1**(n-k))*(Y1**k)*(D**((k-1)/2))

def find_small_x(D):
    sx=1
    sy=1
    while(sx**2-D*sy**2 != 1):
        sy=1
        for i in range(1, sx+1):
            if (sx**2-D*sy**2 == 1):
                return sx, sy
            sy += 1
        sx+=1
    return sx, sy

## user parameter ##
D = 5
k = 4

## calc ##
X1, Y1 = find_small_x(D)
xk = 0
yk = 0
for i in range(0, k+1):
    if (i%2 == 0):
        xk += binomial_even(k, i)
    else:
        yk += binomial_odd(k, i)

print("찾은 가장 작은 해 = (%d, %d)"%(find_small_x(D)))
print("x^2 - %dy^2 = 1의 %d번째 해 (%d, %d)" %(D, k, xk, yk))
print("검산: %d^2 - %d*%d^2 = %d"%(xk, D, yk, xk**2-D*yk**2))

예) 펠 방정식 x^2 - 5y^2 = 1의 양의 정수해를 4개 구하시오.
    코드에서 보면, 차근차근 대입해서 (9, 4)를 구하고
    k값을 1, 2, 3, 4 이렇게 넣어주면 된다.

딸-깍