C++で実装するSchönhage-Strassenアルゴリズム ― 2つの数値を高速に乗算するプログラム

Schönhage-Strassenアルゴリズムとは

Schönhage-Strassen（ショーンハーゲ・シュトラッセン）アルゴリズムは、2つの数値を乗算するために用いられるアルゴリズムです。非常に大きな整数の乗算において漸近的に高速に動作することで知られており、実際には 2²¹⁵〜2²¹⁷（10進数で約10,000〜40,000桁）を超えるような巨大な数値を扱う場合に、Karatsuba法やToom-Cook法といった従来の乗算手法を上回る性能を発揮し始めます。

アルゴリズムの流れ

この実装では、次の3つのステップで乗算を行います。

1. 桁数のカウント

関数 noOfDigit() を使って、入力された2つの整数がそれぞれ何桁かを求めます。

2. 線形畳み込みの計算

2つの数値の各桁同士の積をすべて計算し、対応する位置にある配列 linearConvolution に加算していきます。これにより、筆算の途中過程に相当する「線形畳み込み」が得られます。

3. 繰り上がりの処理

線形畳み込みの各要素を下位の桁から順に走査し、10以上の値を繰り上がりとして上位の桁へ渡しながら、最終的な積を組み立てます。

以上の手順を擬似コードで表すと、次のようになります。

開始
    関数 noOfDigit(x)
        変数 n を宣言し、n = 0 とする
        while (x > 0)
            x = x / 10
            n を 1 増やす
        return n
終了

開始
    schonhageStrassenMultiplication のアルゴリズム:
    schonhageStrassenMultiplication(a, b, n, m)
        配列 linearConvolution[n + m - 1] を定義する
        for i = 0 to (n + m - 2)
            linearConvolution[i] = 0
        p = a
        for i = 0 to m-1
            a = p
            for j = 0 to n-1
                linearConvolution[i + j] += (b mod 10) * (a mod 10)
                a /= 10
            b /= 10
        for i = (n + m - 2) down to 0
            linearConvolution[i] を出力する
        product = 0、nextCarry = 0、base = 1 とする
        for i = 0 to (n + m - 2)
            linearConvolution[i] += nextCarry
            product = product + (base * (linearConvolution[i] % 10))
            nextCarry = linearConvolution[i] / 10
            base *= 10
        2つの数値の積を出力する
終了

C++による実装例

続いて、このアルゴリズムをC++で実装したサンプルコードを示します。標準入力から2つの整数を読み込み、その積を計算して出力します。

#include <iostream>
using namespace std;

// 桁数を数える関数
int noOfDigit(long x) {
    int n = 0;
    while (x > 0) {
        x /= 10;
        n++;
    }
    return n;
}

// Schonhage-Strassen法による乗算
void schonhageStrassenMultiplication(long a, long b, int n, int m) {
    int linearConvolution[n + m - 1];
    for (int i = 0; i < (n + m - 1); i++)
        linearConvolution[i] = 0;
    long p = a;
    for (int i = 0; i < m; i++) {
        a = p;
        for (int j = 0; j < n; j++) {
            linearConvolution[i + j] += (b % 10) * (a % 10);
            a /= 10;
        }
        b /= 10;
    }
    cout << "The Linear Convolution is: ( ";
    for (int i = (n + m - 2); i >= 0; i--) {
        cout << linearConvolution[i] << " ";
    }
    cout << ")";
    long product = 0;
    int nextCarry = 0, base = 1;
    for (int i = 0; i < n + m - 1; i++) {
        linearConvolution[i] += nextCarry;
        product = product + (base * (linearConvolution[i] % 10));
        nextCarry = linearConvolution[i] / 10;
        base *= 10;
    }
    cout << "\nThe Product of the numbers is: " << product;
}

int main(int argc, char **argv) {
    cout << "Enter the numbers:";
    long a, b;
    cin >> a >> b;
    int n = noOfDigit(a);
    int m = noOfDigit(b);
    schonhageStrassenMultiplication(a, b, n, m);
}

実行結果

Enter the numbers:1234 5679
The Linear Convolution is: ( 5 16 34 61 63 55 36 )
The Product of the numbers is: 7007886

この実行例では、1234 × 5679 の計算を行っています。まず各桁の積から構成される線形畳み込み「5 16 34 61 63 55 36」が出力され、その後、下位の桁から順に繰り上がりを処理することで、最終的な積「7007886」が正しく求められていることがわかります。