mk-mode BLOG

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

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

C++ - 多倍長浮動小数点数の加減算!

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

こんばんは。

これまで、多倍長整数の演算については話をしていましたが、今回は多倍長の浮動小数点同士の加減算について考えてみました。

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2

1. 考え方

整数部を 0 とする A = 0.123456789E5, B = 0.123456789E3 という浮動小数点を以下のように考える。(1桁を1個の配列で扱う場合)

ADD_BIG_FLOAT_1

(べき指数部分を別の変数に持たせる考え方もある。その場合は、整数部を使用する考え方もできる)

そして、計算結果格納用配列を Z[0], Z[1], ・・・, Z[9] とする。

加算の場合、べき指数の大きい方に合わせる。
減算の場合は、被減数の方が大きいこと。(計算結果が負にならないこと)
また、べき指数の小さい方の端数は処理しない。すなわち、切り捨てとする。(もう一方の数字が有効数字範囲外で不明のため、計算結果がどうなるか判別不能だから)

以下は、加算の場合。(減算も同様)

  1. べき指数の大きい方を結果配列に格納。
    ADD_BIG_FLOAT_2
  2. 2数のべき指数の差の部分を結果配列に格納。(直接加算処理を行わない部分)
    ADD_BIG_FLOAT_3
  3. 2数の加算すべき桁の加算結果を配列に格納。
    ADD_BIG_FLOAT_4
  4. 正規化(繰り上がり処理)を行う。
    ADD_BIG_FLOAT_5
  5. 正規化(浮動小数点位置調整)を行う。
    ( 0.099999999E5 となった場合 0.999999990E4 とするということ。今回の例ではこの処理を行なっても変化は無い。
    また、小数点以下第1位は「4. 正規化(繰り上がり処理)」で繰り上がり処理を行なっていないので、この部分もここで処理。)
    ADD_BIG_FLOAT_6

2. C++ ソース作成

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

  • 被加減数・加減数の桁数は合わせている。
  • テストなので、被加減数・加減数は適当な乱数で生成している。
  • 今後汎用性を高めるための処置として、色々と部品化している。
  • ソースファイルは1つ。(もう少し大きくなるか複雑になれば、ファイルを分割するかもしれないが)
AddBigFloat.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
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
/*********************************************
 * 浮動小数点数加減算
 *********************************************/
#include <cstdlib>   // for rand()
#include <iostream>  // for cout
#include <math.h>    // for pow()
#include <stdio.h>   // for printf()
#include <time.h>    // for time()

#define D_MAX  11  // 有効桁数(小数点以下+べき指数)

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    int A[D_MAX];  // 被加減数配列
    int B[D_MAX];  // 加数配列
    int C[D_MAX];  // 減数配列

    public:
        Calc();                                    // コンストラクタ
        void calc();                               // 計算
    private:
        void add(int *, int *, int *, int, int);   // 加減算   ( 固定長整数 )
        void fadd(int *, int *, int *, int, int);  // 加減算   ( 浮動小数点数 )
        void setv(int, int*, int);                 // 値セット ( 全配列に同じ整数をセット )
        void copy(int *, int*, int);               // 値コピー ( 全配列をそのままコピー )
        void norm(int *, int);                     // 正規化   ( 固定長整数 )
        void fnorm(int *, int);                    // 正規化   ( 浮動小数点数 )
        void fdisp(int *, int);                    // 結果出力 ( 指数表示 )
};

/*
 * コンストラクタ
 */
Calc::Calc()
{
    // 使用する被加減数・加減数を設定(テストなので適当な乱数を使用)
    //   ( A + B, A - C で A > C となることが条件だが
    //     稀に条件に合致しない値が生成されるかもしれない )
    srand((unsigned int)time(NULL));                     // 乱数の種を初期化
    A[0] = rand() % (D_MAX - 1) + 1;                     // Aのべき指数部
    B[0] = rand() % (D_MAX - 1) + 1;                     // Bのべき指数部
    C[0] = A[0] - rand() % A[0] - 1;                     // Cのべき指数部
    for (int i = 1; i < D_MAX; i++) A[i] = rand() % 10;  // Aの小数部
    for (int i = 1; i < D_MAX; i++) B[i] = rand() % 10;  // Bの小数部
    for (int i = 1; i < D_MAX; i++) C[i] = rand() % 10;  // Bの小数部
}

/*
 * 計算
 */
