Wireframe

오일러 프로젝트 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초 내외의 시간이 걸렸다.

Swift 풀이

이 문제의 Swift 풀이에서 어려운 점이라면 역시, 큰 수 계산에 관한 것이다. 이전까지 큰수의 덧셈과 곱셈 그리고 큰수와 정수의 곱셈과 덧셈을 구현해서 썼는데, 이번에는 “최대값”을 찾아야 하기 때문에 큰 수의 대소 및 동등 비교를 추가로 구현해야 한다.
하지만 이 역시 우리가 십진수에서 대소 비교하는 것과 동일한 방식으로 구현해주면 된다.

  1. 자리수가 많으면 많은쪽이 크다.
  2. 자리수가 같으면 앞쪽부터 숫자(계수)가 큰쪽이 크다.

동등 비교 역시, 자리수가 같고, 모든 항의 계수가 같다라고 비교하면 된다. 또한 BigNumber의 각 항(노드)들은 생성 원리에 의해서 작은 차수부터 생성되기 때문에 비교적 간단하게 구현할 수 있다.

extension BigNumber {
  static func == (left: BigNumber, right: BigNumber) -> Bool {
    if left.nodes.count != right.nodes.count { return false }
    for (x, y) in zip(left.nodes, right.nodes) {
      if x.value != y.value { return false }
    }
    return true
  }
  static func > (left: BigNumber, right: BigNumber) -> Bool {
    if left.nodes.count != right.nodes.count { return left.nodes.count > right.nodes.count }
    for i in stride(from: left.nodes.count - 1, through: 0, by: -1) {
      if left.nodes[i].value > right.nodes[i].value { return true }
    }
    return false
  }
}

이 이후에는 이전 문제에서 만든 BigFraction 을 그대로 쓰면 된다.

자연수의 제곱근을 연분수로 표현하는 함수

파이썬과 대동소이한 구현.

func squaredInfo(of n: Int) -> [Int]? {
  let m = Int(sqrt(Double(n)))
  var (a, b, c) = (m, -m, 1)
  var cache = Set<String>()
  var result = [Int]()
  func next(_ a1:Int, _ b1: Int, _ c1: Int) -> (Int, Int, Int)? {
    let c2 = (n - b1 * b1) / c1
    if c2 == 0 { return nil }
    let a2 = (m - b1) / c2
    let b2 = -b1 - (a2 * c2)
    return (a2, b2, c2)
  }
  
  while case let s = "\(a),\(b),\(c)", !cache.contains(s) {
    if let (x, y, z) = next(a, b, c) {
      cache.insert(s)
      result.append(a)
      (a, b, c) = (x, y, z)
    } else { return nil }
  }
  return result
}

연분수 전개를 계산하는 함수

func calc(_ x: [Int], _ i: Int = 1) -> BigFraction {
  var xs = x
  let a = xs.removeFirst()
  var i = i
  if i == 0 { return BigFraction(a) }
  var pos = (i - 1) % xs.count
  var f = BigFraction(xs[pos])
  while i > 1 {
    pos = (pos > 0) ? pos - 1 : xs.count - 1
    f = xs[pos] + f.reversed()
    i -= 1
  }
  f = a + f.reversed()
  return f
}

역시 파이썬과 같은 방식으로 테스트 해볼 수 있다.  테스트는 각자가 해보고… 드디어 본 문제 풀이.

do {
  var result: (BigNumber, Int) = (BigNumber(integer: 0), 0)
  for d in 2...1000 {
    guard let s = squaredInfo(of: d) else { continue }
    var i = 1
    while true {
      let g = calc(s, i)
      let (x, y) = (g.n, g.d)
      if x*x == (y*y*d) + BigNumber(integer:1) {
        if x > result.0 { result = (x, d) }
        break
      }
      i += 1
    }
  }
  print(result)
}
Exit mobile version