プロジェクトオイラー 問６５

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。


問６５「e についての連分数である近似分数の100項目の分子の桁の合計を求めよ」
2の平方根は無限連分数として書くことができる。
            1
√2 = 1 + -------------------
                1
          2 + ---------------
                    1
              2 + -----------
                        1
                  2 + -------
                      2 + ...
無限連分数である√2 = [1;(2)]と書くことができるが、
(2)は2が無限に繰り返されることを示す。
同様に、√23 = [4;(1,3,1,8)]。
平方根の部分的な連分数の数列から良い有理近似が得られることが分かる。

√2の近似分数について考えよう。
    1   2
1 + - = -
    2   3
      1     7
1 + ----- = -
        1   5
    2 + -
        2
      1         17
1 + --------- = --
          1     12
    2 + -----
            1
        2 + -
            2
      1             41
1 + ------------- = --
          1         29
    2 + ---------
             1
        2 + -----
                1
            2 + -
                2
従って、√2の近似分数からなる数列の最初の10項は：
1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408, 1393/985, 3363/2378, ...

もっとも驚くべきことに、数学的に重要な定数、
e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]。
eの近似分数からなる数列の最初の10項は：
2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71, 1264/465, 1457/536, ...
10項目の近似分数の分子の桁を合計すると1+4+5+7=17である。
eについての連分数である近似分数の100項目の分子の桁の合計を求めよ。
－－－－－


注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。





私の解答例は以下です。畳んでいます。

def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a

def h(n):
	if n<0: return 0
	L = [2]
	for i in xrange(n-1):
		L.append({True:(i//3+1)*2, False:1}[i%3==1])
	return L

def g(L):
	b, c = L.pop(), 1
	for a in L[::-1]:
		x, y = a*b+c, b
		s = gcd(x, y)
		b, c = x//s, y//s
	return b, c

def f(n):
	L = h(n)
	x, y = g(L)
	s = str(x)
	M = [int(s[i]) for i in xrange(len(s))]
	return sum(M), x, y

n = 100
print f(n)

４つの関数を使っています。

１．関数gcd(a, b)
・「ユークリッドの互除法」を使用して、引数a,bの最大公約数を返します。
　説明は問５を参照してください。

２．関数h(n)
・n項までのネイピア数の連分数リストを返します。
　但し問題文ではセミコロンで示している区切り文字も全てカンマ区切りにします。
　e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]
　　→[2, 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]

・if n<0: return 0
　nが0未満は未定義なので0を返します。

・L = [2]
　初期値は連分数の固定部分なので2です。

・for i in xrange(n-1):
　ループ変数iには0以上n-1未満のn個の整数が順に設定されます。

・L.append({True:(i//3+1)*2, False:1}[i%3==1])
　「1,2k,1, ...」を発生させ、リストLに設定します。
　%演算子は割り算の余りです。
　[i%3==1]で３つごとに余り１になればTrue、他はFalseです。
　これを前の｛｝内の辞書要素のうち当てはまる値を設定します。
　//演算子は割り算の商です。True、つまり、３で割り余り１になる場合、
　「(i//3+1)*2」で上記式の2kにあたる部分の値を設定し、その他は1を設定します。

３．関数g(L):
・リストLとして連分数の固定部分aと循環部分b1,b2,...を受け取り、
　その近似分数x/yの分子x、分母yを返します。式で書くと以下の通りです。
      1             x
a + ------------- = -
           1        y
    b1 + --------
               1
         b2 + ---
              ...

・b, c = L.pop(), 1
　リストのpop関数はリストの要素を取り出して、それをリストから削除します。
　引数未指定の場合は最後の要素を取り出して、それをリストから削除します。
　初期値として、bに連分数の循環部分の最後、cに1を設定します。

・for a in L[::-1]:
　リストLを逆順に１つずつ取り出しループ変数aに設定します。

・x, y = a*b+c, b
    c   a*b+c   x
a + - = ----- = -
    b     b     y
上記式を実装しました。

・s = gcd(x, y)
　計算したx，yが約分できることがあるので、関数gcdで最大公約数を求めておきます。

・b, c = x//s, y//s
　先ほど求めたx,yをその最大公約数で割り、できる限り簡単な分数に直します。
　これを次のb,cとしてループ変数aのループを最後まで回します。

・return b, c
　連分数リストLから算出したその近似分数の分子と分母を返します。

４．関数f(n):
・ネイピア数のn項目までの連分数の近似分数の分子の各桁の合計、及びその分子、分母を返します。

・L = h(n)
　関数h(n)でn項までのネイピア数の連分数リストを取得しLとします。

・x, y = g(L)
　関数g(L)で上記リストLの近似分数の分子と分母を取得し、x,yとします。

・s = str(x)
　L = [int(s[i]) for i in xrange(len(s))]
　上記で求めた分子xを文字型にして、１桁ずつ数値にしたリストをMとします。

・return sum(M), x, y
　分子の各桁の合計値、分子、分母を返します。

５．関数の外
・n = 100
　print f(n)
　問題に合うように100を関数fに渡し戻り値を表示します。
　
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
272,
6963524437876961749120273824619538346438023188214475670667,
2561737478789858711161539537921323010415623148113041714756
（追記）
・20130127　ネタバレ注意の追記など