2012年11月19日月曜日

プロジェクトオイラー 問65

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


問65「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)


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

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


2.関数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]で3つごとに余り1になればTrue、他はFalseです。
 これを前の{}内の辞書要素のうち当てはまる値を設定します。
 //演算子は割り算の商です。True、つまり、3で割り余り1になる場合、
 「(i//3+1)*2」で上記式の2kにあたる部分の値を設定し、その他は1を設定します。


3.関数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を逆順に1つずつ取り出しループ変数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から算出したその近似分数の分子と分母を返します。

4.関数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を文字型にして、1桁ずつ数値にしたリストをMとします。


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


5.関数の外
・n = 100
 print f(n)
 問題に合うように100を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
272,
6963524437876961749120273824619538346438023188214475670667,
2561737478789858711161539537921323010415623148113041714756

(追記)
・20130127 ネタバレ注意の追記など

2012年11月11日日曜日

プロジェクトオイラー 問64

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

問64「N ≤ 10000 について奇数の周期をもつ連分数はいくつあるか?」
平方根は連分数の形で表したときに周期的であり, 以下の形で書ける:
           1
N = a0 + --------------------
                1
         a1 + ---------------
                     1
              a2 + ----------
                          1
                   a3 + -----
                         ...

例えば, √23を考えよう.
                             1              1
√23 = 4 + 23 - 4 = 4 + ------------ = 4 + ---------------
                             1                 √23 - 3
                          --------         1 + --------
                          √23 - 4                  7
となる.
この操作を続けていくと,
              1
√23 =  4 + --------------------------
                   1
             1 + ---------------------
                        1
                  3 + ----------------
                             1
                       1 + -----------
                                   1
                             8 + -----
                                  ...
を得る.
操作を纏めると以下になる:
a0 = 4, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a1 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a2 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a3 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)
a4 = 8, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a5 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a6 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a7 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)

よって, この操作は繰り返しになることが分かる.
表記を簡潔にするために,
√23 = [4;(1,3,1,8)] と表す.
(1,3,1,8)のブロックは無限に繰り返される項を表している.

最初の10個の無理数である平方根を連分数で表すと以下になる.
√2=[1;(2)], period=1
√3=[1;(1,2)], period=2
√5=[2;(4)], period=1
√6=[2;(2,4)], period=2
√7=[2;(1,1,1,4)], period=4
√8=[2;(1,4)], period=2
√10=[3;(6)], period=1
√11=[3;(3,6)], period=2
√12= [3;(2,6)], period=2
√13=[3;(1,1,1,1,6)], period=5
N ≦ 13で奇数の周期をもつ平方根は丁度4つある.

N ≦ 10000 について奇数の周期をもつ平方根が何個あるか答えよ.
-----



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






私の解答例は以下です。畳んでいます。
def g(n, b, c, L, M):
	if (b, c) in M: return L
	if n==b*b: return L
	M.append((b, c))
	C, t = (n-b*b)/c, int(n**(1.0/2))
	u = (b+t)%C
	L.append((b+t)//C)
	return g(n, t-u, C, L, M)

def f(n):
	t = 0
	for i in xrange(n+1):
		s = g(i, int(i**(1.0/2)), 1, [], [])
		if len(s)%2: t += 1
	return t

n = 10000
print f(n)


以下のように考えました。
  c                     √n-B
------ ・・・(A),   A + ----- ・・・(B)
√n-b                     C
として、(A)を連分数(B)にすることを考えます。
なお、(B)で計算したAの値を循環部分として保存し、
B, Cの値を次世代のb, cとして再帰的に計算します。
まず(A)を有理化します。
  c          c*(√n+b)     c*(√n+b)
------ = --------------- = ----------- ・・・(A)'
√n-b     (√n-b)*(√n+b)    (n-b*b)
連分数展開後の(B)の形に合わせるようにて分母Cを決めます。
    (n-b*b)
C = ------- ・・・(C)
       c
(A)'に(C)を当てはめ、
(√n+b)
------- ・・・(D)
   C
(D)から、連分数(B)のAに相当する部分をくくりだします。
そのために、√nの整数部分tとして、
t = int(√n) ・・・(E)
(b+t)   √n-t
----- + ------ ・・・(F)
  C       C
ですが、連分数(B)のAは整数のため、
Aは(F)の第1項の整数部分で、余りuは第2項に含めます。
A = (b+t)//C ・・・(G)
u = (b+t)%C  ・・・(H)
(G)(H)を(F)に反映して
 (b+t)        √n-t +u
{----- - u} + ---------
   C               C
             √n - (t-u)
= (b+t)//c + -----------
                 C
以上より、A,B,Cをまとめておくと以下の通りです。
(E)(H)のt,uを使って、
A = (b+t)//C ・・・(G)
B = t-u      ・・・(I)
C = (n-b*b)/c・・・(C)

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

1.関数g(n, b, c, L, M)
  c                        √n-B
------ を連分数展開して、A + ------ (分数部分の絶対値は1未満)
√n-b                         C
に式変形したときの、循環部分のAのリストLを取得します。
リストMにはbとcの組、タプル(b,c)を貯めます。
例)√23 = [4;(1,3,1,8)]の場合、[1,3,1,8]を返します。
・関数内に自分自身を呼び出す記述のある再帰関数です。

