※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.