Skip to content

Instantly share code, notes, and snippets.

@metasta

metasta/README.md

Last active October 8, 2018 06:20
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save metasta/fb00e40900fc13e72ed5acdcade8f597 to your computer and use it in GitHub Desktop.
Save metasta/fb00e40900fc13e72ed5acdcade8f597 to your computer and use it in GitHub Desktop.
Download ZIP
πの連分数展開と Python で遊んだ
Raw
README.md

πの連分数展開

※2017年2月頃に調べたことや書いたコードの断片があったので今回まとめた.

経緯

ある日こんな記事を見つけた.

Rubyでn桁の円周率を求める - ほんまの走り書き技術メモ

len = ARGV[0].to_i
B = 10 ** len
B2 = B << 1
pi = (len * 8 + 1).step(3, -2).inject(B) {|a, i| (i >> 1) * (a + B2) / i} - B
puts "3.#{pi}"

（pi.rb）

# time ruby pi.rb 100
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679

real    0m0.045s
user    0m0.027s
sys     0m0.007s

なぜこれでπを計算できるのか全く分からない. 実におもしろい.

解読

上記コード中の pi が inject 内でどう反復処理されるのか, 具体的に len = 4 の場合を想定してみる.

n回目の処理による pi の値を pi_n と書くと、

pi_n 値
pi_0 B ( = 10000)
pi_1 (33>>1) * (pi_0 +B*2) / 33 => 16/33 (2 + 1)B
pi_2 (31>>1) * (pi_1 +B*2) / 31 => 15/31 (2 + 16/33 (2 + 1))B
pi_3 (29>>1) * (pi_2 +B*2) / 29 => 14/29 (2 + 15/31 (2 + 16/33 (2 + 1)))B
... ...
pi_16 ( 3>>1) * (pi_15+B*2) / 3 => 1/3 (2 + 2/5 (2 + 3/7 (2 + 4/9 (... + 16/33 (2 + 1)...)))))B

あ, 連分数か！

