Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created Dec 9, 2017
Embed
What would you like to do?
C++ source code to compute Pi with arctan's formula.
#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; // 配列サイズ
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++) {
// 各項の分数計算 ( 1つ前の計算結果に追加で除算・乗算 )
for (int i = 0; i < ks; i++) {
// 本来 x * x で除算したい(した方が高速だ)が、
// x が 8桁以下でも x * x が8桁超になるケースがあるため2回に分けている。
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桁超になるケースがあるため2回に分けている。
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");
}
#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 // 1つの int で8桁扱う
/*
* 計算クラス
*/
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
/***********************************************************
* 円周率計算 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 <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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment