円周率の計算

（電子計算機と人間からの続き）

現代科学シリーズの「電子計算機と人間」では，円周率の計算のためのFORTRANプログラムが書かれていたが，高校１年生の自分にはなぜそれで円周率が求まるのかがさっぱり分からなかった。50年ぶりに，図書館でその本を手にとると，たしかに分かりやすくはなかったが，アークタンジェントのマクローリン展開に１を代入した単純な公式なのであった。
\begin{equation}
\begin{aligned}
\arctan(x) & = x - \frac{x^3}{3!} +  \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots \\
\frac{\pi}{4} & =  1 - \frac{1}{3!} +  \frac{1}{5!} - \frac{1}{7!} + \cdots
\end{aligned}
\end{equation}
この公式は収束がとても悪い見本のような式なので，1000項でようやく3.14が0.1%精度で求まる始末である。

そこで，いくつかの代表的な公式をJuliaプログラムに落としてみた。BigFloatの精度がeps(BigFloat) = $10^{-77}$程度なので，有効数字は77桁弱である。そのBigFloatの使い方が，いまひとつわかっておらず適当に入れているが，たぶん冗長だと思われる。

円周率の計算で，有効数字77桁程度の精度を出すには，昔よくみかけたマチンの公式で 54項，なんだかすごいのだけれど今回初めて使った ラマヌジャンの公式で ９項，スマートなガウス＝ルジャンドルの方法で５項が必要であった。もっと楽しそうなJuliaでのπの話はこちら。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)

function fact(f,n)
# calculation of factorial 1!..n!
  for i in 2:n
    f[i]=BigFloat(i)*f[i-1]
  end
# println(f[n])
end

function pie1(n)
# Gregory-Leibniz Formula
# pi=4*arctan(1)
  sum=BigFloat(0)
  for i in 1:n
    sum=sum+(-1)^(i-1)/BigFloat(2*i-1)
  end
  println(4*sum-pi)
end

function pie2(n)
# Machin's Formula
# pi = 16*arctan(1/5)-4*arctan(1/239)
  sum=BigFloat(0)
  sun=BigFloat(0)
  for i in 1:n
    sum=sum+(-1)^(i-1)/((2*i-1)*BigFloat(5)^(2*i-1))
    sun=sun+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
  end
  println(16*sum-4*sun-pi)
end

function pie3(n)
# Takano-Kikuo's Formula
# pi = 48*arctan(1/49)+128*arctan(1/57)-30*arctan(1/239)+48*arctan(1/110443)
  sum1=BigFloat(0)
  sum2=BigFloat(0)
  sum3=BigFloat(0)
  sum4=BigFloat(0)
  for i in 1:n
    sum1=sum1+(-1)^(i-1)/((2*i-1)*BigFloat(49)^(2*i-1))
    sum2=sum2+(-1)^(i-1)/((2*i-1)*BigFloat(57)^(2*i-1))
    sum3=sum3+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
    sum4=sum4+(-1)^(i-1)/((2*i-1)*BigFloat(110443)^(2*i-1))        
  end
  println(48*sum1+128*sum2-20*sum3+48*sum4-pi)
end

function pie4(n)
# Ramanujan Formula
  sum=BigFloat(1103)
  for k in 1:n
    sum=sum+f[4*k]/(f[k]^4*BigFloat(396)^BigFloat(4*k))*BigFloat(1103+26390*k)
  end
  println(BigFloat(9801)/(sqrt(BigFloat(8))*sum)-pi)
end

function pie5(n)
# Chudnovsky formula 
  sum=BigFloat(13591409)
  for k in 1:n
    sum=sum+f[6*k]/(f[3*k]*f[k]^3*BigFloat(-640320)^BigFloat(3*k))*BigFloat(13591409+545140134*k) 
  end
  println(BigFloat(640320)^1.5/(BigFloat(12)*sum)-pi)
end

function pie6(n)
# Gauss-Legendre algorithm
  z=BigFloat(2)
  a=BigFloat(1)
  b=1/sqrt(z)
  t=1/(z*z)
  p=BigFloat(1)
  for k in 1:n
    o=a
    a=(o+b)/z
    b=sqrt(o*b)
    t=t-p*(a-o)*(a-o)
    p=z*p
  end
  println((a+b)*(a+b)/(4*t)-pi)
end

fact(f,160)
@time pie1(10000000)
@time pie2(54)
@time pie3(22)
@time pie4(9)
@time pie5(5)
@time pie6(5)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-9.999999999999975000000000000312499999999990468750000000541015625007904590936847e-08
  5.580711 seconds (60.04 M allocations: 3.131 GiB, 4.69% gc time)
3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
  0.055656 seconds (76.18 k allocations: 3.601 MiB, 7.32% gc time)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
  0.071662 seconds (123.15 k allocations: 5.740 MiB)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
  0.057192 seconds (93.11 k allocations: 4.357 MiB)
-3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
  0.068310 seconds (106.84 k allocations: 5.084 MiB)
0.0
  0.036871 seconds (53.77 k allocations: 2.544 MiB)

（ネイピア数の計算に続く）