C++ source code to compute Pi with Machin's formula.

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;
}