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

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", "Klingenstierna", "Eular",    "Eular2",
        "Gauss",  "Stormer",        "Stormer2", "Takano"
    };                             // 公式名
    const char *str_pre = "PI_";   // 保存ファイル名・プリフィックス
    const char *str_ext = ".txt";  // 保存ファイル名・拡張子
    const int cst_KS[] = {2, 3, 2, 3, 3, 3, 4, 4};              // 公式の項数
    const int cst_KK[8][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: Eular
        { 16, 1,  5,  -4, 1,  70,   4, 1,  99,  0, 0,      0},  // 4: Eular(2)
        { 48, 1, 18,  32, 1,  57, -20, 1, 239,  0, 0,      0},  // 5: Gauss
        { 24, 1,  8,   8, 1,  57,   4, 1, 239,  0, 0,      0},  // 6: Stormer
        {176, 1, 57,  28, 1, 239, -48, 1, 682, 96, 1,  12943},  // 7: Stormer(2)
        { 48, 1, 49, 128, 1,  57, -20, 1, 239, 48, 1, 110443}   // 8: Takano
    };

    // ==== 各種変数設定
    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;                     // 配列サイズ
    for (int i = 0; i < ks; i++)          // 計算項数
        n[i]  = (l / (log10(kk[i * 3 + 2]) - log10(kk[i * 3 + 1])) + 1) / 2 + 1;
}

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

    // ==== 配列宣言
    // 各項の(2k-1)部分を除いた部分、各項の(2k-1)部分含む部分、各項の総和、総和
    int a[ks][l1 + 2], q[ks][l1 + 2], sk[ks][l1 + 2], s[l1 + 2];

    // ==== 配列初期化
    for (int i = 0; i < l1 + 2; i++) {
        s[i] = 0;
        for (int j = 0; j < ks; j++)
            a[j][i] = sk[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 i = 0; i < ks; i++) {
        for (int k = 1; k <= n[i]; k++) {
            // 各項の分数の部分
            ldiv(a[i], kk[i * 3 + 2], a[i]);
            ldiv(a[i], kk[i * 3 + 2], a[i]);
            if (kk[i * 3 + 1] != 1) {  // 分子が 1 でない時
                lmul(a[i], kk[i * 3 + 1], a[i]);
                lmul(a[i], kk[i * 3 + 1], a[i]);
            }

            // 各項の 1 / (2 * k - 1) の部分
            ldiv(a[i], 2 * k - 1, q[i]);

            // 各項の総和に加減算
            if (k % 2 != 0)
                ladd(sk[i], q[i], sk[i]);
            else
                lsub(sk[i], q[i], sk[i]);
        }
    }

    // 各項の総和の加減算
    for (int i = 0; i < ks; i++) {
        if (kk[i * 3] >= 0) {  // 加算
            ladd(s, sk[i], s);
        } else {               // 減算
            lsub(s, sk[i], 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;        // 計算桁数、配列サイズ
    int n[4];         // 計算項数
    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_2.cpp

/***********************************************************
 * 円周率計算 by Arctan 系公式
 * ( 各項の Arctan を個別に計算後に加減算する方法 )
 *
 *    1: Machin
 *    2: Klingenstierna
 *    3: Eular
 *    4: Eular(2)
 *    5: Gauss
 *    6: Stormer
 *    7: Stormer(2)
 *    8: Takano
 *
 * コンパイル方法：
 *   $ 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:Eular, 4:Eular(2),\n");
        printf("5:Gauss, 6:Stormer, 7:Stormer(2), 8:Takano : ");
        scanf("%d", &f);
        if(f < 1 || f > 8) 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;
}