오일러 프로젝트 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

오일러 프로젝트 64

Share this:

이것이 좋아요:

관련