N88-BASICで円周率 (3回目)

2022/12/27(火)
N88-BASICで円周率 (3回目)
 
レオンハルト・オイラー(Leonhard Euler)の展開式
Tan^-1(x) = x/(1+x²)Σ[n=0～∞]Π[k=1～n]
{(2kx²)/(2k+1)(1+x²)} (n = 0の項は1とする)
より
 
π = 4Tan^-1(1)
= 2Σ[n=0～∞]Π[k=1～n]{k/(2k+1)}
(n = 0の項は1)
= 2(1 + 1/3 + 1･2/(3･5) + 1･2･3/(3･5･7) + …)
= 2(1 + 1/3(1 + 2/5(1 + k/(2k+1)(1 + … ))))
= 2 + 1/3(2 + 2/5(2 + k/(2k+1)(2 + … ))))
 
π = 2 + 1/3(2 + 2/5(2 + k/(2k+1)(2 + … ))))
これを基数倍して整数部分を取出していく
 
以後、基数を10として解説します
(2桁ずつ取出す時は基数100になる)
 
[/]:商、%:余、[]:ガウス記号として
k/b { a }
= k[a / b] + k/b { (a % b) }
を利用して帯分数に変換
(1/3{20} → 1･6+1/3{2})
以下、これの繰返し
 
d + k/(2k+1)(a + … )
のk/(2k+1) × aを1未満に変換して
dに加えていく
 
つまり
10πを
d_k-1 + k/(2k+1)(d_k + … )
= 10d_k-1 + k[10a/(2k+1)]
+ k/(2k+1)(10d_k % 10 + … )
{k = n～1}
に変換していく
d₀ = 31となるので
3を表示してd₀ = 1
としてさらに10(10π)して
を繰返す
10倍する毎に{k = n～1}の
nを減らすことができる
 
 
桁数hと項数nの関係
 
nと桁数hとの関係
3.1のときh = 2桁
10^-x = 10^-2 = 0.01の桁まで求めておくとすると
x = hとなる
 
k/(2k+1) ＜ 1/2 (k≧1)を示す
2k / (2k+1) ＜ 1
2k ＜ 2k+1 (k≧1)
 
2Π[k=1～n]{k/(2k+1)} ＜ 2･(1/2)ⁿ ＜ (1/10)^x 
(1/2)･2ⁿ ＞ 10^x 、2^n-1 ＞ 10^x 
log₂2^n-1 ＞ log₂10^x 、n-1 ＞ xlog₂10
n ＞ xlog₂10 + 1
 
k/(2k+1) ≧ 1/3 (k≧1)を示す
3k / (2k+1) ≧ 1
3k ≧ 2k+1 (k≧1)
 
2Π[k=1～n]{k/(2k+1)} ≧ 2･(1/3)ⁿ ≧ (1/10)^x 
(1/2)･3ⁿ ≦ 10^x 、3ⁿ ≦ 2･10^x 
log₃3ⁿ ≦ log₃2･10^x 、n ≦ log₃2 + xlog₃10
 
項数nと10^-xの関係(桁数h = x)
n ＞ xlog₂10 + 1
n ≦ xlog₃10 + log₃2
この関係式は
必要最低限の条件ではなく
無駄を含む十分な条件なので
効率が悪い事は、ご了承ください
 
log₂10 = 3.321928094… ＜ 34/10
log₃10 = 2.095903274… ＞ 2
log₃2  = 0.630929753… ＞ 0
(log_ab = ln b/ln a)
より
 
n ＞ 34x/10 + 1 … 余分に見積もる
(必要な項数の計算に使用、x:桁数)
n ≦ 2x　　　　 … 少なく見積もる
(破棄する項数の計算に使用、x:基数10^x)
 
NL-BASIC(N88-BASIC互換?)で走らせ
1040桁まで正しく表示されている事を確認しました
それ以降の桁の確認はしていません
 
上記式やプログラムのミスなどにより
正しく計算出来ていない可能性はありますので
使用には十分注意して下さい
 
N88-BASIC互換?VL,NL,XL-BASICと
blg～.zip(pi001.bas)は
以下のリンクからダウンロードできます

https://ulprojectmail.blogspot.com
Readme.txtを読んで遊んで下さい