mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

C++ - 円周率計算(by マチンの公式)!

[ プログラミング, 数学 ] [ C言語, 円周率 ]

こんばんは。

先日は、コンピュータで大きな桁数を計算する概念・アルゴリズムを紹介しました。

今回は、その概念を応用して円周率 を計算してみました。

使用した公式は、マチン(Machin)の公式です。
当然、プログラミン言語そのものが保有している三角関数は使用しません。級数展開して計算します。(多桁の円周率計算では、全く無意味ですから)

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2
  • 多桁計算で使用する1つの配列のサイズは8桁としている。
    (当方の環境で扱える int 型の範囲は -2147483648 〜 2147483647 であることから)

1. マチン(Machin)の公式について

数式が多いので で記載。

CALC_PI_MACHIN_1

CALC_PI_MACHIN_2

ちなみに、 と置き換えてもよい。

なお、ここではマチンの公式の証明はしない。(証明方法は何種類かあります)

また、上記の方法で算出した計算項数では、収束の速い方は無駄に多く計算していることになる。
それを防ぎたかったら、各 Arctan 毎に必要な項数分だけ計算して、最後に合算する方法を取ると良いだろう。

2. C++ ソース作成

例として、以下のようにソースを作成した。
計算する桁数は、小数点以下 1,000 桁としている。
ソース内の L の値を変更すれば、任意の桁数を計算可能。

CalcPiMachin.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
/*********************************************
 * 円周率計算 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;
}

3. C++ ソースコンパイル

以下のコマンドで C++ ソースをコンパイルする。

1
$ g++ CalcPiMachin.cpp -o CalcPiMachin

何も出力されなければ成功です。

4. 実行

実際に実行して検証してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
$ ./CalcPiMachin
      3. 14159265 35897932 38462643 38327950 28841971 69399375 10582097
49445923 07816406 28620899 86280348 25342117 06798214 80865132 82306647
09384460 95505822 31725359 40812848 11174502 84102701 93852110 55596446
22948954 93038196 44288109 75665933 44612847 56482337 86783165 27120190
91456485 66923460 34861045 43266482 13393607 26024914 12737245 87006606
31558817 48815209 20962829 25409171 53643678 92590360 01133053 05488204
66521384 14695194 15116094 33057270 36575959 19530921 86117381 93261179
31051185 48074462 37996274 95673518 85752724 89122793 81830119 49129833
67336244 06566430 86021394 94639522 47371907 02179860 94370277 05392171
76293176 75238467 48184676 69405132 00056812 71452635 60827785 77134275
77896091 73637178 72146844 09012249 53430146 54958537 10507922 79689258
92354201 99561121 29021960 86403441 81598136 29774771 30996051 87072113
49999998 37297804 99510597 31732816 09631859 50244594 55346908 30264252
23082533 44685035 26193118 81710100 03137838 75288658 75332083 81420617
17766914 73035982 53490428 75546873 11595628 63882353 78759375 19577818
57780532 17122680 66130019 27876611 19590921 64201989

当方の非力なマシンでも1秒かかりませんでした。
自分の環境と相談して、計算する桁数を大きくしてみるのもよいでしょう。
(当方の非力なマシンで10万桁を計算してみたら、15秒くらいだった)

但し、結果が何万桁以上になる場合は、多桁(多倍長)乗算に上記の方法ではなく「Karatsuba法」や「Toom-Cook法」や「FFT(高速フーリエ変換)」を使うのが一般的なようだ。

5. 参考サイト


上記の参考サイトにあるように、円周率の計算は「マチンの公式」を含めた数々の arctan 系の公式の他に、Ramanujan 系、連分数系、AGM 系、Borwein 系、BBP 系の公式を利用する方法や、多角形を利用した方法もあります。

他の方法で試してみるのもよいですし、より計算速度を上げる方法を考えてみるのよいでしょう。

ちなみに、当方はより高速な「チャドノフスキー(Chudnovsky)の公式」に興味があります。

※ちなみに最近の当方の C++ アルゴリズムについての記事は、古い C によるアルゴリズムに関する書物を参考に C++ に移植した形態となっています。

以上。

Comments