오일러 프로젝트 66

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

오일러 프로젝트 66 더보기

오일러 프로젝트 57

오일러 프로젝트 57 번 문제는 2의 제곱근을 연분수 형태로 전개하는 것을 단계별로 진행할 때, 각 단계를 기약분수 형태로 만들어서 분자와 분모의 숫자 자리수를 비교하는 문제이다. 일단 자리 수가 엄청나게 커지는 관계로 쉬운 문제는 아니다.

오일러 프로젝트 57 더보기

오일러 프로젝트 64

사실 64번까지 진도가 나가려면 한참 멀었는데, 이미 써둔지 일년 넘은 글이 계속 썩고 있어서 관련된 오일러 프로젝트 64번 문제 풀이와 함께 발행한다.

2의 제곱근은 다음과 같은 연분수의 꼴로 전개 가능하다.

 1 + \frac {1} {2 + \frac {1} {2 + \frac {1} {2 + \frac {1} {2 + ... } }}}

이때 정수 부분은 1이고, 연분수에서는 다시 분자는 1, 정수부분은 2가 나오는 패턴이 반복된다. 모든 자연수의 제곱근을 연분수 전개꼴로 풀어내면 연분수 내의 정수 부분은 1개 이상의 숫자가 반복되는 패턴이 된다. 이를 [1;(2)]와 같은 식으로 줄여서 쓴다. 이 때 임의의 자연수 N에 대해서 제곱근을 연분수 전개했을 때 이러한 반복패턴을 찾는 방법을 알아보자.

특정한 차수에서 연분수 항(?)은 보통 아래와 같은 식으로 계산된다.

A _{n} + \frac {\sqrt{N} + B _{n}} {C _{n}}

여기서 A _{n} 은 정수부가 되고, B _{n}은 제곱근과 붙어다니는 정수항, C _{n} 은 나머지 무리수 항의 분모가 된다.

여기서 다음 차수를 계산할 때는 정수부를 제외한 값의 역수를 취하여 다시 위의 꼴로 변환하게 된다. 따라서 다음의 모양의 식에서 시작한다.

\frac { C_{n} } { \sqrt{N} + B _{n}}

이 항의 분모를 정수로 만들기 위해서는 곱셈공식을 이용해서 근호가 있는 항을 분자로 옮긴다.

\frac { C _{n} ( \sqrt{N} - B _{n} ) } { N - B _{n} ^ 2 }

이 때 C _{n} N - B _{n} ^ 2과 약분되기 때문에, 위 식은 다시

\frac { ( \sqrt{N} - B _{n} ) } { \left( N - B _{n} ^ 2 \right) \div C _{n}}

따라서 C_{n+1}은 아래와 같이 얻을 수 있다.

C_{n+1} = \frac{N - B_n^2}{C_n}

다시 식은 조금 더 간단히 아래와 같이 정리된다.

\frac { ( \sqrt{N} - B _{n} ) } { C _{n+1}}

여기서 가분수형태의 분자에서 정수부를 취하기 위해서는 \sqrt{N} 을 넘지 않는 최대의 정수 RN을 이용해서 (이는 보통 A _{0} 이다.) 정수부를 떼어낼 수 있고, 이것이 곧 A_{n+1}이 된다.

A_{n+1} = (( RN - B _{n} ) \div C _{n+1})

따라서 이 식은 최종적으로 다음과 같이 정리된다.

A_{n+1} + \frac { ( \sqrt{N} - B _{n} - A_{n+1} \times C _{n+1} ) } { C _{n+1}} 

최종적으로 점화식을 다시 정리하면

A_{n+1} = (( RN - B _{n} ) \div C _{n+1}) \\ B _{n+1} = - B _{n} - A_{n+1} \times C _{n+1} \\ C _{n+1} = \frac {N - B _{n} ^ 2 } { C _{n} } \\

위 점화식을 이용한 오일러프로젝트 64번 문제의 풀이는 아래와 같다.

def continued_fraction(n:int) -> [int]:
    NR = int(n **0.5)
    if NR * NR == n:
        return [NR]
    a = NR
    b = -NR
    c = 1
    cache = [(a, b, c)]
    result = [a]

    def get_next_fraction(b, c):
        cn = (n - b*b) // c
        an = (NR - b) // cn
        bn = -b - an * cn
        return (an, bn, cn)

    while True:
        p = get_next_fraction(b, c)
        if p in cache:
            return result
        cache.append(p)
        a, b, c = p
        result.append(a)

def p64():
    xs = [continued_fraction(i) for i in range(2, 10001)]
    print(sum((1 for x in xs if len(x) % 2 == 0)))

%time p64()
# 1322
# Wall time: 956 ms