C++ - 多桁乗算（標準（筆算）法）！

以前、多桁（データ型を超える整数）の加減乗除アルゴリズムを C++ に実装してみました。

• C++ - 多桁計算！
• C++ - 多桁計算（その２）！

ただ、乗算は一方が多桁でもう一方が小さい桁数と限定していました。

今回は、「標準（筆算）方式」による多桁同士の乗算アルゴリズムを C++ で実装してみました。

0. 前提条件

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

また、当方の環境で扱える int 型の範囲は以下のとおり。

• int : -2,147,483,648 〜 2,147,483,647

1. 多桁計算について

「筆算方式」の考え方は、小学校で習ったとおりであるが、今回は配列は下位から取るようにしている。

2. C++ ソース作成

例として、以下のようにソースを作成した。概要は以下のとおり。

• 上記の説明では1個の配列で4桁を扱うようにしているが、以下のソースでは1個の配列で1桁を扱うようにしている。
• 計算可能な桁数は2のべき乗桁としている。
  （後に筆算方式以外のアルゴリズムを利用することを考えて）
• 繰り上がり処理は、最後にまとめて行う（但し、非常に大きい乗算桁数（数千万桁以上）では桁あふれを起こすので注意）
• 計算に使用する被乗数・乗数は、手入力は困難なため、乱数を使用している。
• 冒頭の // #define TEST は、乗算回数をカウントしたり、処理時間を計測するテストを行うため。
  テストを行うならコメントを解除する。
  （テストの処理時間計測では clock 関数を使用しているため、あまり精度はよくない）

File: multiply_normal.cpp

/*********************************************
 * 多倍長乗算 ( by 標準(筆算)方式 )
 *   - 多倍長 * 多倍長
 *   - 最下位の桁を配列の先頭とする考え方
 *********************************************/
#include <cstdlib>   // for rand()
#include <iostream>  // for cout
#include <math.h>    // for pow()
#include <stdio.h>   // for printf()

//#define TEST       // テスト ( 乗算回数・処理時間 ) するならコメント解除
#define D_MAX 1024   // 計算可能な最大桁数 ( 2のべき乗 )
#define D     1024   // 実際に計算する桁数 ( D_MAX 以下の 4 の倍数 )

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    int A[D];  // 被乗数配列
    int B[D];  // 乗数配列
#ifdef TEST
    int cnt_mul;       // 乗算回数
    clock_t t1, t2;    // 計算開始CPU時刻、計算終了CPU時刻
    double tt;         // 計算時間
#endif

    public:
        Calc();                                         // コンストラクタ
        void calcNormal();                              // 計算 ( 標準(筆算)法 )
    private:
        void multiplyNormal(int *, int *, int, int *);  // 乗算 ( 標準(筆算)法 )
        void doCarry(int *, int);                       // 繰り上がり・借り処理
        void display(int *, int *, int *);              // 結果出力
};

/*
 * コンストラクタ
 */
Calc::Calc()
{
    /* ====================================== *
     * テストなので、被乗数・乗数は乱数を使用 *
     * ====================================== */

    int i;  // LOOP インデックス

    // 被乗数・乗数桁数設定
    for (i = 0; i < D; i++) {
        A[i] = rand() % 10;
        B[i] = rand() % 10;
    }
}

/*
 * 計算 ( 標準(筆算)法 )
 */
void Calc::calcNormal()
{
    int a[D_MAX];      // 被乗数配列
    int b[D_MAX];      // 乗数配列
    int z[D_MAX * 2];  // 計算結果用配列
    int i;             // LOOPインデックス

#ifdef TEST
    t1 = clock();  // 計算開始時刻
    for (int l = 0; l < 1000; l++) {
        cnt_mul = 0;  // 乗算回数リセット
#endif
    // 配列初期設定
    //for (i = 0; i < D / LIMIT; i++) {
    for (i = 0; i < D; i++) {
        a[i] = A[i];
        b[i] = B[i];
    }

    // 最大桁に満たない部分は 0 を設定
    //for (i = D / LIMIT; i < ASIZE; i++) {
    for (i = D; i < D_MAX; i++) {
        a[i] = 0;
        b[i] = 0;
    }

    // 乗算 ( 標準(筆算)法 )
    multiplyNormal(a, b, D_MAX, z);
#ifdef TEST
    }
    t2 = clock();  // 計算終了時刻
    tt = (double)(t2 - t1) / CLOCKS_PER_SEC;  // 計算時間
#endif

    // 繰り上がり・借り処理
    doCarry(z, D_MAX * 2);

    // 結果出力
    display(a, b, z);
}

/*
 * 乗算 ( 標準(筆算)法 )
 */
