오일러 프로젝트 66

오일러 프로젝트 66번 문제는 펠 방정식의 최소해에 관한 문제이다. 나이브하게 접근했다가는 결코 풀어낼 수 없을 수준으로 시간이 많이 걸린다. 그리고 그나마 빠른 해법 역시 구현하기 위해 여러 가지 지식과 스킬이 동원된다.  접근 방법을 알고 있더라도 구현이 만만치 않은, 1번부터 현재까지는 최고 난이도의 문제라 할 수 있다. 통상 오일러 프로젝트의 문제들은 순서에 맞게 풀어나갈 필요는 없지만, 이 문제 만큼은 57번, 64번, 65번을 풀 때의 지식이 그대로 요구된다.

다음과 같은 2차 디오판토스 방정식이 있습니다.

x2 – Dy2 = 1

(역주: 디오판토스 방정식은 정수해만 허용하는 부정 다항 방정식입니다. 그 중에서도 여기 나온 형태는 펠의 방정식이라고 하는데, x/y로 D의 제곱근의 근사값을 구하는데 이용됩니다.)

예를 들어 D=13일 때, x를 최소화하는 해는 6492 – 13×1802 = 1 이 됩니다.  D가 제곱수일 때는 위 식을 만족하는 자연수 해는 없다고 볼 수 있습니다.

D = {2, 3, 5, 6, 7}에 대해서 x를 최소화하는 자연수 해를 찾아보면, 아래와 같은 결과를 얻게 됩니다.

32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1

위에서 보듯이 D ≤ 7 인 경우에 대해 x를 최소화하는 해를 구하면, x의 값이 가장 큰 것은 D=5 일 때 입니다. D ≤ 1000인 경우에 x를 최소화하는 해를 구하면, 가장 큰 x의 값을 갖는 D는 얼마입니까?

접근

매우 어려운 문제이다. 그리고 디오판토스 방정식에 대해서는 교과과정에서는 뭘 심도있게 다루거나 하지는 않고 (그저 엄청 쉬운 예제 몇 개 던져주는게 전부였던 듯) 하기 때문에 관련해서는 많은 구글링이 필요했다. 실제로 문제를 풀어서 답을 보는데까지 한 달 넘게 고민했던 것 같다.

문제의 힌트는 역주에 있다. 문제에 주어진 패턴의 방정식을 펠의 방정식이라 하는데, (1, 0)과 같은 사소한 해를 시작으로 무한히 많은 해를 가지고 있음이 증명되어 있다. 펠의 방정식은 기하학적으로는 쌍곡선을 표현한다.

펠의 방정식의 해의 근 x, y에서 x/y는 D의 제곱근을 유리근사한 것에 해당하며, 이는즉 D의 제곱근을 연분수로 확장하고 이를 전개한 식으로부터 이 방정식의 해를 찾을 수 있다는 의미이다. (물론 모든 연분수 전개가 펠 방정식의 해에 해당하지는 않는다.) 어떤 펠 방정식의 최소해는 엄청 클 수 있다. 유명한 예로 x2 – 313y2 = 1이 있는데 이 해는 (32188120829134849, 1819380158564160)이다. (이게 x가 최소가 되는 해이다.)

영문 위키의 펠 방정식 항목에는 1에서 128까지의 펠 방정식의 최소해가 소개되어 있다.

이 문제를 푸는데 있어서 필요한 정보는 여기까지이고, 그렇다면 어떤 과정을 통해서 풀어야 할까?

  1.  임의의 자연수 D에 대해서 제곱근을 연분수로 표현하기 위한 순환마디를 구하는 함수가 요구된다.
  2. 1에서 구한 정수부와 순환마디를 이용해서 n차례 확장한 연분수를 전개하는 함수가 요구된다.
  3. 1과 2를 통해서 임의의 D에 대해 1, 2, 3…. 차로 확장하면서 계산된 분자와 분모를 위 식에 대입하여 방정식의 해인지 확인한다. 만약 해가 된다면, 그 때의 x가 해당 방정식의 최소해 (단, 사소한 해를 제외하고)가 된다.

이 과정을 거쳐서 D에 대한 펠 방정식의 최소해를 구하고, 이 값의 최대값과 그 때의 D를 구하면 된다.

참고로 D가 완전제곱수인 경우에는 (1, 0)의 사소한 해만이 유일한 해가 된다.

Python 구현

위에서 언급한 세가지 함수를 각각 작성해보자

제곱근의 연분수 표현식 생성하기

64번 문제로부터 살짝 변형된 함수를 하나 만든다.

def squared_info(n: int) -> [int]:
  m = int(n**0.5)
  x = m, -m, 1
  cache, result = set(), [m]

  ## 점화식을 이용해서 진행하기
  def process(a1, b1, c1) -> (int, int, int):
    c2 = (n - b1**2) // c
    if c2 is 0: ## N이 완전제곱수이다.
      return None
    a2 = (m - b1) // c2
    b2 = -b1 - (a2 * c2)
    return a2, b2, c2

  while True:
    x = process(*x)
    if x is None:
      return None
    if x not in cache:
      cache.add(x)
      result.append(x[0])
    
## 테스트
print(squared_info(23)) ## [4,1,3,1,8]

이 함수의 리턴에서 첫번째 원소는 정수부, 두번째부터는 순환마디를 표현한다. 이를 통해서 제곱근의 연분수 전개를 계산하는 함수를 작성해보자. 유리수 표현을 위해서 오늘도 fractions 모듈이 수고해주시겠다.

연분수 전개

from fractions import Fraction

def calc(ns:[int], k=1) -> Fraction:
  """ ns: squared_info의 리턴, k: 순번 
  k = 0 이면 정수부만 출력
  k =1, 2, 3 ... 이면 변환과정을 거친다.
  """
  if ns is None: return Fraction(1, 1)

  a, *xs = ns
  if k is 0: return Fraction(a, 1)

  ## i 번째 순서가 가리키는 순환마디 내 위치
  pos = (i - 1) / len(xs)
  f = Fraction(xs[pos], 1)
  while i > 0:
    f = xs[pos] + 1 / f
    i -= 1
    pos = pos - 1 if pos > 0 else len(xs) - 1
  ## 최종적으로 정수부를 더하는 과정을 추가한다. 
  f = a + 1 / f
  return f

이 함수를 이용해서 57번 문제에서 소개된 2의 제곱근의 연분수 전개꼴에 대해서 검증해보면 좋다.

k = squared_info(2)
for i in range(1, 11):
  print(calc(k, i))

## 3/2
## 7/5
## 17/12
## 41/29
## 99/70
## 239/169
## 577/408
## 1393/985
## 3363/2378
## 8119/5741

최종 풀이

준비물이 모두 갖춰졌으니 다음과 같이 풀어보자.

def main():
  result = (0, 0)
  for d in range(2, 1001):
    s = squared_info(d)
    if s in None: continue
    i = 1
    while True:
      z = calc(s, i)
      x, y = z.numerator, z.denominator
      if x**^ - d*y**2 == 1:
        ## 해를 찾음
        if x > result[0]:
           result = (x, d)
        break
    print(result[1])

약 2.8초 내외의 시간이 걸렸다.