Chapter 4: Fast Multiplication Algorithms

1. Karatsuba Algorithm

In 1962 Anatolii Karatsuba discovered a landmark algorithm that overturned Kolmogorov's conjecture that $n$-digit multiplication necessarily requires $O(n^2)$ time. The idea is simple and beautiful: reduce $n$-digit multiplication to three $n/2$-digit multiplications instead of four.

1.1 Principle

1.2 Theorem and proof

1.3 Complexity analysis

1.4 Pseudocode

function karatsuba(A[0..n-1], C[0..n-1]):
    // Base case: small integers use schoolbook
    if n <= THRESHOLD:
        return schoolbook_multiply(A, C)

    // Split
    m = ceil(n / 2)
    A0 = A[0..m-1]         // lower m limbs
    A1 = A[m..n-1]         // upper limbs
    C0 = C[0..m-1]         // lower m limbs
    C1 = C[m..n-1]         // upper limbs

    // Three recursive multiplications
    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                  // cross term

    // Combine (shift + add)
    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 Concrete numeric example

Visualization of the Karatsuba algorithm

2. Toom-Cook Algorithm

The Toom-Cook algorithm is a natural generalization of Karatsuba: it represents integers as polynomials and performs multiplication via polynomial interpolation. While Karatsuba uses a 2-way split with 3 evaluation points, Toom-$k$ uses a $k$-way split with $2k - 1$ evaluation points.

2.1 Principle of Toom-3

2.2 Proof of correctness

2.3 Concrete interpolation formula

2.4 Complexity analysis

2.5 General Toom-$k$

3. FFT Multiplication

In Toom-Cook, increasing $k$ pushes the asymptotic exponent toward $1$, but the interpolation overhead grows. Using the Fast Fourier Transform (FFT), one can attain essentially $O(n \log n)$ time for multiplication.

3.1 Multiplication and convolution

3.2 The convolution theorem

3.3 The FFT multiplication algorithm

3.4 Floating-point FFT issues and remedies

3.5 NTT (Number-Theoretic Transform)

3.6 The Schönhage-Strassen algorithm

3.7 The Harvey-van der Hoeven result (2019)

4. Comparison and Algorithm Switching

4.1 Complexity comparison

4.2 Automatic algorithm switching

/* Example from mpn/x86_64/gmp-mparam.h */
#define MUL_TOOM22_THRESHOLD     20   /* schoolbook -> 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 Asymmetric multiplication

Visualization of FFT multiplication

5. GPU Parallelism and Multiprecision Multiplication

A GPU (Graphics Processing Unit) is a processor specialized in massively parallel processing on thousands to tens of thousands of cores. The GPU has delivered revolutionary performance for deep learning and scientific computing — can it also be applied to multiprecision computation? This section examines the possibilities and challenges of GPU-based multiprecision multiplication, taking the GPU's architectural characteristics into account.

5.1 GPU architecture and the SIMT model

NVIDIA GPUs are based on the SIMT (Single Instruction, Multiple Threads) model. A warp of 32 threads executes the same instruction simultaneously, and thousands of threads on an SM (Streaming Multiprocessor) cooperate. This structure is best suited to problems with high data parallelism — where the same operation is applied to a large amount of data.

Two primary approaches exist for porting multiprecision computation to a GPU:

1. One number = many threads: parallelize a single huge multiprecision number's operations across many threads. FFT/NTT multiplication is the canonical example.
2. One number = one thread: each thread handles an independent moderate-precision multiprecision number (a few hundred to a few thousand bits). Monte Carlo methods and cryptographic operations are typical.

5.2 FFT/NTT multiplication on GPU

FFT is intrinsically data-parallel and a good fit for the GPU. The cuFFT library executes double-precision FFT very fast and can serve as the foundation for multiprecision multiplication.

For multiprecision integer multiplication, however, NTT (Number-Theoretic Transform) is more appropriate. NTT is an FFT over the finite field $\mathbb{F}_p$, with no floating-point rounding errors. The basic structure of an NTT implementation on the GPU is:

