第4章: 高速乗算アルゴリズム

1. Karatsuba法

1962年、Anatolii Karatsuba は、$n$ 桁の乗算に必ず $O(n^2)$ かかるという Kolmogorov の予想を覆す画期的なアルゴリズムを発見した。 そのアイデアは単純かつ美しい：$n$ 桁の乗算を 4 回の $n/2$ 桁乗算ではなく 3 回の $n/2$ 桁乗算に帰着させるのである。

1.1 原理

1.2 定理と証明

1.3 計算量の解析

1.4 擬似コード

function karatsuba(A[0..n-1], C[0..n-1]):
    // 基底ケース：小さな整数は筆算で計算
    if n <= THRESHOLD:
        return schoolbook_multiply(A, C)

    // 分割
    m = ceil(n / 2)
    A0 = A[0..m-1]         // 下位 m リム
    A1 = A[m..n-1]         // 上位リム
    C0 = C[0..m-1]         // 下位 m リム
    C1 = C[m..n-1]         // 上位リム

    // 3 回の再帰的乗算
    z0 = karatsuba(A0, C0)              // A0 * C0
    z2 = karatsuba(A1, C1)              // A1 * C1
    z1 = karatsuba(A0 + A1, C0 + C1)   // (A0+A1)(C0+C1)
    z1 = z1 - z0 - z2                  // 交差項

    // 結合（シフトと加算）
    result = z0
    result += z1 << (m * BITS_PER_LIMB)    // z1 * B^m
    result += z2 << (2*m * BITS_PER_LIMB)  // z2 * B^(2m)
    return result

1.5 具体的な数値例

Karatsuba法の可視化

2. Toom-Cook法

Toom-Cook法は Karatsuba法の自然な一般化であり、 整数を多項式として表現し、多項式補間の技法を用いて乗算を行う。 Karatsuba法が 2分割・3点評価であるのに対し、 Toom-$k$ 法は $k$ 分割・$(2k-1)$ 点評価を用いる。

2.1 Toom-3法の原理

2.2 正当性の証明

2.3 補間の具体的な公式

2.4 計算量の解析

2.5 一般のToom-$k$ 法

3. FFT乗算

Toom-Cook法では分割数 $k$ を増やすほど漸近的な指数が $1$ に近づくが、 補間のオーバーヘッドが増大する。FFT（高速フーリエ変換）を用いると、 本質的に $O(n \log n)$ に近い計算量で乗算を実現できる。

3.1 乗算と畳み込み

3.2 畳み込み定理

3.3 FFT乗算のアルゴリズム

3.4 浮動小数点FFTの問題と対策

3.5 NTT（数論変換）

3.6 Schönhage-Strassen アルゴリズム

3.7 Harvey-van der Hoeven の結果（2019年）

4. 方法の比較とアルゴリズムの切り替え

4.1 計算量の比較

4.2 アルゴリズムの自動切り替え

/* mpn/x86_64/gmp-mparam.h の例 */
#define MUL_TOOM22_THRESHOLD     20   /* 筆算 → Karatsuba */
#define MUL_TOOM33_THRESHOLD     65   /* Karatsuba → Toom-3 */
#define MUL_TOOM44_THRESHOLD    166   /* Toom-3 → Toom-4 */
#define MUL_TOOM6H_THRESHOLD    222   /* → Toom-6 */
#define MUL_TOOM8H_THRESHOLD    309   /* → Toom-8 */
#define MUL_FFT_THRESHOLD       4736  /* → FFT (Schönhage-Strassen) */

4.3 非対称乗算

FFT乗算の可視化

5. GPU 並列計算と多倍長乗算

GPU（Graphics Processing Unit）は、数千〜数万のコアによる大規模並列処理を得意とするプロセッサである。 深層学習や科学技術計算で革命的な性能を発揮した GPU は、多倍長計算にも適用できるのか？ 本節では、GPU のアーキテクチャ特性を踏まえた上で、多倍長乗算への応用の可能性と課題を検討する。

5.1 GPU のアーキテクチャと SIMT モデル

