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 という浮動小数点を以下のように考える。（１桁を１個の配列で扱う場合）

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

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

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

以下は、加算の場合。（減算も同様）

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

2. C++ ソース作成

• 被加減数・加減数の桁数は合わせている。
• テストなので、被加減数・加減数は適当な乱数で生成している。
• 今後汎用性を高めるための処置として、色々と部品化している。
• ソースファイルは１つ。（もう少し大きくなるか複雑になれば、ファイルを分割するかもしれないが）

File: add_big_float.cpp

/*********************************************
 * 浮動小数点数加減算
 *********************************************/
#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 の個数

    // 小数点以下第１位が桁あふれする場合、右に１桁シフト＆べき指数加算
    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;
}

• Gist - C++ source code to compute big-digit floats.

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

（-Wall は警告出力、-O2 最適化のオプション）

$ g++ -Wall -O2 -o add_big_float add_big_float.cpp

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

4. 実行

$ ./add_big_float
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

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

---

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

以上。