void Calc::multiplyNormal(int *a, int *b, int tLen, int *z)
{
    int i, j;  // ループインデックス

    // 計算結果初期化
    for(i = 0; i < tLen * 2; i++) z[i] = 0;

    // 各配列を各桁とみなして乗算
    for (j = 0; j < tLen; j++) {
        for (i = 0; i < tLen; i++) {
            z[j + i] += a[i] * b[j];
#ifdef TEST
            cnt_mul++;  // 乗算カウント
#endif
        }
    }
}

/*
 * 繰り上がり・借り処理
 */
void Calc::doCarry(int *a, int tLen) {
    int cr;  // 繰り上がり
    int i;   // ループインデックス

    cr = 0;
    for(i = 0; i < tLen; i++) {
        a[i] += cr;
        if(a[i] < 0) {
            //cr = -(-(a[i] + 1) / (int)pow(10, LIMIT) + 1);
            cr = -(-(a[i] + 1) / 10 + 1);
        } else {
            cr = a[i] / 10;
        }
        a[i] -= cr * 10;
    }

    // オーバーフロー時
    if (cr != 0) printf("[ OVERFLOW!! ] %d\n", cr);
}

/*
 * 結果出力
 */
void Calc::display(int *a, int *b, int *z)
{
    int i;  // LOOPインデックス

    // 上位桁の不要な 0 を削除するために、配列サイズを取得
    int aLen = D_MAX, bLen = D_MAX, zLen = D_MAX * 2;
    while (a[aLen - 1] == 0) if (a[aLen - 1] == 0) aLen--;
    while (b[bLen - 1] == 0) if (b[bLen - 1] == 0) bLen--;
    while (z[zLen - 1] == 0) if (z[zLen - 1] == 0) zLen--;

    // a 値
    printf("a =\n");
    for (i = aLen - 1; i >= 0; i--) {
        printf("%d", a[i]);
        if ((aLen - i) % 10 == 0) printf(" ");
        if ((aLen - i) % 50 == 0) printf("\n");
    }
    printf("\n");

    // b 値
    printf("b =\n");
    for (i = bLen - 1; i >= 0; i--) {
        printf("%d", b[i]);
        if ((bLen - i) % 10 == 0) printf(" ");
        if ((bLen - i) % 50 == 0) printf("\n");
    }
    printf("\n");

    // z 値
    printf("z =\n");
    for (i = zLen - 1; i >= 0; i--) {
        printf("%d", z[i]);
        if ((zLen - i) % 10 == 0) printf(" ");
        if ((zLen - i) % 50 == 0) printf("\n");
    }
    printf("\n\n");

#ifdef TEST
    printf("Counts of multiply / 1 loop = %d\n", cnt_mul);     // 乗算回数
    printf("Total time of all loops     = %f seconds\n", tt);  // 処理時間
#endif
}

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

        // 乗算 ( 標準(筆算)法 )
        objCalc.calcNormal();
    }
    catch (...) {
        cout << "例外発生！" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

• Gist - C++ source code to multiply big-digit values.

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

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

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

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

4. 実行

$ ./multiply_normal
a =
5541050205 0676367860 0766496541 6816500114 2004376325 
9767889290 6255254084 0240121540 8107405889 9787092091 
5419878423 0874067594 2003014938 3144522637 4398888312 
7310798560 7723593059 2397340736 1413730264 5938118523 
4930451302 5419027365 4464726284 4804570926 7856292055 
6998051538 5419752372 6773282553 8808373348 1147925911 
4359623467 7539209376 4653165686 1586483987 2549340871 
3021915969 1733401776 8399208132 2374597797 5435289449 
3231696299 4733994450 9183503889 6019289550 6163237384 
1227614331 9992503157 1308593784 1932936594 0929424801 
1950071971 0154922124 5193227045 8846648127 7190389019 
7738758027 8658264055 2312943257 6472467201 4300012776 
8190966928 1447873191 5142312908 4227930383 4209210656 
3795585493 1003705337 6408294366 3725395154 3867977667 
2088616184 7174172213 2970374580 1779335759 6608705902 
2427117686 6025440760 1301612932 7987359459 5970899460 
8046752633 3980870333 6828242461 4139003942 9948802522 
0620475091 9136269879 5660460399 9031389202 4282408151 
0692857294 5420047984 6336838598 4572973744 3149168608 
0328966779 1710456127 7745469624 3867344243 6645626135 
8419139297 2271203029 6373
b =
1396721814 7429756767 9697316781 5799720913 8264945924 
4078224610 6827742845 2280304370 0503322098 7712258675 
0868106360 9153608487 2203105281 1296790491 3736052352 
3603808771 3400665829 7092057744 8109671179 4731531344 
3030270662 8830811711 2659042336 9560268360 3635996665 
4439774037 4364813709 6286856908 2604898197 3917156409 
3643090100 1850620740 4921392254 2053171203 9180636269 
2516952585 7247460426 4661757895 1691002203 5994933942 
3774155041 9930026893 1701001730 3020409083 4387572767 
2892156396 8494440010 8636358355 3901596440 6681549762 
3712160872 3463720621 7944173552 7388558576 9726708891 
6702305518 3009506306 1062776798 5121255133 1714952137 
9438940915 6880693386 7752483286 3454240702 0741628501 
7609630135 1485681831 5729573433 4022884541 9500149345 
9210492650 0111751027 4399813889 9483783897 0746831499 
0567300092 4668876475 9172290045 8910594972 4181323728 
7410986527 3076964838 8889816294 3754256945 6100329502 
4194808714 9521709228 6724065836 1340110964 7865665010 
2838147426 9059336788 0834918375 1781651644 0088515395 
8991435367 5179991062 2180962809 1488047557 9575103060 
4793267825 3098666971 2556
z =
7739305698 0040071707 2189003016 6876576420 6222934577 
2023571781 7526420167 7391080311 4489897836 5606748456 
8717128396 3330797577 7919515520 8809157391 6002244534 
7621606312 5432389533 8866984447 8836029874 7409726049 
4007356506 3340487227 9506520886 5423598238 9290074291 
2974566830 6442700798 4333021269 4431263626 6374539628 
6851995363 4864882970 5572176762 3085210352 4275581437 
0897142486 7452815352 0137113684 7907862599 4876133849 
0555712981 8780385401 5355569540 8154357653 2471295504 
9948804529 3912328019 6503502489 7665018794 4772954562 
9490632556 4789219622 8100182409 9768422449 2193203215 
2014712520 6011177626 0750743785 4179584096 4152952580 
0306309481 6426383088 4574821896 2120453447 6352845378 
1951677793 9159184937 5572755307 5303031964 7543670100 
5278150467 3261231786 8254099723 7458674159 9839559097 
0426560067 4606133160 9122200444 2711896648 1817977034 
5524752649 8581262782 1602383986 9796468261 1996639974 
2576401764 9436438874 1313980272 5810664985 4229675065 
8756237666 5532150529 3593433033 9405077272 1609682116 
1412446780 1520707085 6922581219 2890317519 6136848317 
5503977733 8847761461 5785599640 3027734783 2267048714 
4548525727 8351825512 2431118909 0134095266 6304941694 
5445017714 0685438098 2538356667 9630620121 3564506833 
4791511337 3608870632 1523183163 7757254037 7900758746 
7454636070 9231962770 4771204380 7741285769 9964102181 
0886535087 9117351045 4341821086 0872546443 4774459105 
1523897188 5238232263 9284079260 9689690063 7880596067 
0489076973 1329053634 9696312126 8798891537 2027232304 
0211601447 5986863179 6199396248 7964245496 2337368324 
3233702362 0325033738 4242202489 2049900204 5559372819 
3157276326 8950937657 0461483163 3231712199 0249152697 
0627420697 7960887162 9208245972 1570702376 5711273790 
4233182037 3921249831 5354876969 9937253686 0092967148 
0764175165 1517752657 8361489911 0434806183 4354010504 
0070277791 9744702000 8959936173 5978939803 1722903884 
7001691849 0262982850 8638795795 0315838055 3558236720 
1363881352 6338316728 0294918962 7054187105 8431713068 
4930188748 7096773537 8631442377 9673526455 0412143942 
3297377182 0841986843 4612914627 7539475285 2246449343 
4643618733 2015199866 9194241922 3316370313 5694752406 
0784066772 8685693975 1607530824 8097551139 9359388

部分的に電卓で計算してみた結果、間違いは無さそう。

5. 検証

上記のコードを利用して、乗算回数が何回になるか、計算にどれくらい時間がかかるかを検証してみた。
乗算回数は予想通り配列数の２乗回に、計算にかかる時間も乗算回数に比例した時間となった。

---

今回は1桁を1個の配列で扱ったが、乗算する桁数があらかじめ分かっているのなら、1個の配列で4桁を扱ったりすると速度が向上する。（但し、桁あふれに注意）
また、今回は繰り上がり処理を最後にまとめて行なっているが、乗算の都度行うとさらに速度が向上する。（但し、繰り上がりにより配列の個数が変動するので厄介）

単純に多桁同士の乗算を計算するだけなら「筆算方式」で充分ですが、もっと桁数を増やしたり計算速度を上げる必要がある場合は、別のアルゴリズムを使用するようです。

以上。