二進数コードで弾き語り

情報系大学生の徒然技術ブログ。

Machinの公式を使ってπを近似する

どうも、drumathです。前回の記事からだいぶ日が経ってしまいましたね。
というのも、大学生活が始まり、なかなか執筆に時間がさけなかったのが原因です。
今日は気分で少し書いてみようと思ったのですが、テーマがまたπの近似という結果になってしまいました…(許して)

Machinの公式でπを求める!

冒頭で、大学生活が始まったと述べましたが、約1か月前の解析学の授業をきっかけにこのテーマを書こうと思いまして、
実はその授業の先生がπの近似値の世界記録を持ってるらしく(詳しい話はよく知りませんが...)、そんな経緯からなのでしょうか、逆三角関数の話をしているときにMachinの公式を教えてくださいました。
そもそも、逆三角関数は(辺の比)→(角度)という写像なので、πが絡む公式はいくつかあるらしく、下記のようなものが授業では紹介されました。


 \displaystyle
\begin{align}

2 \tan^{-1} \frac{1}{2} - \tan^{-1} \frac{1}{7} = \frac{\pi}{4} \\

2 \tan^{-1} \frac{1}{3} + \tan^{-1} \frac{1}{7} = \frac{\pi}{4} \\

4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} = \frac{\pi}{4} \\

\end{align}

実は、一番最後の「239!?」ってなるやつがMachinの公式です。

Machinの公式を示してみる

今からこの公式を使うわけですが、本当にこうなるのかとりあえず示します。
これには、


 \displaystyle
\begin{align}
\tan(\tan^{-1}x)=x
\end{align}
という性質を利用します。

まず、Machinの公式の両辺の正接をとります。すると


 \displaystyle
\begin{align}

\tan (4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} ) &= \tan \frac{\pi}{4} \\ &=1

\end{align}
となりますね。よって、証明のゴールは左辺を変形していって1にすれば良いわけですね。やってみましょう!


 \displaystyle
\begin{align}

\tan (4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} ) について、\\
\tan^{-1} \frac{1}{5} = a,
\tan^{-1} \frac{1}{239} = b
とする。\\

\end{align}


 \displaystyle
\begin{align}

また、\\
\tan 2a = \frac{2\tan a}{1-\tan^{2} a}  &= \frac {2 \cdot \frac{1}{5} } { 1-(\frac {1} {5} )^{2} } \\ &= \frac{5}{12}\\ \\
\tan 4a = \frac{2\tan 2a}{1-\tan^{2} 2a} &= \frac{120}{119}  \\ \\

よって\\
\tan(4a-b) = \frac{\tan 4a - \tan b}{1+ \tan 4a \tan b} &= \frac{ \frac{120}{119} - \frac{1}{239} } { 1 + \frac{120}{119} \cdot \frac{1}{239} } \\ &= 1

\end{align}
とまぁ、最後にごりっと計算すると1になるわけです。すごいな。

πを計算するには

さて、証明ができたので、さっそく使っていきましょう
この公式を使ってどうやってπを求めるかというと、、、
積分だ!
待ってたぜぃって感じですが、どちらかというと積分は最後にちょこっとするだけで(というか積分ほぼやらないか…)、メインはマクローリン展開です。

まず、下の式を与えておきます。


 \displaystyle
\begin{align}

(\tan^{-1} x)^{'} = \frac{1}{1+x^{2}}

\end{align}
証明はめんどくさいのでやりません。ところで、彼(↑)はマクローリン展開すると

 \displaystyle
\begin{align}

\frac{1}{1+x^{2}} = 1-x^2+x^4-...

\end{align}
こうなりますよね。証明はしません。めんどくさいので。
これらがわかると下の式が導けますね。

 \displaystyle
\begin{align}

\tan^{-1} x &= \int (\tan^{-1} x)^{'} dx \\ &= \int \frac{1}{1+x^{2}} dx
\\ &= \int (1-x^2+x^4-...)dx
\\ &= x - \frac{x^3}{3} + \frac{x^5}{5} - ...
\\ &= \lim_{n \to \infty} \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1}

\end{align}

おぉ、アークタンジェントがシグマで表せました!
ここで、


 \displaystyle
\begin{align}

f(x) = \lim_{n \to \infty} \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1}

\end{align}
とします。すると、Machinの公式より

 \displaystyle
\begin{align}

4f \left( \frac{1}{5} \right) - f \left( \frac{1}{239} \right) = \frac{\pi}{4} \\ \\

\therefore \pi = 16f \left( \frac{1}{5} \right) - 4f \left( \frac{1}{239} \right)

\end{align}

プログラム書いてみた

さて、本記事も終盤です。実際にプログラミングで近似します。
その前に、プログラムのイメージがしやすいように(?)上でやったf(x)をちょっと改変したg(x,n)を以下のように定義します。


 \displaystyle
\begin{align}

g(x, n) = \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1}

\end{align}
あんまり見た目は変わりませんが、プログラムの構造がわかりやすいと思ったので…
gを使ってπを表すと、以下のようになります。

 \displaystyle
\begin{align}

\pi = \lim_{n \to \infty} 16g \left( \frac{1}{5}, n \right) - 4g \left( \frac{1}{239}, n \right)

\end{align}
となりますね。つまり、gという関数を実装して、二つ連結してループで回す、という感じです。

以下、私が書いたRubyのコードです。

def maclaurin(x, n)
  sum = 0.0
  0.upto(n) do |k|
    sum += (-1.0)**k * x**(2.0 * k + 1.0) / (2.0 * k + 1.0)
  end
  sum
end

def machin(n)
  4 * maclaurin(1.0 / 5.0, n) - maclaurin(1.0 / 239.0, n)
end

gets.chomp.to_i.times do |i|
  puts "#{i}: #{4 * machin(i)}"
end

ちなみに、rubocopからは「引数名が短い!」って怒られてますw

破壊力ぅ...

どうしてMachinの公式がπの近似に使われるか。
それはこの公式の破壊力がすさまじいからです(語彙力)
なんとサンプル数(ここでいうn)を10にすると、もうRubyではこれ以上の近似はできません!
収束がめっちゃ早い。

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

終わりに

今回、英語の宿題を残して爆速で記事を書きましたが、久しぶりに楽しかったです。
またヒマになったら記事を上げていこうともうので、よかったらご覧になってください。
最後まで読んでいただき、ありがとうございました~。