C++ source code to compute Pi with Machin's formula.
/********************************************* | |
* 円周率計算 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; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment