C++ source code to compute Pi with arctan's formula.(old ver.)

calc.cpp

#include <math.h>
#include <string.h>
#include "Calc.h"

using namespace std;

/*
 * コンストラクタ
 */
Calc::Calc(int x, int y)
{
    // ==== 各種定数定義
    const char *cst_FORMULA[] = {              // 公式名
        "Machin",                              //   1: Machin
        "Klingenstierna",                      //   2: Klingenstierna
        "Euler",                               //   3: Euler
        "Euler2"                               //   4: Euler(2)
    };
    const char *str_pre = "PI_";               // 保存ファイル名・プリフィックス
    const char *str_ext = ".txt";              // 保存ファイル名・拡張子
    const int cst_KS[] = {2, 3, 2, 3};         // 公式の項数
    const int cst_KK[4][12] = {                // 公式内係数
        {16, 1,  5, -4, 1, 239,   0, 0,   0, 0, 0, 0}, // 1: Machin
        {32, 1, 10, -4, 1, 239, -16, 1, 515, 0, 0, 0}, // 2: Klingenstierna
        {20, 1,  7,  8, 3,  79,   0, 0,   0, 0, 0, 0}, // 3: Euler
        {16, 1,  5, -4, 1,  70,   4, 1,  99, 0, 0, 0}, // 4: Euler(2)
    };

    // ==== 各種変数設定
    formula = (char *)cst_FORMULA[x-1];        // 公式名取得
    // 結果出力ファイル名生成
    fname[0] = '\0';
    strcat(fname, str_pre);
    strcat(fname, formula);
    strcat(fname, str_ext);
    printf("\n[ Output file : %s ]\n\n",fname);
    ks = cst_KS[x-1];                          // 計算対象公式の項数
    for (int i = 0; i < ks * 3; i++)
        kk[i] = cst_KK[x-1][i];                // 計算対象公式の係数
    l  = y;                                    // 計算桁数
    l1 = (l / 8) + 1;                          // 配列サイズ
    n  = (l / log10(kk[2]) + 1) / 2 + 1;       // 計算項数
}

/*
 * 計算
 */
void Calc::calc()
{
    // ==== 計算開始時刻
    t1 = clock();

    // ==== 配列宣言
    int s[l1 + 2], a[ks][l1 + 2], q[l1 + 2];

    // ==== 配列初期化
    for (int i = 0; i <= l1 + 1; i++) {
        s[i] = q[i] = 0;
        for (int j = 0; j < ks; j++)
            a[j][i] = 0;
    }

    // ==== 計算初期値
    for (int i = 0; i < ks; i++) {
        if (kk[i * 3] >= 0) {
            a[i][0] = kk[i * 3] * kk[i * 3 + 2];
        } else {
            a[i][0] = kk[i * 3] * kk[i * 3 + 2] * -1;
        }
        // 分子が 1 で無ければ、さらに分子の値で除算しておく
        if (kk[i * 3 + 1] != 1)
            ldiv(a[i], kk[i * 3 + 1], a[i]);
    }

    // ==== 計算
    for (int k = 1; k <= n; k++) {
        // 各項の分数計算 ( １つ前の計算結果に追加で除算・乗算 )
        for (int i = 0; i < ks; i++) {
            ldiv(a[i], kk[i * 3 + 2] * kk[i * 3 + 2], a[i]);
            if (kk[i * 3 + 1] != 1)  // 分子が 1 でない時
                lmul(a[i], kk[i * 3 + 1] * kk[i * 3 + 1], a[i]);
        }

        // 各項の加減算
        for (int i = 1; i < ks; i++) {
            if (kk[i * 3] >= 0) {  // 加算
                if (i == 1)
                    ladd(a[i - 1], a[i], q);
                else
                    ladd(q, a[i], q);
            } else {               // 減算
                if (i == 1)
                    lsub(a[i - 1], a[i], q);
                else
                    lsub(q, a[i], q);
            }
        }

        // 1 / (2 * k - 1) の部分
        ldiv(q, 2 * k - 1, q);

        // 総和計算
        if ((k % 2) != 0)
            ladd(s, q, s);
        else
            lsub(s, q, s);
    }

    // ==== 計算終了時刻
    t2 = clock();

    // ==== 計算時間
    tt = (double)(t2 - t1) / CLOCKS_PER_SEC;

    // ==== 結果出力
    display(tt, s);
}

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

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