π = 2 + 1/3 (2 + 2/5 ( 2 + 3/7 ( ... ( 2 + k/(2k+1) ( ...

この連分数展開と同じ形になってる.

冒頭のコードでは最後に pi から B を引くことで結果的に「(πの近似値 - 3) * 10^len」なる len 桁の整数が得られる. これを "3." のあとに続けて出力する仕組み.

他の連分数展開

仕組みがわかったので遊ぶ. 他の連分数展開でやってみることにした.

π = 3 + 12 / (6 + 32 / (6 + 52 / (6 + 72 / ... / (6 + (2k-1)2 / (...

ここでは仮に「Euler の連分数展開」と呼んでおく.（本当に Euler が発見した式なのかよくわからなかったけど）

この連分数展開（pi_cfrac_euler.py）を実行すると

$ python pi_cfrac_euler.py 20
3.14159312402678768936
$

あれ, 5桁しか合ってない...

精度

この計算法の精度は反復の回数によって決まる. 冒頭のコードではその反復回数を決定するのは (len * 8 + 1).step(3, -2) の部分で, これは [8len + 1, 8len - 1, 8len - 3, ..., 7, 5, 3] なる配列（項数 4len）である. len 桁の精度を出すのに必要な回数が, 冒頭のコードでは 4len 回程度らしい.

「Euler の連分数展開」の収束がどれほど遅いのか、ループ数で評価する.（pi_cfrac_euler_convergence.py）

結果

$ python pi_cfrac_euler1_convergence.py
 loop	pi_value
   10	3.141342106491591265862325425192
   20	3.141561384695304299598041554733
   30	3.141583391796114369868440248953
   40	3.141588746734475866731712289897
   50	3.141590653390850510506784254053
   60	3.141591496102305613909284494475
   70	3.141591924689733015218628331336
   80	3.141592165289509200425817852290
   90	3.141592310643697988118051605269
  100	3.141592403583551514369911113433
  200	3.141592622339597990649466773864
  300	3.141592644330508262825083198933
(...)
 9400	3.141592653589492245425566659517
 9500	3.141592653589501650766836390220
 9600	3.141592653589510668294052788678
 9700	3.141592653589519317791533075938
 9800	3.141592653589527617844692135901
 9900	3.141592653589535585923954082408
10000	3.141592653589543238462018383362
$

（真の値：3.1415926535897932384626433832795...）

100 項で 6 桁まで近似できるが, その後は 10000 項計算しても 12桁止まり. 元のアルゴリズムで 10000 項計算すれば 10000 / 4 = 2500 桁以上の精度が得られるので差は歴然.

さらに他の連分数展開

こんな形の式もある:

4/π = 1 + 12 / (3 + 22 / (5 + 32 / (7 + 42 / ...

これは Brouncker が発見したらしい.

1000桁の精度を達成するために何回のループが必要か, 元のアルゴリズムと Brouncker の連分数展開とを比べた.（pi_cfrac_comparison.py）

結果

$ python pi_cfrac_comparison.py
to reach 1000 digit precision:
pi   = 2+1/3(2+2/5(2+3/7(2+4/9(2+...: 3318 loops
4/pi =1+1^2/(3+2^2/(5+3^2/(7+4^2+...: 1309 loops
$

Brouncker 速い！

Brouncker の収束の速さをループ数で評価する.（pi_cfrac_brouncker_convergence.py）

（上述の Euler の連分数展開の収束とやってることは同じだが, 精度が10桁上がるごとにそれに要したループ数を表示する）

結果

$ python pi_cfrac_brouncker_convergence.py
   15	3.1415926535
   29	3.14159265358979323846
   42	3.141592653589793238462643383279
   55	3.1415926535897932384626433832795028841971
   68	3.14159265358979323846264338327950288419716939937510
   81	3.141592653589793238462643383279502884197169399375105820974944
   94	3.1415926535897932384626433832795028841971693993751058209749445923078164
  107	3.14159265358979323846264338327950288419716939937510582097494459230781640628620899
  120	3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825
  133	3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
$

133回のループで100桁の精度に到達する.

  • 註1：収束は速いが計算は遅い.
  • 註2: Brouncker の式の発見者は Brouncker だが, 証明したのは Euler.
Raw
pi.rb
len = ARGV[0].to_i
B = 10 ** len
B2 = B << 1
pi = (len * 8 + 1).step(3, -2).inject(B) {|a, i| (i >> 1) * (a + B2) / i} - B
puts "3.#{pi}"
Raw
pi_cfrac_brouncker_convergence.py
length = 100
B = 10 ** (length + 1)
B3 = B * 3
BB = B * B
BB4 = BB << 2
# Brouncker
# 4/pi = 1 + 1**2 / (3 + 2**2 / (5 + 3**2 / (7 + 4**2 + ...
def pi_cfrac_brouncker(n):
pi_inv = B
for i in range(2 * n - 1, -1, -2):
pi_inv = B * i + (((i + 1) >> 1) ** 2) * BB / pi_inv
pi = BB4 / pi_inv
pi -= B3
return "3.%d" % pi
precision = 10
pi1000 = '3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989'
match = pi1000[:(precision+2)]
for i in range(50000):
if (precision > length):
break;
pi = pi_cfrac_brouncker(i)
while match in pi:
print("%5d\t%s" % (i, match))
precision += 10
match = pi1000[:(precision+2)]
if (precision > length):
break;
Raw
pi_cfrac_comparison.py
B = 10 ** 1001
B2 = B * 2
B3 = B * 3
BB = B * B
#
# pi = 2 + 1/3 (2 + 2/5 (2 + 3/7 (2 + 4/9 (2 + ...
def pi_cfrac_orig(n):
pi = B
for i in range(2 * n + 1, 1, -2):
pi = (B2 + pi) * (i >> 1) / i
pi -= B
return "3.%d" % pi
# Brouncker
# 4/pi = 1 + 1**2 / (3 + 2**2 / (5 + 3**2 / (7 + 4**2 + ...
def pi_cfrac_brouncker(n):
pi_inv = B
for i in range(2 * n - 1, -1, -2):
pi_inv = B * i + (((i + 1) >> 1) ** 2) * BB / pi_inv
pi = (BB << 2) / pi_inv
pi -= B3
return "3.%d" % pi
pi1000 = '3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989'
print("to achieve 1000-digit precision:")
i = 1
while(1):
pi_i = pi_cfrac_orig(i)
if(pi1000 in pi_i):
print("pi = 2+1/3(2+2/5(2+3/7(2+4/9(2+...: %d loops" % i)
break
i += 1
i = 1
while(1):
pi_i = pi_cfrac_brouncker(i)
if(pi1000 in pi_i):
print("4/pi =1+1^2/(3+2^2/(5+3^2/(7+4^2+...: %d loops" % i)
break
i += 1
Raw
pi_cfrac_euler.py
# Euler
# pi = 3 + 1**2 / (6 + 3**2 / (6 + 5**2 / (6 + 7**2 / ... (6 + (2k-1)**2 / ...
import sys
n = int(sys.argv[1])
B = 10 ** n
pi = B * 3
for i in range(8 * n + 1, -1, -2):
pi = i ** 2 * B * B / pi + B * 6
pi -= B * 6
print("3.%d" % pi)
Raw
pi_cfrac_euler_convergence.py
def pi_cfrac_euler(n):
B = 10 ** 30
pi = B * 3
for i in range(2 * n - 1, -1, -2):
pi = i ** 2 * B * B / pi + B * 6
pi -= B * 6
return "3.%d" % pi
print(" loop\tpi_value")
for i in range(10, 101, 10):
print("%5d\t%s" % (i, pi_cfrac_euler(i)))
for i in range(200, 10001, 100):
print("%5d\t%s" % (i, pi_cfrac_euler(i)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment