BigMath::PIより早く円周率を計算する

ベンチマークを取ってみる

まずはライバルとなるBigMath::PIの円周率100桁はどれくらい時間がかかるか計測してみます。
（本当はそれぞれ複数回試行して平均値を出したり桁数をもっと大きくしないといけないんですが…）

bigmath_pi.rb

require 'bigdecimal/math'
require 'benchmark'

result = Benchmark.realtime do
  puts BigMath::PI(100).to_s
end

puts "Time: #{result}s"

$ ruby bigmath_pi.rb
0.3141592653589793238462643383279502884197169399375105...
Time: 0.0004811960000097315s

ふむふむ、早いですね。
参考までに普通のMath::PIを計測してみると…

math_pi.rb

require 'benchmark'

result = Benchmark.realtime do
  puts Math::PI
end

puts "Time: #{result}s"

$ ruby math_pi.rb
3.141592653589793
Time: 3.3455988159403205e-05s

0.000033sってことでさすがに早いですね。

次に、自分でBigDecimalを使ってガウス＝ルジャンドルのアルゴリズムを書いてみました。

pi.rb

require "bigdecimal"
require "benchmark"

result = Benchmark.realtime do
  prec = 100
  conv = 6
  a = BigDecimal("1")
  b = BigDecimal("1") / BigDecimal("2").sqrt(prec)
  t = BigDecimal("1") / 4
  p = BigDecimal("1")

  for n in 1..conv do
    an = (a + b) / 2
    b = (a * b).sqrt(prec)
    t -= p * (an - a) * (an - a)
    p *= 2
    a = an
  end

  puts (a + b) * (a + b) / (4 * t)
end

puts "Time: #{result}s"

$ ruby pi.rb
0.3141592653589793238462643383279502884197169399375105...
Time: 0.03566918699652888s

うーむ、桁数の割にめちゃくちゃ遅い。

計算方法を見てみる

中身が知りたくなったので、それぞれの計算式を見てみました。

Math::PI

普通のMath::PIは、math.hに直接書かれているM_PIの数字を利用しているようです。
なので計算は特にしていません。

math.h

# define M_PI  3.14159265358979323846  /* pi */

BigMath::PI

BigMath::PIはBigMathモジュールにあるPIメソッドで円周率を計算しています。
どうやらマチンの公式を使っているようです。

がんばって早くする

案1）計算精度を調整する

収束させようとすればするほど平方根の桁数が多くなっているみたいだったので、指定桁数までで丸めることにしてみました。

pi.rb

require "bigdecimal"
require "benchmark"

result = Benchmark.realtime do
  prec = 100
  conv = 6
  b = BigDecimal::limit(prec)
  a = BigDecimal("1")
  b = BigDecimal("1") / BigDecimal("2").sqrt(prec)
  t = BigDecimal("1") / 4
  p = BigDecimal("1")

  for n in 1..conv do
    an = (a + b) / 2
    b = (a * b).sqrt(prec)
    t -= p * (an - a) * (an - a)
    p *= 2
    a = an
  end

  puts (a + b) * (a + b) / (4 * t)
end

puts "Time: #{result}s"

$ ruby pi.rb
0.3141592653589793238462643383279502884197169399375105...
Time: 0.0011667709986795671s

だいぶ早くなったものの、そもそも平方根の計算が遅いみたいで、BigMath::PIには勝てないという結果に。

案2）GMPを使ってみる

平方根の計算が遅いなら、GMP(GNU Multi-Precision Library) を使って桁数の多い小数点計算を最適化すればいいじゃない。ということで試してみました。

gmp_pi.rb

require 'gmp'
require "benchmark"

result = Benchmark.realtime do
  prec = 400
  conv = 6
  a = GMP::F(1, prec)
  b = GMP::F(1, prec) / GMP::F(2, prec).sqrt
  t = GMP::F(1, prec) / 4
  p = GMP::F(1, prec)

  for n in 1..conv do
    an = (a + b) / 2
    b = (a * b).sqrt
    t -= p * (an - a) * (an - a)
    p *= 2
    a = an
  end

  puts (a + b) * (a + b) / (4 * t)
end

puts "Time: #{result}s"

$ ruby gmp_pi.rb
0.3141592653589793238462643383279502884197169399375105...
Time: 0.00019221799999513678s

勝った！
でもgem使っててなんか卑怯な気が…

まとめ

Chudnovskyの公式を使ってみたりもしたのですが、桁数が少なすぎて案2とあまり差が出ませんでした。
このくらいの桁数ならGMPの計算パワーがすべてということですね。

参考

下記のサイトを参考にさせてもらいました。ありがとうございます。

円周率πの計算：ガウス＝ルジャンドルのアルゴリズム - Hatada's Home Page
Ruby - 円周率計算（Chudnovsky の公式使用）！ - mk-mode BLOG