void Calc::calc()
{
    // 各種宣言
    int size = D_MAX;  // 配列サイズ
    int a[size];       // 被加減数配列
    int b[size];       // 加数配列
    int c[size];       // 減数配列
    int z[size + 1];   // 計算結果配列

    // 初期設定
    copy(A, a, size);  // 被加減数配列
    copy(B, b, size);  // 加数配列
    copy(C, c, size);  // 減数配列

    // 計算する被加減数・加数・減数を出力
    printf("A =\n"); fdisp(a, size);
    printf("B =\n"); fdisp(b, size);
    printf("C =\n"); fdisp(c, size);

    // 加算 ( Z = A + B )
    fadd(a, b, z, size, 1);

    // 結果出力
    printf("A + B =\n"); fdisp(z, size);

    // 減算 ( Z = A - C )
    fadd(a, c, z, size, -1);

    // 結果出力
    printf("A - C =\n"); fdisp(z, size);
}

/*
 * 加減算
 * ( 固定長整数同士の加減算, 正規化含む )
 *
 *   a     [I ]  被加減数配列
 *   b     [I ]  加減数配列
 *   size  [I ]  配列数(A, B)
 *   z     [ O]  出力データ配列(Z = A + B or Z = A - B)
 *   id    [I ]  1: 加算(A + B(, -1: 減算(A - B)
 */
void Calc::add(int *a, int *b, int *z, int size, int id)
{
    int i;          // LOOPインデックス

    // 計算
    if (id >= 0) {  // 加算
        for (i = 0; i < size; i++)
            z[i] = a[i] + b[i];
    } else {        // 減算
        for (i = 0; i < size; i++)
            z[i] = a[i] - b[i];
    }

    // 正規化
    norm(z, size);
}

/*
 * 加減算
 * ( 浮動小数点数同士の加減算, 正規化含む )
 *
 *   a     [I ]  被加減数配列
 *   b     [I ]  加減数配列
 *   size  [I ]  配列数(A, B)
 *   z     [ O]  出力データ配列(Z = A + B or Z = A - B)
 *   id    [I ]  1: 加算(A + B), -1: 減算(A - B)
 */
void Calc::fadd(int *a, int *b, int *z, int size, int id)
{
    int k;  // べき指数の差
    int i;  // LOOPインデックス

    // 計算
    if (a[0] >= b[0]) {  // A のべき指数 >= B のべき指数
        k = a[0] - b[0];
        copy(a, z, k + 1);
        add(&a[k + 1], &b[1], &z[k + 1], size - k - 1, id);
    } else {             // A のべき指数 <  B のべき指数
        z[0] = b[0];
        k = b[0] - a[0];
        // 位の異なる部分
        if (id >= 0) {   // 加算
            copy(&b[1], &z[1], k);
        } else {         // 減算
            for (i = 1; i <= k; i++) z[i] = -b[i];
        }
        // 位の同じ部分の加減算
        add(&a[1], &b[k + 1], &z[k + 1], size - k - 1, id);
    }

    // 正規化
    norm(&z[1], size - 1);  // 固定長整数
    fnorm(z, size);         // 浮動小数点数
}

/*
 * 値セット
 * ( 全配列に同じ整数をセット )
 *
 *   v     [I  ]  セットする整数
 *   z     [  O]  セット先配列
 *   size  [I  ]  配列サイズ
 */
void Calc::setv(int v, int *z, int size) {
    int i;  // LOOPインデックス

    for (i = 0; i < size; i++) z[i] = v;
}

/*
 * 値コピー
 * ( 全配列をコピー )
 *
 *   a     [I  ]  コピー元配列
 *   z     [  O]  コピー先配列
 *   size  [I  ]  配列サイズ
 */
void Calc::copy(int *a, int *z, int size) {
    int i;  // LOOPインデックス

    for (i = 0; i < size; i++)
        z[i] = a[i];
}

/*
 * 正規化
 * ( Fixed Normalize(A) by Base-Value )
 *
 *   a     [IO]  正規化前データ・正規化後データ配列
 *   size  [I ]  配列数
 */
void Calc::norm(int *a, int size) {
    int i;       // LOOPインデックス
    int cr = 0;  // 繰り上がり

    for (i = size - 1; i >= 1; i--) {
        cr        = (int)(a[i] / 10);
        a[i]     -= cr * 10;
        a[i - 1] += cr;
    }

    for (i = size - 1; i >= 1; i--) {
        if (a[i] < 0) {
            a[i-1] -= 1;
            a[i]   += 10;
        }
    }
}

