Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.