・if (b, c) in M: return L
 まず再帰の終了条件です。
 引数b,cの組が今までにでてきていれば、循環部分がひと回りしたことになるので処理を終了し循環部分のAのリストLを呼び出し元に戻します。

・if n==b*b: return L
 この条件の場合、入力値の分母が0でゼロ割り算エラーのため連分数にならないのでこの時点でのリストLを返却します。

・M.append((b, c))
 ここまでくれば、引数b,cは連分数の循環しない部分なので、タプル(b,c)をリストMに貯めます。

・C, t = (n-b*b)/c, int(n**(1.0/2))
 上記説明の(C)(E)を実装しました。
 **演算子は累乗で、平方根を求める部分を1/2乗として記述しました。
 整数で1/2とするとint型のため値が0になるので、1.0とfloat型にすることで
 1/2がfloat型の0.5になります。分母の2を2.0にしても同じことです。
 共通モジュールのmathのsqrt関数を使っても同じことです。
 import math
 t = math.sqrt(n)

・u = (b+t)%C
 上記説明の(H)を実装しました。

・L.append((b+t)//C)
 循環部分リストLには、上記説明の(G)を実装しました。

・return g(n, t-u, C, L, M)
 連分数展開を続けます。
 現世代のB,Cが次世代のb,cなので、上記説明の(C)(I)を実装します。


2.関数f(n)
・問題に合う、n以下で奇数の周期を持つ平方根の個数を返します。

・t = 0
 個数カウンタtを初期クリアしておきます

・for i in xrange(n+1):
 ループ変数iとして、0以上n+1未満(n以下)の整数を発生させます。

・s = g(i, int(i**(1.0/2)), 1, [], [])
 関数gで連分数展開したときの循環部分リストを取得しsとします。
 第1引数は連分数展開する平方根の元の値なので、ここではループ変数iです。
 第2引数は√iからくくりだせる最大整数なので、ここでは√iの整数部分です。
 第3引数には1、第4第5引数には空のリストを設定します。

・if len(s)%2: t += 1
 %演算子は割り算の余りです。
 len(s)は関数gで取得した連分数循環部分リストsの要素数です。
 この要素数が奇数ということは2で割った余りが1で、1は論理的にTrueなので、上記if文が成り立ちます。このとき、個数カウンタtを1つカウントアップします。

・return t
 変数iのループが終了したら、個数カウンタを呼び出し元に返します。

3.関数の外
・n = 10000
 print f(n)
 問題に合うように10000を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1322


(追記)
・20121208 式(C)の分子分母が逆だったので修正。説明文の中に2箇所あり。
 誤)C = c/(n-b*b)・・・(C)
 ↓
 正)C = (n-b*b)/c・・・(C)
・20130127 ネタバレ注意の追記など

2012年11月2日金曜日

プロジェクトオイラー 問63

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

問63 「n 乗して得られる n 桁の正整数は何個あるか?」
5桁の数 16807 = 75は自然数を5乗した数である.
同様に9桁の数 134217728 = 89も自然数を9乗した数である.

自然数をn乗して得られるn桁の正整数は何個あるか?

-----

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






私の解答例は以下です。畳んでいます。
def f():
	import math
	return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]

L = f()
print len(L),L


自然数aをb乗してb桁の正整数ということは、桁数の区切りは10の累乗値なので以下の式が成り立ちます。
0 ≦ 10(b-1) ≦ ab < 10b ・・・(A)
a>1、b>1 ・・・(B)

まず、aの変域を求めます。(A)の右部分の2式より、
1 ≦ a < 10 ・・・(C)

次にbの変域を求めます。
(A)の左部分をbについて解きます。
まず、左部分の2式の対数をとって、
(b-1)*log10 ≦ b*(log a)

次にbを左辺に寄せて、
(b-1)/b ≦ (log a)/log10

対数の底を10にして、
1-(1/b) ≦ log10a

1-log10a ≦ 1/b

b ≦ 1/(1-log10a)

(B)よりb>0なので、bの変域は
1 ≦ b ≦ 1/(1-log10a) ・・・(D)

上記(C)(D)の範囲で二重ループさせます。

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

1.関数f()
・import math
 数学用関数モジュールmathをインポートして使えるようにしておきます。

・return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]
 内包表記内に二重ループを入れています。
 まず前の「for a in xrange(1,10)」で、aが1以上10未満の整数でループします。
 そしてループ変数aを後ろのfor文へ渡します。
 後ろの「for b in xrange(1, int(1/(1-math.log(a, 10)))+1)」で、
 bが1以上1/(1-math.log10a)未満の整数でループします。
 math.log()関数は対数関数です。第2引数は対数の底です。
 int関数で切り捨てて整数にします。
 このa,bの値を先頭の式、(a,b,a**b)に渡します。
 これは問題に合うa, b, aのb乗の値のタプルです。
 このタプルを内包表記でリストにため、そのまま呼び出し元に返します。

2.関数の外
・L = f()
 print len(L),L
 関数f()で問題に合うa, b, aのb乗のタプルを格納したリストを取得しLとします。
 求める個数、とそのリストを表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
49

(追記)
・20130127 ネタバレ注意の追記など