超越数π
超越数(transcendental number)は有理係数の代数方程式の解ではない複素数である。
複素数は大部分が超越数である。自然対数の底e(Napier/Euler number)と円周率πは超越数である。
今回はπについて情報科学の観点からprogrammingを利用して調べる。
1.πを含む幾何公式
円の面積、A=πr^2, 円周はC=2πr, r=半径, (1)
ここで記号”^”はべき乗を示す。πを左辺に移すと、
π=A*r^-2, π=C*r^-1. (1a)
記号”*”は乗算を示す。
2.π生成級数(series)
複数の級数が知られている。15世紀より知られているMDL(Mādhava,Leibniz)級数がprogramming的に扱い易い。自然数をn(n=1,2,3...)と書く。
Σ{(-1)^n/(2n+1)}={1 -1/3 +1/5 -1/7 +…}=π/4, (2)
π=4Σ{(-1)^n/(2n+1)}, (2a)
記号Σはn=1から∞まで総和を計算すること。
3.質問
式(1a)と(2a)の違いは何ですか？
半径rが式(2a)に無い。
級数に半径なんて無いから、当然と思いますか？
でも円の面積、周の計算式も半径は自分たちで設定したものでしょう？
例えば、半径を直径に変えても面積、周は計算できます。
それなら、級数Σ{(-1)^n/(2n+1)}の中に半径に相当するprogramming上の操作が無いのでしょうか？
無いという証明は困難ですが、逆の「有る」ということは、何か1つでも例を見つければ良いです。
4.Computer計算の特徴
Computerでは無限を取り扱えません。これは何時も考えておかなければならないことです。Σ記号は∞まで加算することなので、級数値を計算する場合、有限の値に収束するように見えても無限大に発散することがあります。例Σ(1/n)→∞。
Computer hardwareから見ると、整数、実数の内部表現(2進数)は違います。
Hardware資源の制約から、整数は4B(byte)、実数は4/8Bです。
取り扱える数の範囲は、
整数：[-2147483648, 2147483647];
実数8B：[4.94065645841246544*10^-324, 1.79769313486231570*10^308] and
[-1.79769313486231570*10^308, -4.94065645841246544*10^-324]です。
実数4Bは表現範囲が狭くなるので、今回の目的には使用しません。
実数8Bの有効数字は上記では22桁もあるように思われるかも知れませんが、実際は15桁、実務上の計算では12桁ぐらいに低下します。
特に絶対値の大きい数に小さい数を加減算すると低下します。丸め誤差といって、加減算器で必ず生じる問題です。
問題を極小化するために、加減算では絶対値の小さい順番に演算する必要があります。
今回、頻繁に現れる総和計算では注意しましょう。
5.MDL級数の変形、計算処理
MDL級数を変形しpointer"i”の更新を±1とします。
Σ{sgn/(n-0.5)^M}={1/0.5^M -1/1.5^M +1/2.5^M -1/3.5^M +…}=π/2(in case of M=1),
sgn={1, -1, 1, -1, 1,…}, M={1,2,3,…}, (2c)
Mはparameterです。M=1のときはπ/2ですが、M=3では(π^3)/4となります。
さらにζ関数を意識して、
Σ{1/n^M}={1/1^M +1/2^M +1/3^M +1/4^M +…}, (2d)
という級数(＠級数と言います)も用意します。式(2d)は式(2c)でsgn={1,1,1,…} & 0.5→0とした場合に等価です。
ζ(1)は∞なので、この場合はM≧2です。M=2,ζ(2)=(π^2)/6.
Σを「十分大きいnまで」総和計算する、と変えます。
実習用computerでは5*10^7ぐらいが適切です。今後10^7をO(7)とも書きます。
O(-7)は小数点以下7桁の0が続くという意味です。
有効数字7桁の意味で(-7)の精度とも言います。
Computerの実数計算は近似計算なので、1.0000001を「O(-6)の精度の」整数1と取り扱うことがあります。
実数から整数への型変換処理では注意が必要です。
整数型変数n=1.0000001はn=1です。しかし、整数型変数n=0.9999999ではn=0です。
6.ゼータ関数
Programming codeの精度検証のために、Riemann zeta functionを式(2d)から5*O(7)のΣ範囲で計算します。
ζ(2)=0.1666 666=(π^2)/6.
正しく計算されています。
7. MDL級数
変形ζ(), MDL級数、式(2d,2c)のMを大きくしていくとπの様々な姿が見えてきます。
M=4までの結果を示します。
表１.ζ関数とMDL級数のparameter Mとの関係
1段目：ζ関数(標準計算式)
2段目：MDL級数 共にΣは5*O(7)まで
---べき乗(M)---------値------------- π^2比
1.0000 18.304749238294 0.539183
1.0000 0.693147170560 14.238830
1.5000 2.612092505974 3.778428
1.5000 0.765147024624 12.898965
2.0000 1.644934046848 6.000000<---ζ()
2.0000 0.822467033424 12.000000<---MDL
2.5000 1.341487257249 7.357211
2.5000 0.867199889012 11.381003
3.0000 1.202056903160 8.210597
3.0000 0.901542677370 10.947462
3.5000 1.126733867317 8.759481
3.5000 0.927553577774 10.640468
4.0000 1.082323233711 9.118907
4.0000 0.947032829497 10.421607
-----------------------------------------------
25.0000 1.000000029804 9.869604
25.0000 0.999999970199 9.869605
===============================================
表1の考察１：
M=2で、級数の"n”の偶奇の値をB,Aとすると、
ζ(),MDLの値が2:1なので、A+B=2(A-B),
ゆえに3B=A=1.23370054 である。
検証：実際に級数計算のpointerをmod2で2分して別々に計算すると、
A=0.411233506712, even sum
B=1.233700540136, odd sum.
A=(π^2)/24.0000006,
B=(π^2)/8.00000006.
となる。O(-6)の精度である。
表1の考察２：
M=25で
A+B=A-B=1, B=0, A=1 となるのは、
奇数sumでは1/1^M +1/3^M +… →1.
偶数では1/2^M +1/4^M+ …→0.
であるから、M→∞ではB=0となるため。
表１から考えられること：
ζ(2)関数の級数計算のpointerをmod2からmod3,mod4,...とするとどうなるのでしょうか？
次のような結果を得ます。
表２：ζ(2)関数の剰余pointer分離
------------------------------------------------
mod3(0)=0.182770445205,<---(π^2)/54.000002
mod3(1)=1.121733007270, 54=2*3*3*3=6*(3*3)
mod3(2)=0.340430594373.
mod4(0)=0.102808374178,<---(π^2)/96.000005
mod4(1)=1.074833067157, 96=2*3*4*4=6*(4*4)
mod4(2)=0.308425132534,<---(π^2)/32.000005
mod4(3)=0.158867472979. 32=2*16=2*(4*4)
mod5(0)=0.065797358674,<---(π^2)/150.000009
mod5(1)=1.050695084217, 150=2*3*5*5=6*(5*5)
mod5(2)=0.291014259621,
mod5(3)=0.145448382836,
mod5(4)=0.091978961500.
---------------------------------------------------
以上をまとめると、
ζ(2)=π^2/6=(π^2/6)/(1*1),
ζ(2)mod2(0)=(π^2/6)/(2*2),
ζ(2)mod3(0)=(π^2/6)/(3*3),
ζ(2)mod4(0)=(π^2/6)/(4*4),
ζ(2)mod5(0)=(π^2/6)/(5*5),
modQ(0)の”Q”は、πのべき乗を除外すれば、式(1a)的には「半径のような」性質である。