2019年3月6日水曜日

ネイピア数の計算

円周率の計算からの続き)

ネイピア数は自然対数の底 $e$ のこと。欧米では Euler's Number とよばれることが多いとあるが,Eulerが関わる数はたくさんあってややこしいので,ここではネイピア数とする。

円周率の計算方法というか表現方法はたくさんあって,その中には与えられた級数の少数の項で良い精度の近似を与えるものがあった。ところで,ネイピア数の表式をざっと調べても,なんだかぱっとしないのである。
\begin{equation}
\begin{aligned}
e^x &= 1 + x + \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots\\
e &= 1 + 1 +  \frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots
\end{aligned}
\end{equation}
を多少変形したものか,連分数で表現したものはあるけれど,ラマヌジャンの公式のようにはっとするものが見当たらない。単に調査不足だけかもしれないし,そもそも指数関数を微分や積分しても元に戻るだけなのが原因かもしれないが,よくわからないのだった。

5項くらいで精度が70桁程度になる公式はないものなのだろうか?Improving the Convergence of Newton's Series Approximation for $e$ とかあったけど,目的のものではない。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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 cf(x,a,b,n)
# continued fraction
  for k = n:-1:1
    x = a[k] + b[k]/x
  end
  return x
end

function napier1(n)
# by series expansion
  o=BigFloat(1)
  sum=o
  for k in 1:n
    sum=sum+o/f[k]
  end
  println(sum-exp(o))
end

function napier2(n)
# by series expansion
  o=BigFloat(1)
  sum=o
  for k in 1:n
    sum=sum+o/((BigFloat(4)^k)*f[k])
  end
  println(sum*sum*sum*sum-exp(o))
end

function napier3(n)
# by continued fraction
  o=BigFloat(1)
  z=BigFloat(2)
  a=ones(BigFloat,1000)
  for k in 1:1000
    a[k]=z*(z*k-z-o)
  end
  a[1]=o
  a[2]=o
  b=ones(BigFloat,1000)
  b[1]=z
  println(cf(z,a,b,n)-exp(o))
end

function napier4(n)
# sqrt(e) by series expansion
  o=BigFloat(1)
  z=BigFloat(2)
  sum=BigFloat(0)
  for k in 0:n
    sum=sum+BigFloat(4*k+3)/(z^(2*k+1)*f[2*k+1])
  end
  println(sum*sum-exp(o))
end

fact(f,160)
@time napier1(57)
@time napier2(41)
@time napier3(26)
@time napier4(24)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
  0.018460 seconds (30.73 k allocations: 1.488 MiB)
-2.763573937630222280123632596096127862757120116619610043207585110453949377011977e-76
  0.029598 seconds (49.70 k allocations: 2.358 MiB)
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
  0.054272 seconds (100.99 k allocations: 5.090 MiB)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-76
  0.032685 seconds (72.55 k allocations: 3.431 MiB)

Juliaの数学的定数と精度へ続く)

0 件のコメント: