C++ - 円周率計算(Arctan 系公式)!

Updated:


今まで、円周率をマチンの公式や Klingensitierna の公式、オイラーの公式で多桁計算する概念、C++ アルゴリズムを紹介しました。

Arctan 系の公式はこれら以外に多数(理論上無数に)ありますが、今回は4つまとめて紹介します。
「Gauß の公式」、「Störmer の公式」2つ、「高野喜久雄の公式」です。

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

また、Arctan 系公式は項数や係数が異るだけなので、1つのプログラムで公式を選択できるようにしています。

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2
  • 多桁計算で使用する1つの配列のサイズは8桁としている。
    (当方の環境で扱える int 型の範囲は -2,147,483,648 〜 2,147,483,647 であることから)
  • 指定する桁数は int 型の範囲としているが、あまり大きいと計算に膨大な時間を要するので注意!

1. 公式について

他の Arctan 系の公式と考え方は同じなので、詳細は前述の過去記事を参照のこと。
過去に紹介したものも含めて Arctan 系公式をまとめてみた。

PI_ARCTAN_1

ちなみに、\(\tan ^{-1} = \arctan\) と置き換えてもよい。

2. C++ ソース作成

例として、以下のようにソースを作成した。

  • 使用する公式番号を入力するようにしている。
  • 計算桁数を入力するようにしている。
  • 計算する項数もプログラム内で計算させるようにしている。
  • 計算結果をテキストファイルに出力するようにしている。
  • ソースは、メイン部分・計算クラス部分・計算クラスプロトタイプ宣言部分の3つに分けるようにした。
  • これまでと異るのは、\(x ^ 2\) が8桁を超えるケースが出てくるので、多桁計算ロジックの都合上、\(x\) での除算を2回行うように変更している。(このため、若干処理に時間がかかるようになってしまったが)

File: calc.h

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
#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

File: calc.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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
#include <math.h>
#include <string.h>
#include "Calc.h"

using namespace std;

/*
 * コンストラクタ
 */
Calc::Calc(int x, int y)
{
    // ==== 各種定数定義
    const char *cst_FORMULA[] = {
        "Machin", "Klingenstierna", "Euler",    "Euler2",
        "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: Euler
        { 16, 1,  5,  -4, 1,  70,   4, 1,  99,  0, 0,      0},  // 4: Euler(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");
}

File: calc_pi_arctan.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
/***********************************************************
 * 円周率計算 by Arctan 系公式
 *
 *    1: Machin
 *    2: Klingenstierna
 *    3: Euler
 *    4: Euler(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:Euler, 4:Euler(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;
}

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

$ g++ calc.h calc.cpp calc_pi_arctan.cpp -o calc_pi_arctan

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

4. 実行

実際に実行して検証してみる。
以下では、「高野喜久雄の公式」を使用して小数点以下 10,000 桁を計算している。

$ ./calc_pi_arctan
1:Machin, 2:Klingenstierna, 3:Euler, 4:Euler(2),
5:Gauss, 6:Stormer, 7:Stormer(2), 8:Takano : 8
Please input number of Pi Decimal-Digits : 10000

[ Output file : PI_Takano.txt ]

** Pi Computation by Arctan method **
   Formula = [ Takano ]
   Digits  = 10000
   Time    = 1.180000 seconds

公式別にテキストファイルができているはずである。

File: PI_Takano.txt

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
** Pi Computation by Arctan method **
   Formula = [ Takano ]
   Digits  = 10000
   Time    = 1.180000 seconds

          3.
00000001: 14159265 35897932 38462643 38327950 28841971 69399375 10582097 49445923 07816406 28620899
00000081: 86280348 25342117 06798214 80865132 82306647 09384460 95505822 31725359 40812848 11174502
00000161: 84102701 93852110 55596446 22948954 93038196 44288109 75665933 44612847 56482337 86783165
00000241: 27120190 91456485 66923460 34861045 43266482 13393607 26024914 12737245 87006606 31558817
00000321: 48815209 20962829 25409171 53643678 92590360 01133053 05488204 66521384 14695194 15116094
00000401: 33057270 36575959 19530921 86117381 93261179 31051185 48074462 37996274 95673518 85752724
00000481: 89122793 81830119 49129833 67336244 06566430 86021394 94639522 47371907 02179860 94370277
00000561: 05392171 76293176 75238467 48184676 69405132 00056812 71452635 60827785 77134275 77896091
00000641: 73637178 72146844 09012249 53430146 54958537 10507922 79689258 92354201 99561121 29021960
00000721: 86403441 81598136 29774771 30996051 87072113 49999998 37297804 99510597 31732816 09631859
00000801: 50244594 55346908 30264252 23082533 44685035 26193118 81710100 03137838 75288658 75332083
00000881: 81420617 17766914 73035982 53490428 75546873 11595628 63882353 78759375 19577818 57780532
00000961: 17122680 66130019 27876611 19590921 64201989 38095257 20106548 58632788 65936153 38182796
        :
===< 途中省略 >===
        :
00009121: 37800297 41207665 14793942 59029896 95946995 56576121 86561967 33786236 25612521 63208628
00009201: 69222103 27488921 86543648 02296780 70576561 51446320 46927906 82120738 83778142 33562823
00009281: 60896320 80682224 68012248 26117718 58963814 09183903 67367222 08883215 13755600 37279839
00009361: 40041529 70028783 07667094 44745601 34556417 25437090 69793961 22571429 89467154 35784687
00009441: 88614445 81231459 35719849 22528471 60504922 12424701 41214780 57345510 50080190 86996033
00009521: 02763478 70810817 54501193 07141223 39086639 38339529 42578690 50764310 06383519 83438934
00009601: 15961318 54347546 49556978 10382930 97164651 43840700 70736041 12373599 84345225 16105070
00009681: 27056235 26601276 48483084 07611830 13052793 20542746 28654036 03674532 86510570 65874882
00009761: 25698157 93678976 69742205 75059683 44086973 50201410 20672358 50200724 52256326 51341055
00009841: 92401902 74216248 43914035 99895353 94590944 07046912 09140938 70012645 60016237 42880210
00009921: 92764579 31065792 29552498 87275846 10126483 69998922 56959688 15920560 01016552 56375678

全ての公式による計算結果が一致することを確認する。

5. 検証

序盤で紹介している8つの公式を使用して、どれがどの程度高速に計算できているのかを検証してみた。
小数点以下 1,000 桁、 10,000 桁、 100,000 桁を5回計算して平均をとってみた。
ちなみに計算に使用したマシンのスペックは以下のとおり。

  • CPU: Core2Duo E8500 (3.16GHz)
  • メモリ: 3.9 GiB

PI_ARCTAN_2

今回の8つの公式の中では、「Störmer の公式(2)」による計算が最も高速で、次が「高野喜久雄の公式」による計算であった。
一般的に発表時期の新しい公式の方が高速であるように感じた。 また、計算する桁数を 10 倍すれば、計算時間は 100 倍になる傾向にある。(そういう計算式だから)

6. 参考サイト


自分の環境と相談して、計算する桁数を大きくしてみるのもよいでしょう。

ちなみに、結果が何万桁以上になる場合は、多桁(多倍長)乗算に上記の方法ではなく「Karatsuba法」や「Toom-Cook法」や「FFT(高速フーリエ変換)」を使うのが一般的なようです。
これも踏まえて、多桁(多倍長)演算を改良してより高速化に挑んでみたいとも思っている次第です。

以上。





 

Sponsored Link

 

Comments