// Sketch of NTT on GPU (CUDA pseudocode)
__global__ void ntt_kernel(uint64_t* a, int n, uint64_t w, uint64_t p) {
    // Run butterflies in parallel across threads
    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();
    }
}

The performance of NTT multiplication on the GPU depends heavily on memory access patterns. At each butterfly stage, the access stride changes as 1, 2, 4, ..., so early stages have non-contiguous access and are prone to bank conflicts. Effective use of shared memory and careful data layout are key to performance.

5.3 Challenges of multiprecision computation on GPUs

GPU multiprecision computation faces several intrinsic challenges arising from the architecture.

(a) Sequential nature of carry propagation

The carry of multiprecision addition is intrinsically sequential. In the worst case, the carry propagates through all $n$ digits (e.g. $999\ldots9 + 1 = 1000\ldots0$). On a CPU this is a simple loop; in the GPU's parallel model it is a poor fit.

The main mitigations are:

• Redundant representation: leave slack in each limb to defer carry propagation; only the final normalization is sequential.
• Parallel prefix addition: Brent-Kung or Kogge-Stone algorithms determine all carries in $O(\log n)$ stages.
• Carry-lookahead principle: compute Generate/Propagate bits in parallel and resolve carries in $O(\log n)$.

(b) Warp divergence

In the SIMT model, performance drops when threads of a warp take different branches (warp divergence). In multiprecision computation, this penalty arises when number sizes are uneven or conditional branches occur. For example, GCD computation has a data-dependent termination condition, making GPU parallelization difficult.

(c) Memory bandwidth bottleneck

Multiprecision multiplication has relatively low arithmetic intensity. An $n$-limb multiplication performs $O(n \log n)$ operations and $O(n)$ memory accesses, giving an arithmetic intensity of about $O(\log n)$. Since GPU compute throughput far exceeds memory bandwidth (about 5:1 byte/FLOP on the A100), small to medium multiprecision numbers tend to be memory-bound.

5.4 Existing libraries and implementations

Several research projects and libraries exist for GPU multiprecision computation.

5.5 Performance characteristics and applicability

Where the GPU has an advantage over the CPU, and where it does not.

5.6 Future outlook for GPU multiprecision computation

At present, GPUs have not reached the stage of fully replacing CPU-based multiprecision libraries (GMP, MPFR). Nevertheless, the following technological advances steadily expand the GPU's domain of applicability.

• Unified Memory: reduces CPU-GPU transfer overhead, making medium-size multiprecision computations natural to offload.
• Warp-level primitives: instructions like __shfl_sync exchange data inside a warp, enabling fast carry propagation across registers.
• Tensor Core applications: research uses Tensor Cores (specialized for matrix multiplication) for partial-product computation in integer multiplication. INT8/INT4 Tensor Cores are available.
• Chiplets and multi-GPU: NVLink/NVSwitch high-speed inter-GPU communication allows distributing very large FFTs across multiple GPUs.

In conclusion, GPU multiprecision computation greatly outperforms CPUs in the niches of "many independent operations" or "a single very large multiplication". On the other hand, efficiently realizing all features of a general-purpose multiprecision library on the GPU remains difficult because of the intrinsic sequentiality of carry propagation and warp divergence. The most realistic approach is the CPU-GPU hybrid — the CPU handles sequential control and small-scale operations while large multiplications and large batches of independent operations are offloaded to the GPU.

6. Distributed Parallelism (MPI/OpenMP)

The previous section's GPU parallelism exploited intra-node parallelism. For computations of billions to trillions of digits (e.g. world-record $\pi$ runs), distributed parallelism across multiple nodes is essential. This section surveys the use of MPI (Message Passing Interface) and OpenMP in multiprecision computation.

6.1 Thread parallelism with OpenMP