NVIDIA の GPU は SIMT（Single Instruction, Multiple Threads）モデルに基づく。 32 スレッドを1単位とする ワープ（warp） が同一命令を同時に実行し、 SM（Streaming Multiprocessor）上の数千スレッドが協調動作する。 この構造は、同じ演算を大量のデータに適用する データ並列性 が高い問題に最も適している。

多倍長計算を GPU に移植する場合、主に二つのアプローチが考えられる：

1. 1 数 = 多スレッド方式：1つの巨大な多倍長数の演算を多数のスレッドで並列化する。FFT/NTT 乗算が典型例。
2. 1 数 = 1 スレッド方式：中程度の精度（数百〜数千ビット）の多倍長数を各スレッドが独立に処理する。モンテカルロ法や暗号演算が典型例。

5.2 GPU 上の FFT/NTT 乗算

FFT は本質的にデータ並列性が高く、GPU との相性が良い。 cuFFT ライブラリは倍精度浮動小数点の FFT を高速に実行でき、これを多倍長乗算の基盤として利用できる。

しかし、多倍長整数の乗算に FFT を用いる場合、NTT（Number Theoretic Transform）のほうが適切である。 NTT は有限体 $\mathbb{F}_p$ 上の FFT であり、浮動小数点の丸め誤差が一切発生しない。 GPU 上での NTT 実装の基本構造は以下の通りである：

// GPU 上の NTT の概略（CUDA 擬似コード）
__global__ void ntt_kernel(uint64_t* a, int n, uint64_t w, uint64_t p) {
    // バタフライ演算をスレッド並列で実行
    for (int len = 1; len < n; len <<= 1) {
        int tid = blockIdx.x * blockDim.x + threadIdx.x;
        if (tid < n / 2) {
            int i = (tid / len) * (2 * len) + (tid % len);
            int j = i + len;
            uint64_t u = a[i];
            uint64_t v = mulmod(a[j], w_table[...], p);
            a[i] = addmod(u, v, p);
            a[j] = submod(u, v, p);
        }
        __syncthreads();
    }
}

GPU 上の NTT 乗算の性能は、メモリアクセスパターンに大きく依存する。 バタフライ演算の各段階で、アクセスするデータのストライドが 1, 2, 4, ... と変化するため、 初期段階ではメモリアクセスが非連続になり バンクコンフリクト が発生しやすい。 共有メモリの活用と適切なデータ配置が性能の鍵を握る。

5.3 多倍長計算における GPU の課題

GPU による多倍長計算には、アーキテクチャに起因するいくつかの本質的な課題がある。

（a）桁上がり伝播の逐次性

多倍長加算の桁上がり（carry）は本質的に逐次的な処理である。 最悪の場合、$n$ 桁の加算で桁上がりが全桁にわたって伝播する（例：$999\ldots9 + 1 = 1000\ldots0$）。 これは CPU 上では単純なループで処理されるが、GPU の並列モデルとは相容れない。

この問題に対する主なアプローチは：

• 冗長表現（redundant representation）：各リムに余裕を持たせ、桁上がりの伝播を遅延させる。最終的な正規化のみ逐次処理。
• 並列プレフィックス加算（parallel prefix sum）：Brent-Kung または Kogge-Stone アルゴリズムにより、$O(\log n)$ 段で全桁の桁上がりを確定。
• 繰り上がり先見加算器の原理：Generate/Propagate ビットを並列計算し、桁上がりを $O(\log n)$ で解決。

（b）ワープダイバージェンス

SIMT モデルでは、ワープ内の 32 スレッドが異なる分岐を取ると性能が低下する（ワープダイバージェンス）。 多倍長計算では、数値の桁数が不揃いの場合や条件分岐が発生する場合に、このペナルティが生じる。 たとえば、GCD 計算では繰り返しの終了条件がデータ依存であるため、GPU 並列化が困難である。

（c）メモリ帯域幅のボトルネック

