오일러 프로젝트 64

오일러 프로젝트 64 번 문제는 자연수의 제곱근을 연분수로 전개할 때, 반복되는 패턴을 찾는 문제이다. 이는 자연수의 제곱근을 연분수로 전개해서 풀어 쓰는 과정을 추적하여, 연분수식 모양내에서의 각 위치의 항의 계수간의 점화식을 찾는 것으로 출발할 수 있다.

모든 제곱근은 아래와 같이 연분수로 나타낼 수 있는데, 이 때 반복되는 부분이 나타납니다.

\sqrt{N} = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + ...}}}

예를 들어서 \sqrt{23}을 풀어 보면,

\sqrt{23} = 4 + (\sqrt{23} - 4) = 4 + \frac{1}{(\frac{1}{\sqrt{23} - 4)}} = 4 + \frac{1}{\frac{\sqrt{23}+4}{7}} = 4 + \frac{1}{1 + \frac{\sqrt{23}-3}{7}}

같은 식으로 계속하면 아래와 같은 모양이 됩니다.

\sqrt{23} = 4 + \frac{1}{1 + \frac{1}{3 +\frac{1}{1 + \frac{1}{8 + ...}}}}

이 과정을 상세히 보면 다음과 같습니다.

위에서 보듯이 4라는 정수부 다음에 1,3,1,8 이라는 숫자가 무한히 반복되는데, 이것을 \sqrt{23} = [4;(1,3,1,8)]과 같이 표시하기로 합니다. 이런 식으로 해서 무리수인 제곱근들을 연분수로 나타내면 다음과 같이 됩니다.

√2=[1;(2)], 주기=1
√3=[1;(1,2)], 주기=2
√5=[2;(4)], 주기=1
√6=[2;(2,4)], 주기=2
√7=[2;(1,1,1,4)], 주기=4
√8=[2;(1,4)], 주기=2
√10=[3;(6)], 주기=1
√11=[3;(3,6)], 주기=2
√12=[3;(2,6)], 주기=2
√13=[3;(1,1,1,1,6)], 주기=5

반복 주기가 홀수인 경우는 N <= 13일 때 모두 4번 이음을 볼 수 있습니다. 그러면 N <= 10000일 때 반복 주기가 홀수인 경우는 모두 몇 번이나 있습니까?

접근

연분수를 전개하는 것은 매우 어려운 일이다. (재귀적으로 전개할 수 있으나, 분모/분자의 숫자가 엄청나게 커져서 제대로 계산하기가 어렵다) 대신에 연분수의 단계를 만들어 나가는 것은 그리 어렵지 않다. 이는 손으로 연분수를 만들어 나가는 과정을 써 보고, 다음 단계의 값을 이전 단계로부터 어떻게 만들 수 있는가를 써보고 이 과정을 점화식으로 재구성한다.

점화식만 재구성할 수 있으면 각 단계의 가수, 분모, 분자를 알 수 있고, 이 값이 반복되는지를 알아보기 위해 맵에 넣은 후 확인해보면 된다.

연분수 추적하기

문제에서 제공한 연분수 풀이 과정을 참고해보자.

  1. a0은 4와 \frac{1}{\sqrt{23} - 4}가 되었다. N = 23이므로, 이 때 4는 N의 제곱근을 넘지 않는 가장 큰 정수가 된다. 이 4를 M이라 하면 남은 분수 부분은 \frac{1}{\sqrt{N} - M}이 되었다.
  2. 위에서 남은 분수 부분의 분모, 분자에 √N + M 을 곱해주면 \frac{\sqrt{N} + M}{N - M^2} 이 된다.  이 값의 가수는 결국 2M \div (N-M^2)의 몫이 된다. 이 몫을 Q라하자. 그러면 남은 식은 Q + \frac{N+M-Q(N-M^2)}{N-M^2} 가 될 것이다.
  3. 따라서 a1은 이 때의 Q인 1이 되고, 남은 식을 다시 같은 식으로 정리하여 가수 영역을 계속 취할 수 있다.
  4. 만약 같은 부분이 계속 나타난다면 다음 차의 분수를 이루는 값들을 비교하면 된다.

초기항 계산

N이 주어질 때, 최초의 정수부분의 값은 해당 N의 제곱근보다 작거나 같은 최대 정수이며, 연분수 항을 추적하기 위한 분수식은 다음과 같은 모양이라고 가정하자.

\sqrt{N}  =  M +  \frac{\sqrt{N} + B}{C}

이때 최초의  A0 와 B, C는 M, -M, 1 로 정의된다. (\sqrt{N} = M + (\sqrt{N} - M) ) 이 최초의 A0는 정수 부분이고, B C 는 다음 단계의 연분수를 만들기 위한 기본자료이다. 제곱근을 연분수로 표현할 때는 분자가 항상 1이 되는 모양이므로, M이후의 값의 역수를 취하고 이를 분모로 하는 분수식으로 쓰게 된다.

