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


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

```ruby:bigmath_pi.rb
require 'bigdecimal/math'
require 'benchmark'

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

puts "Time: #{result}s"
```

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

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

```ruby:math_pi.rb
require 'benchmark'

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

puts "Time: #{result}s"
```

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

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

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

```ruby: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"
```

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

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

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

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

```c:math.h
# define M_PI  3.14159265358979323846  /* pi */
```

#### BigMath::PI
`BigMath::PI`はBigMathモジュールにあるPIメソッドで[円周率を計算](https://github.com/ruby/ruby/blob/ruby_2_2/ext/bigdecimal/lib/bigdecimal/math.rb#L173-L216)しています。
どうやらマチンの公式を使っているようです。

## がんばって早くする
#### 案1）計算精度を調整する
収束させようとすればするほど平方根の桁数が多くなっているみたいだったので、指定桁数までで丸めることにしてみました。

```ruby: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"
```

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

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

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

```ruby: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"
```

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

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

## まとめ

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

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

[円周率πの計算：ガウス＝ルジャンドルのアルゴリズム  - Hatada's Home Page](http://home.a00.itscom.net/hatada/math/pi/pi01mj.html)
[Ruby - 円周率計算（Chudnovsky の公式使用）！ - mk-mode BLOG](http://www.mk-mode.com/octopress/2015/05/08/ruby-pi-computation-by-chudnovsky-bsa-with-gmp/)