Project Euler Problem 57 - 本当は怖いHPC

何年ぶりだろうか、Project Eulerを再開した。今回はProblem 57。

問題は、

${ \displaystyle A_1= 1 + 1/2 = 3/2 }$

${ \displaystyle A_2 = 1 + 1/(2 + 1/2) = 7/5 }$

${ \displaystyle A_3 = 1 + 1/(2 + 1/(2 + 1/2)) = 17/12 }$

${ \displaystyle A_4 = 1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 }$

という数列を考えたとき、1000番目までの数の中に、分子の桁数が分母の桁数より多いような分数はいくつ登場するか、ということ。最初に登場するのは ${A_8= 1393/985}$ である。

まず、最初の1を加算しているところが邪魔なので、1を除いた分数を ${B_n}$ と置く。

${ \displaystyle A_n = 1 + B_n }$

すると、

${ \displaystyle B_n = \frac{1}{2+B_{n-1}} }$

となることが容易にわかる。次に、

${ \displaystyle B_n = \frac{N_n}{D_n} }$

とおいてやると、簡単な式変形により

${ \displaystyle N_{n+1} = D_n }$

${ \displaystyle D_{n+1}= 2D_n + N_n }$

とわかる。あとは、これを計算するだけだ。実際に計算すべき分数は ${1 + \frac{N_n}{D_n}}$ であることと、約分をする必要があることに注意してコードを書けば、

import sys

if len(sys.argv) >= 2:
    Len = int(sys.argv[1])
else:
    Len = 1000

D = [0] * (Len+1)
N = [0] * (Len+1)

def gcd(m, n):
    if m < n: return gcd(n, m)

    r = m % n
    while r > 0:
        m, n = n, r
        r = m % n
    return n

def reduce(a, b):
    d = gcd(a, b)
    return a/d, b/d

count = 0

for n in range(1, Len+1):
    if n == 1:
        N[1] = 1
        D[1] = 2
    else:
        N[n] = D[n-1]
        D[n] = N[n-1] + 2 * D[n-1]
    nr, dr = reduce(N[n] + D[n], D[n])

    if len(str(nr)) > len(str(dr)) :
        count += 1

print count

これを実行すれば答えが表示される。 ${N_n}$ も消せるけど、まぁいいや。