ニュートン法で黄金比を求める[Python]

　まだ慣れないPythonで書いたのでミスや書き方がPythonっぽくないかもしれないが許してつかあさい。

1.ニュートン法とは

平たく言えば漸化式にもとづいて計算していくと答えが出てくるというもので，細かい説明をはしょると次の式を用います。

X0にf(x) = 0の解に近い値を入れて計算をしていくと，f(x) = 0 の解に収束するという素晴らしい漸化式です(関数によって必ずしも収束するとは言い切れない)。

　例えばf(x) = x^2 - 2 を代入した次の式は有名。

大学入試の相加平均・相乗平均とかの問題に出てきそうな形してますね。X0 = 1.5 を代入すると，X1=1.41666……，X2=1.41421568……と，どんどんx^2 - 2 = 0の解である√2に近づいていきます。

2.黄金比とは

いろいろな書き方がありますが，黄金比とは，関数f(x)

の正の解とします。つまり，黄金比φは

と表せる事ができるわけです。黄金比については面白い性質が沢山あったりするわけですが，プログラムを書いて，さらにその話を書くにはちょっと余白が狭すぎます。

3.求める方法

1，2で説明したことをくっつければいいわけですから，1のf(x)に，2で示したf(x)を代入します。すると

という漸化式が得られます。φ=1.6180339……ですから，X0 = 1.5ぐらいで始めれば答えは出そうです。

　まぁちなみに手元の電卓を使ってX0=1.5を代入すると，X1 = 1.625，X2 = 1.618055……となっていますので正しそうです。

　ここで電卓でやればいいじゃん，という声が出てきそうですね。ですが今回は任意の桁数を計算したいのです。1000桁とか表示してくれる電卓があればいいのですがそうも行きません。

　なので一工夫を加えます。まずXnを次のように置きます。

ここでan,bnは整数です。なぜこんなことをしたのかというと，コンピュータは小数を扱うのが苦手です。なので小数を整数の組み合わせ，つまり分数で扱おうという作戦です。
　これを代入すると次のようになります。

4.ソースコード

前置きが長かったです。つまり言いたいことはこういうこと。本当にプログラミング言語は便利な言葉やなぁ。

# coding: utf-8
# 2015/10/15 黄金比を求めるプログラム

# python2.xは下の一行がいる
from __future__ import print_function
import math

if __name__ == "__main__":
   # 分母分子を表す変数を用意
   b0 = 3
   a0 = 2
   bn = 0
   an = 0

   # 桁数をN=1000とする
   N = 1000

   # 誤差は2次収束する
   # つまり正しい桁数が1回の計算で2倍になっていく
   # よってN桁正しい値を得るためにはlog[2]Nとなる
   # ここでは多めにlog[1.5]Nにした
   for i in range(0,int(math.log(N,1.5))):
      an = a0*(2*b0-a0)
      bn = a0*a0+b0*b0
      
      a0 = an
      b0 = bn
   
   # 整数部は消す
   bn = (bn - an)*10
   print("1.")
   
   # 一桁ずつ表示する
   for j in range(N):
      print(int(bn/an),end="")
      if( j % 10  == 9):
         print(" ",end="")
      if( j % 50 == 49):
         print()
      bn = (bn % an)*10
   

　実行結果