OpenMP is a multithreading framework for shared-memory environments. Where OpenMP is effective in multiprecision computation:

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);  /* fall back to schoolbook */
        return;
    }

    /* Split: a = a1 * B + a0, b = b1 * B + b0 */
    size_t half = n / 2;
    mpz_t a0, a1, b0, b1, z0, z1, z2;
    /* ... split processing ... */

    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 {
        /* fall back to sequential execution */
        karatsuba_parallel(z0, a0, b0, depth + 1);
        karatsuba_parallel(z2, a1, b1, depth + 1);
        /* ... */
    }

    /* result = z2 * B^2 + z1 * B + z0 */
    /* ... combine ... */
}

6.2 Distributed-memory parallelism with MPI

For computations larger than a single machine's memory or to accelerate across many nodes, distributed-memory parallelism with MPI is required. The main MPI patterns for multiprecision multiplication:

• Distributed FFT/NTT: split the huge limb arrays across nodes and combine FFT butterflies with inter-node communication (MPI_Alltoall). Per-stage communication is $O(N/P)$ ($P$ being the number of processes), but all-to-all latency tends to dominate.
• CRT-based distributed multiplication: assign different NTT primes to different nodes, run NTT multiplication independently within each node, and combine via CRT. Communication is needed only at the final combine, giving very high communication efficiency.
• Distributed Binary Splitting: assign high-level branches of the recursion tree to different nodes and combine the results in the final stage.

6.3 y-cruncher's parallelization strategy

y-cruncher, which has set many world records in $\pi$ computation, is a practical model of parallelization for multiprecision computation. It uses the following multi-level parallelism.

• SIMD: AVX-512 vectorization of NTT butterflies
• Multithreading: core-level parallelism via OpenMP / a custom thread pool
• Disk swapping: data that does not fit in memory is paged to disk and processed using out-of-core algorithms
• NUMA optimization: data placement aware of memory affinity

Summary

This chapter studied the main algorithms for fast multiplication of multiprecision integers.

• Karatsuba algorithm (1962): 2-way split, 3 multiplications, $O(n^{1.585})$. The most fundamental fast multiplication, reducing $n$-digit multiplication to three $n/2$-digit multiplications. The idea of cutting from 4 to 3 multiplications relies on obtaining the cross term $A_0 C_1 + A_1 C_0$ from the single multiplication $(A_0 + A_1)(C_0 + C_1) - A_0 C_0 - A_1 C_1$.
• Toom-Cook algorithm: $k$-way split with polynomial interpolation, $O(n^{\log_k(2k-1)})$. A generalization of Karatsuba understandable in the framework of "evaluation → pointwise multiplication → interpolation → recomposition." Larger $k$ lowers the exponent, but the interpolation overhead grows; in practice Toom-3 and Toom-4 are the limit.
• FFT multiplication: by the convolution theorem, multiplication reduces to pointwise multiplication in the frequency domain, computing the convolution in $O(n \log n)$.
• NTT (number-theoretic transform): an FFT over a finite field that completely eliminates rounding error. Essential for integer multiplication.
• Schönhage-Strassen (1971): an FFT over the Fermat number ring achieving $O(n \log n \log \log n)$. Core of GMP's large multiplication.
• Harvey-van der Hoeven (2019): the theoretically optimal $O(n \log n)$, but with such a large constant factor that it is not practical.
• Algorithm switching: automatically picking the optimal algorithm based on input size is critically important in the design of practical multiprecision libraries.
• GPU parallelism: GPUs greatly outperform CPUs on large batches of independent multiprecision operations and on very large NTT multiplications, but generic application remains challenging due to sequential carry propagation and warp divergence. A CPU-GPU hybrid is the realistic answer.
• Distributed parallelism: at the trillion-digit scale, MPI-based distributed FFT/NTT and CRT-based distributed multiplication are essential. OpenMP thread parallelism is effective for the recursive splits of Karatsuba and for Binary Splitting; y-cruncher applies a multilayered combination of SIMD, multithreading, disk swapping, and NUMA awareness.

The next chapter builds on these multiplication algorithms to handle advanced multiprecision integer operations: GCD, modular arithmetic, modular exponentiation, and primality testing.