Julia: マチンの公式による円周率の計算

マチンの公式による円周率の計算

https://ja.wikipedia.org/wiki/マチンの公式

マチンの公式

計算

マチンの公式を

の形にし，arctan(x) をグレゴリー級数に直して，それぞれ最初の方の項だけを計算して差を取る。

wikipedia では 各項の和を取る範囲を違ったものにしているが，以下のプログラムでは同じにした。

また指定した桁が確保できるまで収束計算を続ける。収束までの反復回数は「0.72 * 指定桁」である（経験則で導いた）。

  function machin(precision=15; trace=false)
      if precision > 15
          f1, f2, f4, f5, f16, f239, N = big(1), big(2), big(4), big(5), big(16), big(239), 0.72precision
          setprecision(BigFloat, ceil(Int, precision/log10(2)))
      else
          f1, f2, f4, f5, f16, f239, N = 1, 2, 4, 5, 16, 239, 11
      end
      f1byf5 = f1 / f5
      f1byf239 = f1 / f239
      x1 = x2 = x = 0
      for n in 0:N
          alpha = f2 * n + f1
          sig = iseven(n) ? f1 : -f1
          x1 += f16 * sig / alpha * f1byf5 ^ alpha
          x2 +=  f4 * sig / alpha * f1byf239 ^ alpha
          new_x = x1 - x2
          trace && println("$(n+1): $new_x")
          x == new_x && return new_x
          x = new_x
      end
      println("not converged!")
  end

  machin (generic function with 2 methods)

浮動小数点数 Float64 では 11 回の反復回数でほぼ正しい結果が得られる。

  println(machin(trace=true))
  println(1pi)

  1: 3.18326359832636
  2: 3.1405970293260603
  3: 3.1416210293250346
  4: 3.1415917721821773
  5: 3.1415926824043994
  6: 3.1415926526153086
  7: 3.141592653623555
  8: 3.1415926535886025
  9: 3.141592653589836
  10: 3.1415926535897922
  11: 3.141592653589794
  12: 3.141592653589794
  3.141592653589794
  3.141592653589793

指定した桁数より多めに計算されるが，結果はかなり正確である。

  println(machin(50))
  println(big(pi))

  3.141592653589793238462643383279502884197169399375101
  3.141592653589793238462643383279502884197169399375101

  length("141592653589793238462643383279502884197169399375101")

  51

  20000 桁の計算では，末尾 4 桁が異なる。

  @time a = machin(20000);

   39.491599 seconds (1.91 M allocations: 11.740 GiB, 0.17% gc time)

  println(string(a)[end-50:end])

  296810620377657883716690910941807448781404907552024

  println(string(big(pi))[end-50:end])

  296810620377657883716690910941807448781404907551788

  @time b = machin(20010);

   39.650767 seconds (2.36 M allocations: 15.209 GiB, 0.22% gc time)

  println(string(b)[end-50:end])

  776578837166909109418074487814049075517820385654186

Julia: ガウス=ルジャンドルのアルゴリズムによる円周率

https://ja.wikipedia.org/wiki/ガウス=ルジャンドルのアルゴリズム により円周率を求める

初期値

反復式

a, b が希望する精度（桁数）になるまで，以下の計算を繰り返す。

収束は二次収束（反復のたびに正しい桁数が直前のものの 2 倍になるので，小数第 n 位まで求めるときには，log2(n) 回程度の反復で十分である。

π の算出

上で求めた a, b, t により，以下の式で近似できる。

以上に基づき，Julia で，任意の精度で計算できるようにプログラムを書いた。収束過程もトレースできる。

Julia では（Python も同じであるが），反復式に示された 4 個の数式を，タプルの演算式として 1 行で書くことができる。 ただし，t の更新に a の更新値が使われているので，式を重複して指定している。

  function Gauß_Legendre(precision=15; trace=false)
      if precision > 15
          one, two, four, n = big(1), big(2), big(4), ceil(Int, log2(precision))
          setprecision(BigFloat, ceil(Int, precision/log10(2)))
      else
          one, two, four, n = 1, 2, 4, 3
      end
  
      a, b, t, p = one, one/sqrt(two), one/four, one
      for i in 1:n
          a, b, t, p = (a+b)/two, sqrt(a*b), t - p*(a - (a+b)/two)^2, two*p
          trace && println("i = $i $((a+b)^2/(four*t))")
      end
      (a+b)^2/(four*t)
  end

  Gauß_Legendre (generic function with 2 methods)

倍精度浮動小数点数 Float64 は有効数字は 16 桁弱なので，log2(16) = 4 回の反復が必要そうであるが，実際には 3 回で十分である。

  println(Gauß_Legendre(trace=true))
  println(big(π))

  i = 1 3.1405792505221686
  i = 2 3.141592646213543
  i = 3 3.141592653589794
  3.141592653589794
  3.141592653589793238462643383279502884197169399375105820974944592307816406286198

  100 桁でも実際は 5 回の反復で正しく求められている。

  println(Gauß_Legendre(100, trace=true))
  println(big(π))

  i = 1 3.14057925052216824831133126897582331177344023751294833564348669334558275803490290782728762155276690067
  i = 2 3.14159264621354228214934443198269577431443722334560279455953948482143476722079526469464344891799130592
  i = 3 3.14159265358979323827951277480186397438122550483544693578733070202638213783892739903141694204346905842
  i = 4 3.14159265358979323846264338327950288419711467828364892155661710697602676450064306171100657772659806854
  i = 5 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862562870321167203607
  i = 6 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807
  i = 7 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807
  3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807
  3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807

  println(Gauß_Legendre(1000))
  println(big(π))

3.1415926535897932384626433832795028841971693993751058...9375195778185778053217122680661300192787661119590921642019898
3.1415926535897932384626433832795028841971693993751058...9375195778185778053217122680661300192787661119590921642019898

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