誰もが小学生の時に友達とこぞって暗記する数学の超有名定数こと円周率ですが、みなさんは円周率を数値計算したいと思ったことは無いですか？私は日ごろからしてみたいと思っていました（？）。
パソコンで数値計算をしようとしたら、超有名ソフト「SUPER PI」や、現在でも速度面で覇権をとっている「y-cruncher」などが有名ですね。当然ながらそういったソフトに勝てるようなソフトは作れないのですが、円周率の数値計算の世界を少し覗いてみましょう。
逆正接関数（arctan）のテイラー展開
数値計算といったらテイラー展開です。過去にネイピア数を求めたときもテイラー展開でしたね。円周率の数値計算の話をすると、arctan（tanの逆関数）やarcsin（sinの逆関数）を使う方法が2番目くらいに出てくるかと思います（1番目は内接多角形・外接多角形の話かなと思います。個人の意見です）。というわけで、早速arctanのテイラー展開に取り組んでいきましょう。
まず、唐突ですがarctanの微分を考えます。
\[y=\arctan x\]
と置くと、定義より
\[x=\tan y\]
両辺微分して
\[\frac{dx}{dy}=\frac{d}{dy}\tan y =\frac{1}{\cos^2 y}\]
よって
\[\frac{dy}{dx}=\cos^2 y=\frac{1}{1+\tan^2 y}\]
ここで$y=\arctan x$を代入すると
\[\frac{dy}{dx}=(\arctan x)'=\frac{1}{1+x^2}\]
arctanの微分って意外とシンプルな表現になるんですね。両辺積分して、arctanの形で表現するようにしておきましょう。
\[\arctan x=\int_{0}^{x} \frac{1}{1+t^2} dt\]
ここで、さらに唐突ですが、無限等比級数の公式を思い出します。
\[ 1+r+r^2+r^3+\cdots=\frac{1}{1-r} \qquad (|r|<1) \]
これに$r=-t^2$を代入すると
\[ 1-t^2+t^4-t^6+\cdots=\frac{1}{1+t^2} \qquad (t^2<1) \]
あら不思議、右辺は先ほどのarctanの積分表記と同じですね。ということで代入して、
\[ \begin{align} \arctan x &=\int_{0}^{x} \frac{1}{1+t^2} dt \notag \\ &=\int_{0}^{x} (1-t^2+t^4-t^6+\cdots) dt \notag \\ &=\left [t-\frac{t^3}{3}+\frac{t^5}{5}-\frac{t^7}{7}+\cdots \right ]_{0}^{x} \notag \\ &= x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\cdots \notag \end{align} \]
無事テイラー展開ができました。
ついでに円周率を求めておきましょう。arctanの定義より、
\[ \arctan 1 = \frac{\pi}{4} \]
となるので、テイラー展開を適用して、
\[ \frac{\pi}{4} = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots \]
という級数が得られます（途中式で$t^2<1$という条件があったから$x=1$を入れられないじゃないかって？収束半径とかごにょごにょ改めて計算すると大丈夫らしいです）。見ての通り、この級数を計算して4倍すれば円周率ですね。
この級数はとてつもなく収束が遅いことで知られています。10桁計算するのに100億項まで計算する必要があるんだとか。そりゃ、分母が等差数列でしか増えていかなくて、+と-が交互に出てくるんじゃ収束は遅いですよね。
マチンの公式
1706年にイギリスの天文学者マチンは、次の公式を発表しました。
\[ \frac{\pi}{4} = 4 \arctan\frac{1}{5}-\arctan\frac{1}{239} \]
これを上のテイラー展開の式に当てはめて考えてみましょう。
\[ \frac{\pi}{4} = 4 \left (\frac{1}{5}-\frac{1}{3 \cdot 5^3}+\frac{1}{5 \cdot 5^5}-\frac{1}{7 \cdot 5^7}+\cdots \right) - \left (\frac{1}{239}-\frac{1}{3 \cdot 239^3}+\frac{1}{5 \cdot 239^5}-\frac{1}{7 \cdot 239^7}+\cdots \right ) \]
$x^{2n}$で増えていく項が分母に付いたので、収束がぐんと早くなることは想像に難くありません。実際、この方法だと10桁を計算するのに10項ほどで事足ります。
せっかくなので、この公式も導出しておきましょう。
$\theta=\arctan(1/5)$と置いて、2倍角の公式を2回使って$\tan 4\theta$を求めます。
\[ \tan 2\theta = \frac{2 \tan \theta}{1-\tan^2 \theta} = \frac{2 \cdot \frac{1}{5}}{1-\frac{1}{25}}=\frac{5}{12} \]
\[ \tan 4\theta = \frac{2 \tan 2\theta}{1-\tan^2 2\theta} = \frac{2 \cdot \frac{5}{12}}{1-\left(\frac{5}{12}\right)^2}=\frac{120}{119} \]
ここで、本命となる$\tan(4\theta-\pi/4)$を計算します。加法定理より
\[ \tan \left(4\theta-\frac{\pi}{4}\right) = \frac{\tan 4\theta - \tan\frac{\pi}{4}}{1+\tan 4\theta \tan\frac{\pi}{4}}=\frac{\frac{120}{119}-1}{1+\frac{120}{119}}=\frac{1}{239}\]
よって
\[ 4\theta-\frac{\pi}{4} = \arctan \frac{1}{239} \]
\[ \frac{\pi}{4} = 4\arctan\frac{1}{5}-\arctan \frac{1}{239} \]
実装
さて、マチンの公式を実装していきましょう。
static void Main()
{
int digits = 100;
int scale = digits + 10; // 十分な精度確保のため余裕を持たせる
var pi = 16 * Arccot(5, scale) - 4 * Arccot(239, scale);
var piText = ToDecimal(pi, scale).Substring(0, digits + 2);
Console.WriteLine(piText);
}
static BigInteger Arccot(int x, int scale)
{
BigInteger unity = BigInteger.Pow(10, scale);
BigInteger xPower = unity / x;
BigInteger sum = xPower;
BigInteger term;
int n = 3;
bool negative = true;
long x2 = (long)x * x;
while(true) {
xPower /= x2;
term = xPower / n;
if(term == 0)
break;
sum += negative ? -term : term;
negative = !negative;
n += 2;
}
return sum;
}
static string ToDecimal(BigInteger value, int scale)
{
string s = BigInteger.Abs(value).ToString().PadLeft(scale + 1, '0');
int decimalPosition = s.Length - scale;
return (value.Sign < 0 ? "-" : "") + s.Insert(decimalPosition, ".");
}
Arccot関数はいわゆる逆余接関数（コタンジェントの逆関数）です。要するに、正接の逆数を与えることになります。上のテイラー展開の式を見てもわかる通り、マチンの公式だと各項の分子は常に1ですので、コタンジェントと見立てれば引数はすべて整数にできるのですね。
また、計算の過程でも、例えば小数点以下第100桁まで求めたければ、ひとまず10^100を用意して、それを各項の分母で割って足していくことで級数の和を計算します。最終的に項の値が0になってしまったら、それ以降は1未満の項で無視できるということで計算を終了させています。これによって、分数を使わずにすべて整数として計算できるようになっています。
また、求めたい桁数に対して、10桁ほど多めに計算するようにしています。これは、テイラー展開の各項の計算における割り算の誤差や、最後に足し合わせた時の繰り上げ/繰り下げ等で多少の誤差が出ることから、少し多めの桁数を計算しているといったところです。本当はちゃんと挟みうちをして誤差を評価すべきでしょうが、まあ、除算や足し算の誤差なんて知れているので、この程度で概ね充分でしょう。
マチンの公式の拡張
さて、上のコードは充分にシンプルで良いものですが、例えば私のパソコンで実行すると、30万桁で約47秒かかります。もう少し良い方法は無いでしょうか？
実は、マチンの公式に類似した公式、すなわちarctanの中身が1/nであるものを足し合わせただけの公式は多数あります。当然分母が大きいほうが収束が早いはずですので、そのようなものがあればより早く計算できることが期待できます。
というわけで次の公式を使って計算してみましょう。
\[ \frac{\pi}{4}=44 \arctan \frac{1}{109}+95 \arctan \frac{1}{239}-12 \arctan \frac{1}{682}+24 \arctan \frac{1}{12943}-44 \arctan \frac{1}{6826318}\]
すべての分母が100超えということで、計算速度に期待できますね。
また、項が5つありますが、各項は独立して計算できるので並列処理してしまいましょう。
static void Main()
{
int digits = 100;
int scale = digits + 10; // 十分な精度確保のため余裕を持たせる
var members = new BigInteger[5];
Parallel.Invoke(
() => members[0] = 176 * Arccot(109, scale),
() => members[1] = 380 * Arccot(239, scale),
() => members[2] = -48 * Arccot(682, scale),
() => members[3] = 96 * Arccot(12943, scale),
() => members[4] = -176 * Arccot(6826318, scale)
);
var pi = members[0] + members[1] + members[2] + members[3] + members[4];
var piText = ToDecimal(pi, scale).Substring(0, digits + 2);
Console.WriteLine(piText);
}
これで30万桁計算すると、約15秒で計算が終わりました。マチンの公式に比べて3倍ほど早いですね。
項数が増えると各項を足し合わせるのに時間がかかるのでは？と思うかもしれませんが、多倍長整数の和は簡単に計算できるのでここはそんなにボトルネックになりません。arccotの計算にかかるが支配的で、その次は10進変換です。
まとめ
ということで、マチンの公式及びその亜種を使っての円周率の計算をやってきました。高校数学くらいで理解できる式で、簡単かつそれなりに計算できるものなのですね。
ちなみに、私のコードで30万桁15秒でしたが、SUPER PIだと52万桁で3秒でした。y-cruncherだと1億桁計算して5秒なので、足元にも及ばないですね…。
