日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。
問６４「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)の第１項の整数部分で、余りuは第２項に含めます。
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)
２つの関数を使っています。
１．関数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)を実装します。
２．関数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とします。
第１引数は連分数展開する平方根の元の値なので、ここではループ変数iです。
第２引数は√iからくくりだせる最大整数なので、ここでは√iの整数部分です。
第３引数には1、第４第５引数には空のリストを設定します。
・if len(s)%2: t += 1
%演算子は割り算の余りです。
len(s)は関数gで取得した連分数循環部分リストsの要素数です。
この要素数が奇数ということは2で割った余りが1で、1は論理的にTrueなので、上記if文が成り立ちます。このとき、個数カウンタtを1つカウントアップします。
・return t
変数iのループが終了したら、個数カウンタを呼び出し元に返します。
３．関数の外
・n = 10000
print f(n)
問題に合うように10000を関数fに渡し戻り値を表示します。
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1322
（追記）
・20121208 式(C)の分子分母が逆だったので修正。説明文の中に２箇所あり。
誤）C = c/(n-b*b)・・・(C)
↓
正）C = (n-b*b)/c・・・(C)
・20130127 ネタバレ注意の追記など
0 件のコメント:
コメントを投稿