多倍長乗算は 演算強度（Arithmetic Intensity）が比較的低い。 $n$ リムの乗算で $O(n \log n)$ 回の演算と $O(n)$ 回のメモリアクセスが発生するため、 演算強度は $O(\log n)$ 程度である。 GPU の演算性能がメモリ帯域幅に対して圧倒的に高い（A100 で約 5:1 のバイト/FLOP 比率）ため、 小〜中規模の多倍長数では メモリ律速 になりやすい。

5.4 既存のライブラリと実装

GPU 上の多倍長計算に関して、いくつかの研究プロジェクトとライブラリが存在する。

5.5 性能特性と適用領域

GPU が CPU に対して優位性を持つ場面と持たない場面を整理する。

5.6 GPU 多倍長計算の将来展望

現時点では、GPU が CPU ベースの多倍長ライブラリ（GMP, MPFR）を全面的に置き換える段階には至っていない。 しかし、以下の技術進化により、GPU の適用領域は着実に拡大している。

• 統合メモリ（Unified Memory）：CPU-GPU 間のデータ転送オーバーヘッドが削減され、中規模の多倍長計算でも GPU を自然に活用できる。
• ワープレベルプリミティブ：__shfl_sync 等のワープ内データ交換命令により、レジスタ間でリムの桁上がり情報を高速に伝播できるようになった。
• Tensor Core の応用：行列乗算に特化した Tensor Core を、整数乗算の部分積計算に転用する研究がある。INT8/INT4 の Tensor Core が利用可能。
• チップレットとマルチ GPU：NVLink/NVSwitch による GPU 間の高速通信により、超大規模な FFT を複数 GPU に分散して実行可能。

結論として、GPU 多倍長計算は 「大量の独立した演算」または 「超大規模な単一乗算」というニッチにおいて CPU を大幅に凌駕する。 一方、汎用的な多倍長ライブラリの全機能を GPU 上で効率的に実現することは、 桁上がり伝播の逐次性やワープダイバージェンスといった本質的な課題により、依然として困難である。 最も現実的なアプローチは、CPU と GPU のハイブリッド方式 — CPU が逐次的な制御と小規模演算を担当し、 大規模な乗算や大量の独立演算を GPU にオフロードする — であろう。

6. 分散並列計算（MPI/OpenMP）

前節の GPU 並列化は単一ノード内の並列性を活用するものであった。 数十億桁〜数兆桁規模の計算（円周率の世界記録など）では、 複数ノードにまたがる分散並列計算が不可欠となる。 本節では、多倍長計算における MPI（Message Passing Interface）と OpenMP の活用を概説する。

6.1 OpenMP によるスレッド並列化

OpenMP は共有メモリ環境でのマルチスレッド並列化フレームワークである。 多倍長計算において OpenMP が有効な場面を以下に整理する。

void karatsuba_parallel(mpz_t result, const mpz_t a, const mpz_t b, int depth) {
    size_t n = MAX(mpz_size(a), mpz_size(b));

    if (n < KARATSUBA_THRESHOLD) {
        mpz_mul(result, a, b);  /* 筆算にフォールバック */
        return;
    }

    /* a = a1 * B + a0,  b = b1 * B + b0 に分割 */
    size_t half = n / 2;
    mpz_t a0, a1, b0, b1, z0, z1, z2;
    /* ... 分割処理 ... */

    if (depth < MAX_PARALLEL_DEPTH) {
        #pragma omp task shared(z0)
        karatsuba_parallel(z0, a0, b0, depth + 1);  /* z0 = a0 * b0 */

        #pragma omp task shared(z2)
        karatsuba_parallel(z2, a1, b1, depth + 1);  /* z2 = a1 * b1 */

        /* z1 = (a0+a1)(b0+b1) - z0 - z2 */
        mpz_add(a0, a0, a1);
        mpz_add(b0, b0, b1);
        karatsuba_parallel(z1, a0, b0, depth + 1);

        #pragma omp taskwait
        mpz_sub(z1, z1, z0);
        mpz_sub(z1, z1, z2);
    } else {
        /* 逐次実行にフォールバック */
        karatsuba_parallel(z0, a0, b0, depth + 1);
        karatsuba_parallel(z2, a1, b1, depth + 1);
        /* ... */
    }

    /* result = z2 * B^2 + z1 * B + z0 */
    /* ... 結合処理 ... */
}

