黄金比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番目で黄金比をとるとき、正確な桁数はf(500) = 209.638155248なので209.638155248/500 = 0.419276310496これで直線の傾きがでた。この値を係数にもつ直線の関数はy = 0.419276310496x

この式で合っているかN=1000桁で試す。
y = 1000のとき、
x = 1000 / 0.419276310496 = 2385.0620103411261
正解は2389。いい線いってるか?ずれすぎな感じもする。


もっと正確さを高めるために2000番目でも傾きを計算して平均をとってみた。
f(2000) = 836.601075998
836.601075998 / 2000 = 0.41830053799900002
sqrt(0.419276310496 * 0.41830053799900002) = 0.4187881400549836

改良式 y = 0.4187881400549836x

y = 1000のとき
x = 1000 / 0.4187881400549836 = 2387.8422150844763


さらに精度を高めるために最小二乗法で係数を計算してみた。まず上のプログラムの出力を整数化せずに浮動小数点のまま出力して"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 = 0.41795541216418891x
y = 1000のときx = 2392.5997149360078
あれ?精度下がった??まあこれでいいや。

桁数 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で決定!間違った非効率なやり方の気がしてきたけど気にしない!

つづく。