黄金比1000桁(2) -- N桁の(正しい)黄金比を得るには何番目の項で比をとればよいか
フィボナッチ数列で黄金比を計算する時、N桁の黄金比を正確に出すには何番目の隣り合う2項で比をとれば計算できるのだろうか。とか思ったのでメモ
d = 1000桁の(整数化した)黄金比
factor = 10^999
fib(n) = n番目のフィボナッチ数列
f(x) = 1000 - log(d - factor*fib(x+1)/fib(x), 10) # n番目の隣り合う項で比を出したときの正確な桁数
とする。あらかじめ、1000桁の整数化した黄金比を、前回のプログラムを使って"out"という名前で保存しておく。
# coding: utf-8 # g.py def fib(): a = b = 1 while True: yield a a,b = b,a+b def f(proc): it = fib() x = it.next() y = it.next() for i in xrange(2400): proc(i+1, x, y) x = y y = it.next() d = int(open("out").read()[:1000]) factor = 10**999 import math def printprecision(i, x, y): # 何桁正確かを表示する a = abs(d - factor*y/x) if a != 0: b = math.log(a, 10) else: # math domain errorを防止 b = 0 print "%d %d" % (i, 1000 - int(b)) #print "%d %f" % (i, 1000-b) def main(): f(printprecision) main()
このプログラムの出力を"ggg.dat"という名前で保存してgnuplotでプロットした。
正確な桁数は比を出す項の番号に比例するらしい。
なので、この直線の傾きを出して、項の番号を引数にとって正確な桁数を返す関数precition(x)を求められれば、precition(a)=NとしてN桁の黄金比を求めるのに必要な項の番号aがわかるはず。
まず、傾きを出す。
500番目で黄金比をとるとき、正確な桁数はなのでこれで直線の傾きがでた。この値を係数にもつ直線の関数は。
この式で合っているかN=1000桁で試す。
y = 1000のとき、
正解は2389。いい線いってるか?ずれすぎな感じもする。
もっと正確さを高めるために2000番目でも傾きを計算して平均をとってみた。
改良式
y = 1000のとき
さらに精度を高めるために最小二乗法で係数を計算してみた。まず上のプログラムの出力を整数化せずに浮動小数点のまま出力して"gggg.dat"という名前で保存した。そして次のプログラムを実行した。
#! /usr/bin/python # coding: utf-8 lines = open("gggg.dat").readlines() x = [] y = [] for line in lines: c = [float(a) for a in line.split()] x.append(c[0]) y.append(c[1]) def sumwith(seq, f): acc = 0 for x in seq: acc += f(x) return acc def square(n): return n**2 def mulp(p): return p[0]*p[1] def minsqr(xs, ys): """最小自乗法""" n = len(xs) A = sumwith(ys, square) B = sumwith(xs, square) C = sumwith(ys, lambda x:x) D = sumwith(zip(xs, ys), mulp) E = sumwith(xs, lambda x:x) den = n*B - E**2 a = (n*D-C*E) / den b = (B*C-D*E) / den return a,b def main(): print minsqr(x, y) main()
実行結果
(0.41795541216418891, 0.66647261522383361)
改良式2
y = 1000のとき
あれ?精度下がった??まあこれでいいや。
桁数 | x(小数点以下打ち切り) | 正解 | 差 |
---|---|---|---|
1000 | 2392 | 2389 | 3 |
2000 | 4785 | 4785 | 0 |
3000 | 7177 | 7177 | 0 |
4000 | 9570 | 9571 | -1 |
5000 | 11962 | 11962 | 0 |
6000 | 14355 | 14355 | 0 |
7000 | 16748 | 16747 | 1 |
8000 | 19140 | 19143 | -3 |
9000 | 21533 | 21532 | 1 |
10000 | 23925 | 23924 | 1 |
ってことで、関数precision(x)はy = 0.41795541216418891xで決定!間違った非効率なやり方の気がしてきたけど気にしない!
つづく。