Julia: グレゴリー・ライプニッツ級数で円周率を計算

グレゴリー・ライプニッツ級数を用いると円周率を計算することができる。

x = 1 のとき

よって，右辺を求めて 4 倍すれば π の近似値が得られる。

ここで，あわてて（式2）に基づいてプログラムするのはよくない。

  function f1(n)
      x = 0
      for i in 1:n
          x += (-1) ^ (i + 1) / (2i - 1)
      end
      4x
  end

  f1 (generic function with 1 method)

一番の問題点は時間がかかりすぎること。

  n = 10^8
  @time f1(n)

    5.474844 seconds (31.41 k allocations: 1.800 MiB, 0.21% compilation time)

  3.141592643589326

時間がかかる主たる原因は，項の符号が + / - 交代にでてくるのを，毎回 (-1) ^ (i + 1) で計算していること。

単純に奇数項は +，偶数項は - とすればよい。やり方は色々あるが，三項演算子を使ってみよう。

  function f2(n)
      x = 0
      for i in 1:n
          x += (isodd(i) ? 1 : -1) / (2i - 1)
      end
      4x
  end

  f2 (generic function with 1 method)

これで，実行速度はほぼ 10 倍になる。

  n = 10^8
  @time f2(n)

    0.598436 seconds (31.04 k allocations: 1.780 MiB, 1.80% compilation time)

  3.141592643589326

別の解決策としては，プラスの項（奇数項）とマイナスの項（偶数項）の和を別々に取ることである。

  function f3(n)
      n = (n ÷ 2) * 2
      plus = minus = 0
      for i in 1:4:2n
          plus += 1/i
          minus += 1/(i + 2)
      end
      4(plus - minus)
  end

  f3 (generic function with 1 method)

  n = 10^8
  @time f3(n)

    0.313891 seconds (41.63 k allocations: 2.349 MiB, 4.08% compilation time)

  3.141592643584488

f(3) は f(2) のほぼ倍速である。

結果は
f2(10^8) = 3.141592643589326
f3(10^8) = 3.141592643584488
であるが，この違いは，演算順序の違いによる丸め誤差の集積の差である。
π の真値は 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
であるから，いずれも 3.1415926 までしか正しくない。
収束は大変遅く，10^k 項までで，有効桁 k の精度である。

  for i in 3:9
      n = 10 ^ i
      println("n = 10^$i $(f2(n))")
  end
  println("pi =     $(big(π))")

  n = 10^3 3.140592653839794
  n = 10^4 3.1414926535900345
  n = 10^5 3.1415826535897198
  n = 10^6 3.1415916535897743
  n = 10^7 3.1415925535897915
  n = 10^8 3.141592643589326
  n = 10^9 3.1415926525880504
  pi =     3.141592653589793238462643383279502884197169399375105820974944592307816406286198