Instantly share code, notes, and snippets.

komasaru/calc.cpp

Created Dec 9, 2017
C++ source code to compute Pi with arctan's formula.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #include #include #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; // 配列サイズ 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++) { // 本来 x * x で除算したい（した方が高速だ）が、 // x が 8桁以下でも x * x が8桁超になるケースがあるため２回に分けている。 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 でない時 // 本来 x * x を乗算したい（した方が高速だ）が、 // x が 8桁以下でも x * x が8桁超になるケースがあるため２回に分けている。 lmul(a[i], kk[i * 3 + 1], a[i]); lmul(a[i], 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"); }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #ifndef INCLUDED_CALC_H #define INCLUDED_CALC_H #include #include // 出力文字列 ( コンソール、テキストファイル共通 ) #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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 /*********************************************************** * 円周率計算 by 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 #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; }