【Project Euler】Problem 65: ネイピア数eに収束する分数

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題： ネイピア数eに収束する分数

原文 Problem 65:Convergents of e

問題の要約：　ネイピア数eの連分数展開の100番目の項まで収束した分数の分子の各桁の和を求めよ

連分数に関しては連分数(Wikipedia)に説明があるので参考にしてください。

まずネイピア数eの連分数展開ですが、とても不思議なことに規則性を持っているようです。その辺を確認するためにpython sympyのcontinued_fraction_iteratorを使ってみます。すると問題の中でも提示されているように$[\dots,1,2k,1,\dots]$のパターンが確認できました。

from sympy.ntheory.continued_fraction import continued_fraction_iterator
from sympy import E, Rational, fraction

cf = continued_fraction_iterator(E)
print([next(cf) for i in range(22)])
#[2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, 14, 1]

問題を解くだけであればsympyにcontinued_fraction_convergentsという連分数から収束する分数を求める関数があるので、これを使えば以下のようにN番目の分数を求めることが出来ます。

from sympy.ntheory.continued_fraction import continued_fraction_iterator, continued_fraction_convergents 
from sympy import E, Rational, fraction

it = continued_fraction_convergents(continued_fraction_iterator(E))
conv = [fraction(next(it)) for i in range(10)]
print(conv[-1])
#(1457, 536)

これだけだと味気ないのでConvergentの部分を自分で書いてみます。アルゴリズムは前出のWikipediaの連分数の性質という章を参考にしました。

連分数の性質
いま、$a0$は整数、それ以外の$a_n$は正の整数であるような数列

 $\ \ \ a_0,a_1,a_2,a_3,\dots$

があるとき、数列$p_n,q_n$を以下のように定める

\left\{
\begin{array}{ll}
p_0=1  \\
p_1=a_0 \\
p_n=a_{n-1}p_{n-1}+p_{n-2} & (n\ge2)
\end{array}
\right.\\

\left\{
\begin{array}{ll}
q_0=0  \\
q_1=1 \\
q_n=a_{n-1}q_{n-1}+q_{n-2} & (n\ge2)
\end{array}
\right.

このとき、連分数は

[a_0;a_1,a_2,\dots,a_{n-1}]=\frac{a_{n-1}p_{n-1}+p_{n-2}}
{a_{n-1}q_{n-1}+q_{n-2}}=\frac{p_n}{q_n}

これをそのままPythonのコードにして、連分数のリストから収束する分数を求める関数cf2frac書いてみました。

import itertools
def cf2frac(cf):
  if len(cf) == 1: return (cf[0],1)
  pn_2, pn_1, qn_2, qn_1, an_1 = 1, cf[0], 0, 1, cf[1]
  for n in itertools.count(2):
    pn, qn = an_1 * pn_1 + pn_2, an_1 * qn_1 + qn_2
    if n >= len(cf): break
    an_1, pn_1, pn_2, qn_1, qn_2 = cf[n], pn, pn_1, qn, qn_1
  return (pn, qn)

N = 10
a = [2]+[2 * (i // 3 ) if (i % 3) == 0 else 1 for i in range(2,1+N) ]
print(a)

(num, den) = cf2frac(a)
print(len(a),(num, den))
#[2, 1, 2, 1, 1, 4, 1, 1, 6, 1]
#10 (1457, 536)