黄金比1000桁(3) -フィボナッチ数列で効率良く求めてみる-

programming python

前回のエントリでは、黄金比をフィボナッチ数列の隣り合う2数の比を使って求めるとき、N桁の精度を求めるには何番目の隣り合う2数をとればいいかを考えた。そして、隣り合う項のx番目で比をとったときの、比の黄金比に対する正確な桁数を求める関数を、precision(x) = 0.41795541216418891xとした。

今回は、それを使って前前回に作ったプログラムよりも効率良くN桁の黄金比を求めるプログラムを書いた。

# coding: utf-8
# egoldenratio.py
def fib():
    a = b = 1
    while True:
        yield a
        a,b = b,a+b

def inverseprecision(n):
    return n / 0.41795541216418891

def gratio(n):
    it = fib()
    x = int(inverseprecision(n))
    for _ in xrange(x-1): it.next()
    x = it.next()
    y = it.next()
    factor = 10**(n-1)
    return factor*y/x

def main():
    import sys
    if sys.argv.__len__() != 2:
        print "usage: %s prec" % __file__
        sys.exit()
    n = int(sys.argv[1])
    print gratio(n)

if __name__ == '__main__':
    main()

実行

$ python egoldenratio.py 1000
1618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088075386891752126633862223536931793180060766726354433389086595939582905638322661319928290267880675208766892501711696207032221043216269548626296313614438149758701220340805887954454749246185695364864449241044320771344947049565846788509874339442212544877066478091588460749988712400765217057517978834166256249407589069704000281210427621771117778053153171410117046665991466979873176135600670874807101317952368942752194843530567830022878569978297783478458782289110976250030269615617002504643382437764861028383126833037242926752631165339247316711121158818638513316203840052221657912866752946549068113171599343235973494985090409476213222981017261070596116456299098162905552085247903524060201727997471753427775927786256194320827505131218156285512224809394712341451702237358057727861600868838295230459264787801788992199027077690389532196819861514378031499741106926088674296226757560523172777520353613936

整形

$ python egoldenratio.py 1000 | python -c "import sys,re; lst = re.findall(r'(\d{10,10}|\d+)', sys.stdin.read()); print '\n'.join(['%4d | ' % (i*50+1) + ' '.join(lst[i*5:i*5+5]) for i in range(len(lst)/5)]); print '       ' + ' '.join(lst[len(lst)/5*5:len(lst)/5*5+5])" 
   1 | 1618033988 7498948482 0458683436 5638117720 3091798057
  51 | 6286213544 8622705260 4628189024 4970720720 4189391137
 101 | 4847540880 7538689175 2126633862 2235369317 9318006076
 151 | 6726354433 3890865959 3958290563 8322661319 9282902678
 201 | 8067520876 6892501711 6962070322 2104321626 9548626296
 251 | 3136144381 4975870122 0340805887 9544547492 4618569536
 301 | 4864449241 0443207713 4494704956 5846788509 8743394422
 351 | 1254487706 6478091588 4607499887 1240076521 7057517978
 401 | 8341662562 4940758906 9704000281 2104276217 7111777805
 451 | 3153171410 1170466659 9146697987 3176135600 6708748071
 501 | 0131795236 8942752194 8435305678 3002287856 9978297783
 551 | 4784587822 8911097625 0030269615 6170025046 4338243776
 601 | 4861028383 1268330372 4292675263 1165339247 3167111211
 651 | 5881863851 3316203840 0522216579 1286675294 6549068113
 701 | 1715993432 3597349498 5090409476 2132229810 1726107059
 751 | 6116456299 0981629055 5208524790 3524060201 7279974717
 801 | 5342777592 7786256194 3208275051 3121815628 5512224809
 851 | 3947123414 5170223735 8057727861 6008688382 9523045926
 901 | 4787801788 9921990270 7769038953 2196819861 5143780314
 951 | 9974110692 6088674296 2267575605 2317277752 0353613936

前前回のプログラムと実行速度を比較してみた。比較はコマンドラインでtimeコマンドのrealの値を使った。
環境はcore2 1.3GHz, mem2GB, Ubuntu。-は6分以上かかったので面倒くさくなったとこ。

桁数 今回のプログラム 前前回のプログラム
1000 0.02 0.3
2000 0.03 2.2
3000 0.03 7.
4000 0.05 16.
5000 0.06 -
6000 0.08 -
7000 0.1 -
8000 0.12 -
9000 0.15 -
10000 0.18 -

最後に10万桁を出してみた。

$ time python egoldenratio.py 100000 > /dev/null

real	0m14.265s
user	0m14.261s
sys	0m0.000s

それなりの速度で計算できた。