/*
 * 正規化
 * ( Floating Normalize(A) )
 *
 *   a     [IO]  正規化前データ・正規化後データ配列
 *   size  [I ]  配列数
 */
void Calc::fnorm(int *a, int size) {
    int exp;     // Exponent
    int i;       // LOOPインデックス
    int k;       // 上位の 0 の個数

    // 小数点以下第1位が桁あふれする場合、右に1桁シフト&べき指数加算
    if (a[1] >= 10) {
        for (i = size - 2; i >= 2; i--) a[i + 1] = a[i];
        a[2] = a[1] % 10;
        a[1] = a[1] / 10;
        a[0]++;
    }

    // 正規化前のべき指数を退避
    exp = a[0];

    // 上位の 0 の個数をカウント
    for (i = 1;  i <= size; i++) if (a[i] != 0) break;
    k = i - 1;

    // 上位の 0 の部分に以降の数字をシフト
    for (i = 1; i <= size - k; i++) a[i] = a[i + k];
    // シフト後の下位に 0 をセット
    setv(0, &a[size - k], k);

    // べき指数セット
    a[0] = exp - k;
}

/*
 * 結果出力 (指数表示)
 *
 *   s     [I ]  結果出力データ配列
 *   size  [I ]  配列数
 */
void Calc::fdisp(int *s, int size)
{
    int i;

    printf("0.");

    // 1行に50桁出力
    for (i = 1; i < size; i++) {
        printf("%d", s[i]);
        if (i % 10 == 0 && i != 0) printf(" ");
        if (i % 50 == 0 && i != 0) printf("\n  ");
    }

    // べき指数
    if (s[0] < 0) {
        printf(" * 10^(%d)\n", s[0]);
    } else {
        printf(" * 10^%d\n", s[0]);
    }
    printf("\n");
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 計算
        objCalc.calc();
    }
    catch (...) {
        cout << "例外発生!" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

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

以下のコマンドで C++ ソースをコンパイルする。 (-Wall は警告出力、-O2 最適化のオプション)

1
$ g++ -Wall -O2 -o AddBigFloat AddBigFloat.cpp

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

4. 実行

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

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
$ ./AddBigFloat
A =
0.3166506397 5517322953 4455927857 1987757519 4847779964
  2009701407 3384179042 9825594691 0333535697 5305734848
  1583155516 8414767456 1213087473 5083511072 8760356093
  9260861528 8708350489 6665600522 8344642795 7927074012
  9017577542 9354938200 7911819312 7248098972 1679272395
  5846930083 3567689478 8166093444 9946739747 1949607766
  5572193647 2586159235 1666758625 2802231597 4352987012
  6402277132 8820043314 1578966667 1421949874 0522528069
  4963951540 5549635723 4406781589 7585091247 2928804406
  9421099378 2452951296 9173264624 4306406443 4899106236
  0298616278 8138870616 9750660503 1182164783 7850657611
  3288276226 8598904374 9424601462 5945982526 4035036685
  0724386443 8944764937 8791306165 6940589243 7551129633
  3302554161 8523573699 2849263128 4727992621 9181658920
  8683721043 8022731294 5873938022 0200773879 3895099410
  0697102906 4284562861 9298929860 1880590791 5553901100
  3924018926 2340510501 8564974295 2466060982 5070241212
  5073681212 9586284294 6545779810 2922400553 9499230979
  5540019122 3638888955 4725668677 6322326450 8962752767
  3045813691 6748091836 8374541947 7701484764 0214133418
   * 10^662

B =
0.2140616910 7899137835 2585011554 6960141832 0220145820
  6664385432 1137151554 9572708510 5783131887 2185802977
  5654825676 8538862262 9837519826 7932437106 7012920146
  5568231433 2127581293 9832244270 0278222365 6881641963
  2717062601 9366188326 1079355544 8811537920 3339947001
  8325470026 4247421315 3488593315 6080528840 4031893126
  6836685937 6094269043 1146445759 5370149654 5264893797
  3209843530 5952968044 7265844146 9701696161 3108626440
  8126933347 1362723359 4837891457 6411204558 4740114469
  5270179217 0031153738 5751882220 4053061214 9342679942
  7025123668 8083400521 9366723918 3005539196 1069301376
  6992596887 8188952883 4933095243 0121590290 1841637037
  6071370815 1185334574 6940025938 7104233325 0588311999
  3507700205 2763308632 3384137526 6616586728 4770201096
  3223856628 2844624599 3661681982 8372657212 1407099641
  5771051272 6753929860 0743762636 9858266503 7605017995
  9236522707 3060006656 1533034926 4389931422 1749761459
  0643268852 4164379631 0784863028 9265881878 1317046910
  1179857359 1649348627 6328224925 9398473401 3399233218
  1814868270 8859628160 2132177628 4388966893 0918890608
   * 10^891

C =
0.6311740612 4765751350 6619690871 8569654717 9747225470
  6688549847 9247214805 5121173611 6977943928 8302245627
  3718925488 7676371765 7617032626 5546358018 8976525811
  4897692074 6322977081 0607450943 2021703920 6838355255
  3543581705 0272644156 1155100727 4648199269 9738112620
  2439159344 2826936884 9712656018 2546551605 0254148726
  1155630263 2887454725 2498902997 7304068006 5343011485
  1393305715 4005503524 1777298359 8920272746 4964614968
  5954705731 6400627812 4633294903 7825297890 1561825738
  9964295908 2964486185 4922227902 8914665157 9777338176
  8383525425 6270836663 3820556393 2982254008 5604819478
  7269537178 6174329929 8980663530 3231405114 1758302294
  4964222575 8383423956 3658688920 5893235515 2982446943
  5910880892 0612055726 4774919623 1723923934 2305722325
  1649132015 5329624137 5525200552 2608511035 5089063497
  3447179144 5836609905 4554752011 7652691507 2855635726
  2738305639 9257436798 4161844327 1467219092 9956610291
  2549723175 0928660569 8705666607 9333425932 7639415674
  8718477059 0445987600 1784836515 9943082700 8043563223
  5325020882 3937437727 7303651904 5799292372 7230186850
   * 10^586

A + B =
0.2140616910 7899137835 2585011554 6960141832 0220145820
  6664385432 1137151554 9572708510 5783131887 2185802977
  5654825676 8538862262 9837519826 7932437106 7012920146
  5568231433 2127581293 9832244270 0278222365 6881641963
  2717062601 9366188326 1079355548 0476601895 8513176536
  2884748598 4124996510 1966392957 6177542913 7873683556
  5092632847 9429626018 4203794241 1201704822 9412568358
  5340718265 6788078773 4869405086 2310311450 0192131337
  4782938575 4809151317 4108631586 6586979987 8289496477
  4388372344 2512143460 2544606176 2522362048 5019574730
  8686058118 7550797993 8862801583 8727475668 6930893728
  3660183140 6211268857 8462965369 4144361619 0042070179
  1861037486 5404833315 2165306633 6743748730 6084669233
  7575516102 8614221105 2672181596 0827580510 9299714065
  4956502872 5908689034 2652744343 1358820000 2795805811
  3277656303 8575577697 9250338750 2741028772 3594061745
  3482537333 2519831920 1883401777 1633795861 1197410837
  8556330509 3570272068 6296159362 2291423496 6552783902
  9672488643 8929274846 8144814134 6235683839 3626546164
  0554248491 0867366954 1083171729 1359995957 3764519227
   * 10^891

A - C =
0.3166506397 5517322953 4455927857 1987757519 4847779964
  2009701407 3384179042 9825588379 2927410931 7792228228
  4674436947 1867587708 8958380784 9585031825 6612300972
  7524744550 9269062187 4209326803 9089755119 4209416395
  8691311996 5774749224 2653704415 0327352649 1908461788
  1337498061 6528482640 4613539901 4129689474 5508046611
  4564918999 0593459497 0540556186 1208788770 4984137299
  9842094586 3303993060 0091705511 5119316986 5975275570
  5933974236 4869570380 4291930196 4527937241 7893562629
  6437500457 9725486332 3023578669 7249090042 8620981602
  7349578453 5159969055 1493270538 8223081819 2988802689
  1009247312 1947324597 6042833079 0691726255 5668402864
  5160453461 6404679333 0596518896 1568803069 4251830652
  6667250930 4472431940 9826318164 2502234238 4942095262
  1794515150 5667578312 1404502111 1391853267 3337834635
  1500871182 5045220556 2075678211 0560435461 9312525575
  1918496317 7230155412 7930000848 0674615145 8971186657
  7553563560 2671211439 0188517071 9866001296 5131246817
  7096747655 1447958998 8622756127 9090575522 2357054061
  6379734358 2488764197 4217793229 2930894318 0338131634
   * 10^662

桁をずらして加減算してみると、正しく計算できていることが分かる。


考え方としては難しいものではありませんが、一度実際に実装してみることでより理解が深まります。

以上。

Comments