C++ - 円周率計算（マチンの公式）！

先日は、コンピュータで大きな桁数を計算する概念・アルゴリズムを紹介しました。

• C++ - 多桁計算！

今回は、その概念を応用して円周率 \(\pi\) を計算してみました。

使用した公式は、マチン(Machin)の公式です。
当然、プログラミン言語そのものが保有している三角関数は使用しません。級数展開して計算します。（多桁の円周率計算では、全く無意味ですから）

0. 前提条件

• Linux Mint 14 Nadia (64bit) での作業を想定。
• g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2
• 多桁計算で使用する1つの配列のサイズは8桁としている。
  （当方の環境で扱える int 型の範囲は -2,147,483,648 〜 2,147,483,647 であることから）

1. マチン(Machin)の公式について

数式が多いので \(\TeX\) で記載。

ちなみに、\(\tan ^{-1} = \arctan\) と置き換えてもよい。

なお、ここではマチンの公式の証明はしない。（証明方法は何種類かあります）

また、上記の方法で算出した計算項数では、収束の速い方は無駄に多く計算していることになる。
それを防ぎたかったら、各 Arctan 毎に必要な項数分だけ計算して、最後に合算する方法を取ると良いだろう。

2. C++ ソース作成

例として、以下のようにソースを作成した。
計算する桁数は、小数点以下 1,000 桁としている。
ソース内の L の値を変更すれば、任意の桁数を計算可能。

File: calc_pi_machin.cpp

/*********************************************
 * 円周率計算 by マチンの公式
 *********************************************/
#include <iostream>  // for cout
#include <stdio.h>   // for printf()

#define L  1000                 // 算出桁数
#define L1 ((L / 8) + 1)        // 配列サイズ
#define L2 (L1 + 1)             // 配列サイズ + 1
#define N  ((L / 1.39794) + 1)  // 計算項数

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 各種変数
    int k, i;           // LOOP インデックス
    int carry, borrow;  // 繰り上がり、借り
    long w;             // 被乗数、被除数ワーク
    long r;             // 剰余ワーク

    public:
        // 計算
        void calc();
        // ロング + ロング
        void ladd(int *, int *, int *);
        // ロング - ロング
        void lsub(int *, int *, int *);
        // ロング / ショート
        void ldiv(int *, int, int *);
        // 結果出力
        void display(int *);
};

/*
 * 計算
 */
void Calc::calc()
{
    // 配列宣言(総和、マシンの公式の前の項、後の項、前の項＋後の項)
    int s[L2 + 1], a[L2 + 1], b[L2 + 1], q[L2 + 1];

    // 配列初期化
    for (k = 0; k <= L2; k++)
        s[k] = a[k] = b[k] = q[k] = 0;

    // マチンの公式
    a[0] = 16 * 5;
    b[0] = 4 * 239;
    for (k = 1; k <= N; k++) {
        ldiv(a, 5 * 5, a);
        ldiv(b, 239 * 239, b);
        lsub(a, b, q);
        ldiv(q, 2 * k - 1, q);
        if ((k % 2) != 0)
            ladd(s, q, s);
        else
            lsub(s, q, s);
    }

    // 結果出力
    display(s);
}

/*
 * ロング + ロング
 */
void Calc::ladd(int a[], int b[], int c[])
{
    carry = 0;
    for (i = L2; i >=0; i--) {
        c[i] = a[i] + b[i] + carry;
        if (c[i] < 100000000) {
            carry = 0;
        } else {
            c[i] -= 100000000;
            carry = 1;
        }
    }
}

/*
 * ロング - ロング
 */
void Calc::lsub(int a[], int b[], int c[])
{
    borrow = 0;
    for (i = L2; i >=0; i--) {
        c[i] = a[i] - b[i] - borrow;
        if (c[i] >= 0) {
            borrow = 0;
        } else {
            c[i] += 100000000;
            borrow = 1;
        }
    }
}

/*
 * ロング / ショート
 */
void Calc::ldiv(int d[], int e, int f[])
{
    r = 0;
    for (i = 0; i < L2 + 1; i++) {
        w = d[i];
        f[i] = (w + r) / e;
        r = ((w + r) % e) * 100000000;
    }
}

/*
 * 結果出力
 */
void Calc::display(int s[])
{
    printf("%7d. ", s[0]);
    for (i = 1; i < L1; i++)
        printf("%08d ", s[i]);
    printf("\n");
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 円周率計算
        objCalc.calc();
    }
    catch (...) {
        cout << "例外発生！" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

• Gist - C++ source code to compute Pi with Machin’s formula.

3. C++ ソースコンパイル

$ g++ calc_pi_machin.cpp -o calc_pi_machin

何も出力されなければ成功です。

4. 実行

$ ./_calc_pi_machin
      3. 14159265 35897932 38462643 38327950 28841971 69399375 10582097 
49445923 07816406 28620899 86280348 25342117 06798214 80865132 82306647 
09384460 95505822 31725359 40812848 11174502 84102701 93852110 55596446 
22948954 93038196 44288109 75665933 44612847 56482337 86783165 27120190 
91456485 66923460 34861045 43266482 13393607 26024914 12737245 87006606 
31558817 48815209 20962829 25409171 53643678 92590360 01133053 05488204 
66521384 14695194 15116094 33057270 36575959 19530921 86117381 93261179 
31051185 48074462 37996274 95673518 85752724 89122793 81830119 49129833 
67336244 06566430 86021394 94639522 47371907 02179860 94370277 05392171 
76293176 75238467 48184676 69405132 00056812 71452635 60827785 77134275 
77896091 73637178 72146844 09012249 53430146 54958537 10507922 79689258 
92354201 99561121 29021960 86403441 81598136 29774771 30996051 87072113 
49999998 37297804 99510597 31732816 09631859 50244594 55346908 30264252 
23082533 44685035 26193118 81710100 03137838 75288658 75332083 81420617 
17766914 73035982 53490428 75546873 11595628 63882353 78759375 19577818 
57780532 17122680 66130019 27876611 19590921 64201989 

当方の非力なマシンでも１秒かかりませんでした。
自分の環境と相談して、計算する桁数を大きくしてみるのもよいでしょう。
（当方の非力なマシンで10万桁を計算してみたら、15秒くらいだった）

但し、結果が何万桁以上になる場合は、多桁（多倍長）乗算に上記の方法ではなく「Karatsuba法」や「Toom-Cook法」や「FFT（高速フーリエ変換）」を使うのが一般的なようだ。

5. 参考サイト

• 円周率を求める公式・アルゴリズム - 円周率.jp -
• πのArctan公式集

---

上記の参考サイトにあるように、円周率の計算は「マチンの公式」を含めた数々の arctan 系の公式の他に、Ramanujan 系、連分数系、AGM 系、Borwein 系、BBP 系の公式を利用する方法や、多角形を利用した方法もあります。

他の方法で試してみるのもよいですし、より計算速度を上げる方法を考えてみるのよいでしょう。

ちなみに、当方はより高速な「チャドノフスキー(Chudnovsky)の公式」に興味があります。

※ちなみに最近の当方の C++ アルゴリズムについての記事は、古い C によるアルゴリズムに関する書物を参考に C++ に移植した形態となっています。

以上。