/*
 * ロング * ショート
 */
void Calc::lmul(int d[], int e, int f[])
{
    cr = 0;
    for (int i = l1 + 1; i >=0; i--) {
        w = d[i];
        f[i] = (w * e + cr) % MAX_DIGITS;
        cr = (w * e + cr) / MAX_DIGITS;
    }
}

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

/*
 * 結果出力
 */
void Calc::display(double tt, int s[])
{
    // ==== コンソール出力
    printf(STR_TITLE);
    printf(STR_FORMULA, formula);
    printf(STR_DIGITS, l);
    printf(STR_TIME, tt);

    // ==== ファイル出力
    out_file = fopen(fname, "w");
    fprintf(out_file, STR_TITLE);
    fprintf(out_file, STR_FORMULA, formula);
    fprintf(out_file, STR_DIGITS, l);
    fprintf(out_file, STR_TIME, tt);
    fprintf(out_file, "\n          %d.\n", s[0]);
    for (int i = 1; i < l1; i++) {
        if (i % 10 == 1) fprintf(out_file, "%08d:", (i - 1) * 8 + 1);
        fprintf(out_file, " %08d", s[i]);
        if (i % 10 == 0) fprintf(out_file, "\n");
    }
    fprintf(out_file, "\n\n");
}

calc.h

#ifndef INCLUDED_CALC_H
#define INCLUDED_CALC_H

#include <stdio.h>
#include <time.h>

// 出力文字列 ( コンソール、テキストファイル共通 )
#define STR_TITLE   "** Pi Computation by Arctan method **\n"
#define STR_FORMULA "   Formula = [ %s ]\n"
#define STR_DIGITS  "   Digits  = %d\n"
#define STR_TIME    "   Time    = %f seconds\n"
// 多桁計算
#define MAX_DIGITS 100000000  // １つの int で８桁扱う

/*
 * 計算クラス
 */
class Calc
{
    // 各種変数
    int ks;           // 公式・項数
    int kk[12];       // 公式・係数
    int l, l1, n;     // 計算桁数、配列サイズ、計算項数
    int cr, br;       // 多桁計算・繰り上がり、借り
    long w;           // 多桁計算・被乗数、被除数ワーク
    long r;           // 多桁計算・剰余ワーク
    clock_t t1, t2;   // 計算開始CPU時刻、計算終了CPU時刻
    double tt;        // 計算時間
    char *formula;    // 公式名
    char *str_pre;    // 結果出力ファイル名・プリフィックス
    char *str_ext;    // 結果出力ファイル名・拡張子
    char fname[32];   // 結果出力ファイル名
    FILE *out_file;   // 結果出力ファイル名・ポインタ

    public:
        Calc(int, int);                  // コンストラクタ
        void calc();                     // 計算
        void ladd(int *, int *, int *);  // ロング + ロング
        void lsub(int *, int *, int *);  // ロング - ロング
        void lmul(int *, int,   int *);  // ロング * ショート
        void ldiv(int *, int,   int *);  // ロング / ショート
        void display(double, int *);     // 結果出力
};

#endif

calc_pi_arctan_old.cpp

/***********************************************************
 * 円周率計算 by Arctan 系公式
 *
 *    1: Machin
 *    2: Klingenstierna
 *    3: Euler
 *    4: Euler(2)
 *
 * コンパイル方法：
 *   $ g++ Calc.h Calc.cpp CalcPiArctan.cpp -o CalcPiArctan
 **********************************************************/
#include <iostream>
#include "Calc.h"

using namespace std;

/*
 * メイン処理
 */
int main()
{
    int f, n;  // 使用公式番号、計算桁数

    try
    {
        // ==== 使用公式番号入力
        printf("1:Machin, 2:Klingenstierna, 3:Euler, 4:Euler(2) : ");
        scanf("%d", &f);
        if(f < 1 || f > 4) f = 1;  // 範囲外なら 1:Machin と判断

        // ==== 計算桁数入力
        printf("Please input number of Pi Decimal-Digits : ");
        scanf("%d", &n);

        // ==== 計算クラスインスタンス化
        Calc objCalc(f, n);

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

    // ==== 正常終了
    return 0;
}