\sqrt{N} = M + \frac{1}{\frac{C}{\sqrt{N} + B}}

이 때의 분자 부분을 처음의 식 (A + \frac{\sqrt{N} + B}{C})과 같은 꼴로 만들어내는 과정을 반복하여 연분수를 확장해 나갈 수 있다. 따라서 이 값을 다음값으로 어떻게 바꿔나갈지를 정리해보자.

\frac{c_0}{\sqrt{N} + b} 는 다음과 같이 전개해나가면 다음 단계를 만들 수 있다.

\begin{array}{rl}\frac{C}{\sqrt{N} + B} = & \frac{C(\sqrt{N}-B))}{N-B^2} = A_1 + \frac{(\sqrt{N} - B) - A_1(N - B^2) \div C)}{(N - B^2) \div C} \\ = & A_1 + \frac{\sqrt{n}-B-A((N - B^2) \div C)}{(N - B^2) \div C} \\ = & A_1 + \frac{\sqrt{N} + B_1}{C_1}\end{array} 

따라서 점화식은 다음과 같이 구해진다.

  • C_{n+1} =  (N - B_n^2) \div C)
  • A_{n+1} = | \sqrt{N} - B \div C_{n+1} | = (A_0 - B_n) \div C_{n+1}
  • B_{n+1} =  -B_n - A_n \cdot C_{n+1}

Swift 풀이

점화식이 나왔으니 깔끔하게 풀어보도록하자. 먼저 연분수 표기를 반복하면서 출력해주는 타입을 하나 정의한다.

struct Squared: Sequence, IteratorProtocol {
  let n : Int
  let m : Int
  var a : Int
  var b : Int
  var c : Int
  init(number n: Int) {
    self.n = n
    m = Int(sqrt(Double(n)))
    (a, b, c) = (m, -m, 1)
  }

  mutating func next() -> Int {
    let c = (n - self.b * self.b) / self.c
    // 만약 제곱수라면...
    if c == 0 { return nil }
    let a = (m - self.b) / c
    let b = -self.b - (a * c)
    (self.a, self.b, self.c) = (a, b, c)
    return a
  }
}

제대로 동작하는지 보고싶다면, 아래와 같이 테스트 해볼 수 있다. 이 타입은 정수부분은 생략하고, 그 이후 부분을 계속 출력할 것이다.

for (i, e) in Squared(number: 23).enumerated() where i < 30 {
    print(e)
}

그런데, 숫자값만 출력해서는 순환하는지 여부를 확인할 수 없다. (a, b, c) 순서쌍을 모두 받아서 이전에 나왔던 조합인지를 판별해야 한다. 하지만 튜플은 캐시에 불리하므로, 캐시를 위한 노드 타입을 하나 추가로 만든다.

struct CalcNode: Hashable, CustomStringConvertible {
  let a: Int, b: Int, c: Int
  var hashValue: Int { return description.hashValue }
  var description: String { return "\(a):\(b):\(c)" }
  static func == (left: CalcNode, right: CalcNode) -> Bool {
    return left.a == right.a && left.b == right.b && left.c == right.c
  }
}

그리고 다시 Int? 타입이 아닌 CalcNode? 타입을 리턴하도록 next() 메소드를 수정한다.

struct Squared: ...
  mutating func next() -> CalcNode? {
  ...
  return CalcNode(a:a, b:b, c:c)
 }

이제 문제를 풀어볼 차례이다. 한 번 나온 조합이 반복될 때까지, 나왔던 노드들을 모두 쌓아두었다가 그 개수를 리턴하는 함수를 만들고, 이 결과로 필터링하면 끝이다.

func lengthOfCircularUnits(of n: Int) -> Int {
  let x = Squared(number: n)
  var c = Set<CalcNode>()
  for a in k {
    if !c.contains(a) { c.insert(a) }
    else { return c.count }
  }
  return 0
}

print((2...10000).filter{ lengthOfCircularUnits(of: $0) % 2 == 1 }.count)

파이썬 풀이

파이썬에서는 튜플로 처리할 수 있기 때문에 간단히 요약하는 수준으로 정리했다.

def process(n):
  m = int(n**0.5)
  x = (m, m, 1)
  def next(a, b, c):
    z = (n - (b**2)) // c
    if z is 0: return None
    x = (m - b) // z
    y = -b - (x*z)
    return x, y, z

  s = set()
  while True:
    if x is None: return 0
    if x not in s: s.add(x)
    else: return len(s)
    x = next(*x)

def solve():
  print(len(1 for x in range(2, 10_000) if process(x) % 2 is 1))