Harvey and Hoeven O(n log n) Integer Multiplication

In 2019, mathematicians David Harvey and Joris van der Hoeven published a breakthrough algorithm achieving the long-conjectured \(O(n \log n)\) time complexity for multiplying two \(n\)-bit integers. This article explains how their method processes arbitrary-length binary numbers by dividing them into structured bit blocks, mapping them into multidimensional polynomial rings, computing their convolution using higher-dimensional Fast Fourier Transforms (FFTs), and propagating carries to reconstruct the final binary product.

1. Partitioning Binary Integers into Polynomials

The algorithm begins by breaking two large binary integers, \(A\) and \(B\), of length \(n\) bits, into smaller chunks. Given a chosen block size \(\beta\), the integers are represented in base \(2^\beta\):

\[A = \sum_{i=0}^{k-1} a_i 2^{i\beta}, \quad B = \sum_{i=0}^{k-1} b_i 2^{i\beta}\]

where each coefficient \(a_i, b_i \in [0, 2^\beta - 1]\) fits into a \(\beta\)-bit register.

These representations are mapped directly to one-dimensional polynomials:

\[A(x) = \sum_{i=0}^{k-1} a_i x^i, \quad B(x) = \sum_{i=0}^{k-1} b_i x^i\]

Computing the integer product \(A \times B\) is equivalent to evaluating the polynomial product \(C(x) = A(x)B(x)\) at \(x = 2^\beta\) and resolving any carries where coefficients exceed \(2^\beta - 1\).

2. Moving to Multidimensional Transforms

Standard fast multiplication methods, such as the Schönhage–Strassen algorithm (\(O(n \log n \log \log n)\)) or Fürer’s algorithm (\(O(n \log n \cdot 2^{O(\log^* n)})\)), compute one-dimensional FFTs. These approaches suffer from coefficient growth and recursive overhead when evaluating large roots of unity.

Harvey and Hoeven eliminate this bottleneck by restructuring the one-dimensional polynomial multiplication into a multidimensional convolution. The single index \(i\) of polynomial coefficients is mapped to a \(d\)-dimensional grid:

\[(i_1, i_2, \dots, i_d) \in \mathbb{Z}_{m_1} \times \mathbb{Z}_{m_2} \times \dots \times \mathbb{Z}_{m_d}\]

This transforms the problem into multiplying multidimensional polynomials over a ring of truncated polynomials or Gaussian integers:

\[R = \mathbb{C}[y] / (y^M + 1)\]

By distributing the total transform length across \(d\) dimensions, each individual dimension remains small. This keeps the precision required for the roots of unity bounded and prevents the progressive blowup of coefficients across recursive stages.

3. Evaluating via Multidimensional FFTs

With the inputs arranged across a multidimensional lattice, the algorithm performs the following steps:

  1. Forward Transform: Apply an FFT along each of the \(d\) dimensions independently to convert spatial coefficient arrays into the frequency domain.
  2. Pointwise Multiplication: Multiply the transformed arrays pointwise. Because the transform domain maps convolution to standard multiplication, this requires multiplying smaller elements in the underlying coefficient ring \(R\).
  3. Inverse Transform: Apply the multidimensional Inverse Fast Fourier Transform (IFFT) to return the result to the coefficient domain.

The choice of ring \(R\) ensures that the arithmetic operations inside the FFT can themselves be computed using the same multiplication algorithm recursively, maintaining the \(O(n \log n)\) bound throughout the recursion tree.

4. Carry Propagation and Binary Reconstruction

After computing the inverse transform, the result is a sequence of integer coefficients \(c_k\) representing \(C(x) = \sum c_k x^k\). Because the transform operates on exact coefficients, each \(c_k\) is the exact sum of products:

\[c_k = \sum_{i+j=k} a_i b_j\]

Each coefficient \(c_k\) can span up to \(2\beta + \lceil \log_2 k \rceil\) bits, which is larger than the original chunk size \(\beta\). To produce the final binary number:

  1. The coefficients are placed into memory shifted by their corresponding powers of two: \(c_k \cdot 2^{k\beta}\).
  2. A linear-time (\(O(n)\)) carry-propagation pass scans the overlapped segments from lowest bit to highest bit, adding overflow values to higher-order positions.

Once carries are resolved, the output is the exact \(2n\)-bit product \(A \times B\) in standard binary format.