6.2 MPI による分散メモリ並列化

単一マシンのメモリを超える規模の計算や、多数のノードで加速したい場合は MPI を用いた分散メモリ並列化が必要となる。 多倍長乗算における MPI の主要なパターンは以下の通りである。

• 分散 FFT/NTT：巨大なリム配列を各ノードに分割配置し、 FFT のバタフライ演算をノード間通信（MPI_Alltoall）と組み合わせて実行する。 通信量は各段 $O(N/P)$（$P$ はプロセス数）だが、全対全通信のレイテンシが支配的になりやすい
• CRT ベースの分散乗算：複数の NTT 素数を各ノードに割り当て、 ノード内で独立に NTT 乗算を実行した後、CRT で合成する。 通信は最終段の合成のみで、通信効率が極めて高い
• Binary Splitting の分散化：再帰木の高レベルの分岐を異なるノードに割り当て、 最終段で結果を合成する

6.3 y-cruncher の並列化戦略

円周率の世界記録を何度も更新してきた y-cruncher は、 多倍長計算の並列化における実践的な手本である。以下の多層的な並列化を行っている。

• SIMD：AVX-512 による NTT バタフライのベクトル化
• マルチスレッド：OpenMP/独自スレッドプールによるコアレベル並列化
• ディスクスワップ：メモリに収まらない巨大データをディスクに退避し、 アウトオブコア（out-of-core）アルゴリズムで処理する
• NUMA 最適化：メモリアフィニティを考慮したデータ配置

まとめ

本章では、多倍長整数の乗算を高速化する主要なアルゴリズムを学んだ。

• Karatsuba法（1962年）：2分割・3回乗算で $O(n^{1.585})$。最も基本的な高速乗算であり、$n$ 桁の乗算を 3 回の $n/2$ 桁乗算に帰着させる。筆算の 4 回を 3 回に減らすアイデアは、交差項 $A_0 C_1 + A_1 C_0$ を $(A_0 + A_1)(C_0 + C_1) - A_0 C_0 - A_1 C_1$ として 1 回の乗算で求めることにある。
• Toom-Cook法：$k$ 分割・多項式補間で $O(n^{\log_k(2k-1)})$。Karatsuba の一般化であり、多項式の「評価→点ごとの乗算→補間→再構成」という枠組みで理解できる。$k$ を増やすほど指数は小さくなるが、補間のオーバーヘッドが増大するため、実用上は Toom-3, Toom-4 が限界。
• FFT乗算：畳み込み定理により、乗算を周波数領域での成分ごとの積に帰着させる。$O(n \log n)$ で畳み込みが計算できる。
• NTT（数論変換）：有限体上のFFTで丸め誤差を完全に排除する。整数乗算に不可欠な技術。
• Schönhage-Strassen（1971年）：フェルマー数環上でのFFTにより $O(n \log n \log \log n)$ を達成。GMP の大規模乗算の中核。
• Harvey-van der Hoeven（2019年）：理論的に最適な $O(n \log n)$ を達成したが、定数項が巨大で実用的ではない。
• アルゴリズムの切り替え：入力サイズに応じて最適なアルゴリズムを自動選択することが、実用的な多倍長ライブラリの設計において極めて重要である。
• GPU 並列計算：大量の独立した多倍長演算や超大規模 NTT 乗算では GPU が CPU を大幅に凌駕するが、桁上がり伝播の逐次性やワープダイバージェンスにより汎用的な適用は依然として課題が多い。CPU と GPU のハイブリッド方式が現実的な解である。
• 分散並列計算：数兆桁規模の計算では MPI による分散 FFT/NTT や CRT ベースの分散乗算が不可欠。OpenMP によるスレッド並列化は Karatsuba の再帰分割や Binary Splitting で有効であり、y-cruncher は SIMD・マルチスレッド・ディスクスワップ・NUMA を多層的に活用する。

次章では、これらの乗算アルゴリズムを基盤として、GCD（最大公約数）、モジュラ演算、べき乗剰余、素数判定など、多倍長整数の高度な演算を扱う。