Divide & Conquer/Fast Multiplication

Lesson 2.43,415 words

Fast Multiplication

Grade-school multiplication is Θ(n2)\Theta(n^2), yet divide and conquer beats it. Karatsuba multiplies nn-bit integers with three half-size products instead of four, giving Θ(nlog23)\Theta(n^{\log_2 3}), and Strassen multiplies matrices with seven block products instead of eight, giving Θ(nlog27)\Theta(n^{\log_2 7}).

╌╌╌╌

Adding two -bit numbers is easy: the grade-school ripple-carry method is , and you cannot beat linear because you must at least read the input. Multiplication is the interesting one. The grade-school algorithm forms shifted partial products and sums them, costing . For a long time that was assumed to be the best possible — until divide and conquer showed otherwise. The same idea then breaks the cubic bound for matrix multiplication. Both stories share a single idea: spend a few cheap additions to buy back one expensive multiplication.1

Multiplying integers: the schoolbook baseline

Write each -bit integer as a high half and a low half. With , split into its least-significant bits and its top bits , and likewise split into and :

Splitting each -bit integer and into high and low -bit halves, so and .

Multiplying out and gives

The naive split forms all four cross-products of the halves : the diagonal terms and become the low and high coefficients, while and both feed the middle coefficient . Four half-size multiplications in all.

The powers of two are just bit-shifts (free, up to the additions), so the cost is dominated by the four half-size products , , , . That yields a recurrence

Unrolling it shows this naïve split buys nothing:

At the leaf term is , which dominates everything: . The master theorem says the same instantly: with , , , the branching exponent exceeds the work exponent , so the recurrence is leaf-heavy and . Four subproblems of half size land us right back in the quadratic regime.

karatsuba.pypython
def schoolbook_multiply(first: int, second: int) -> int:
  """
    The naive split into four half-size products, kept as the baseline.\n
    Writing x = a + 2^t * b and y = c + 2^t * d gives\n
    xy = ac + 2^t (ad + bc) + 2^(2t) bd, four recursive products, which the\n
    master theorem pins at Theta(n^2) — no better than grade school.\n
  """
  # factor out the sign so the recursion runs on magnitudes.
  sign: int = 1
  if first < 0:
    sign, first = -sign, -first
  if second < 0:
    sign, second = -sign, -second

  def recurse(left: int, right: int) -> int:
    # single-word operands multiply directly.
    if left < 2 or right < 2:
      return left * right

    # split each operand at the midpoint of the wider one's bits.
    half_bits: int = max(left.bit_length(), right.bit_length()) // 2
    shift: int = 1 << half_bits
    low_left, high_left = left & (shift - 1), left >> half_bits
    low_right, high_right = right & (shift - 1), right >> half_bits

    # all four cross-products of the halves — the wasteful part.
    low_product: int = recurse(low_left, low_right)
    high_product: int = recurse(high_left, high_right)
    middle: int = recurse(low_left, high_right) + recurse(high_left, low_right)

    # recombine: high * 2^2t + middle * 2^t + low.
    return (high_product << (2 * half_bits)) + (middle << half_bits) + low_product

  return sign * recurse(first, second)

Karatsuba: three products instead of four

The fix, due to Karatsuba (1960), is to compute the middle coefficient without a third and fourth multiplication.2 Notice the algebraic identity

We were going to compute and anyway. So once we have those two products, the middle term needs only one more multiplication, , plus a few additions. Three half-size products now suffice where four were needed.

Karatsuba reuses the products and to recover the middle coefficient from a single extra product (the solid blue box), replacing the naive four products with three.

The low and high coefficients are just and ; the middle coefficient is recovered from the single extra product plus the two reused products, as .

Algorithm 1:Karatsuba(x,y)\textsc{Karatsuba}(x, y) — multiply two nn-bit integers
  1. 1
    if n=1n = 1 then
  2. 2
    return xyx \cdot y
    base case: one-bit (or word) product
  3. 3
    tn/2t \gets \floor{n / 2}
  4. 4
    split x=a+2tbx = a + 2^{t} b and y=c+2tdy = c + 2^{t} d
    low and high halves
  5. 5
    acKaratsuba(a,c)\mathit{ac} \gets \textsc{Karatsuba}(a, c)
    recursive product 1
  6. 6
    bdKaratsuba(b,d)\mathit{bd} \gets \textsc{Karatsuba}(b, d)
    recursive product 2
  7. 7
    mKaratsuba(ab,dc)m \gets \textsc{Karatsuba}(a - b, d - c)
    recursive product 3 — the trick
  8. 8
    midm+ac+bd\mathit{mid} \gets m + \mathit{ac} + \mathit{bd}
    recovers ad+bcad + bc
  9. 9
    return ac+2tmid+22tbd\mathit{ac} + 2^{t}\,\mathit{mid} + 2^{2t}\,\mathit{bd}
    shift and add
karatsuba.pypython
def karatsuba(first: int, second: int) -> int:
  """
    Product of two integers via Karatsuba's three-product recursion.\n
    Handles negative operands by factoring out the sign; the recursion runs\n
    on magnitudes, splitting at the midpoint of the wider operand's bits.\n
  """
  # factor out the sign, then recurse on the magnitudes.
  sign: int = 1
  if first < 0:
    sign, first = -sign, -first
  if second < 0:
    sign, second = -sign, -second

  return sign * _karatsuba_magnitude(first, second)


def _karatsuba_magnitude(left: int, right: int) -> int:
  """
    Karatsuba on non-negative magnitudes.\n
    Base case multiplies single-word operands directly; otherwise it forms\n
    the three products ac, bd, and (a - b)(d - c) and recovers the middle\n
    coefficient as (a - b)(d - c) + ac + bd = ad + bc before shifting.\n
  """
  if left < 2 or right < 2:
    return left * right

  # split at the midpoint of the wider operand so neither half vanishes.
  half_bits: int = max(left.bit_length(), right.bit_length()) // 2
  shift: int = 1 << half_bits
  low_left, high_left = left & (shift - 1), left >> half_bits
  low_right, high_right = right & (shift - 1), right >> half_bits

  low_product: int = _karatsuba_magnitude(low_left, low_right)        # a*c
  high_product: int = _karatsuba_magnitude(high_left, high_right)     # b*d

  # the single extra product is the whole trick: it lets us skip a fourth call.
  cross: int = _karatsuba_magnitude(
    abs(low_left - high_left),
    abs(low_right - high_right),
  )

  # sign of (a - b)(d - c): the magnitude above dropped it, so restore it.
  cross_sign: int = 1
  if (low_left < high_left) != (high_right < low_right):
    cross_sign = -1
  middle: int = cross_sign * cross + low_product + high_product       # = ad + bc

  return (high_product << (2 * half_bits)) + (middle << half_bits) + low_product

A worked example in base 10

The identity is base-agnostic, so it is easiest to watch on ordinary decimal digits. Multiply by . Split each at the middle two digits, using as the base (so counts digits here, not bits):

Karatsuba needs three products of the two-digit halves:

The middle coefficient falls out of the identity with no further multiplication:

Reassemble by shifting each coefficient to its place value and adding:

Karatsuba on in base . Three two-digit products feed the low, middle, and high coefficients; each is shifted to its place value and summed to give .

A cross-check confirms it: . The schoolbook method would have formed four two-digit products (, , , ); Karatsuba used three and recovered the fourth's contribution from a subtraction.

Cost

There are now recursive calls, each on half-size inputs, with work to split, shift, and add (the subtractions , may produce one extra bit, which does not change the asymptotics):

The recursion tree makes the cost legible. Level has nodes, each doing work, for a row total of . The rows grow geometrically downward, so the leaves dominate:

Each Karatsuba node spawns half-size children, so the rows grow by downward and the leaf level dominates the cost.

Summing the geometric row totals,

using (the identity ). The master theorem confirms it in one line: with , , , the branching exponent beats the work exponent , so this is leaf-heavy and . Trading one multiplication for a handful of additions drops integer multiplication below quadratic.

Strassen: seven products instead of eight

The same spend additions to save a multiplication idea breaks the cubic bound for matrix multiplication. The schoolbook method costs : each of the entries of is a length- dot product. Divide and conquer splits each matrix into four blocks, and the product is read off block by block exactly as for scalars:

Each matrix splits into four blocks. Naively, every output block is a sum of two block products, so the four outputs cost eight half-size multiplications.

Taken at face value, the four output blocks need eight block multiplications, giving . The master theorem returns : more subproblems exactly cancel their smaller size, so blocking alone buys nothing — the integer-multiplication story again, one dimension up.

Strassen's 1969 insight is that seven products suffice.3 Form seven recursive -size multiplications,

and recover the output blocks by addition only:

The eighth multiplication is gone, traded for a handful of extra block additions that cost only . The recurrence becomes

and now the leaves win. With , , , the branching exponent exceeds the work exponent , so by the master theorem

Naive blocking fans out into half-size products (root-heavy ties at ); Strassen drops one to , and the smaller fan-out — not the block size — pulls the exponent down to .

The constant is large, so the crossover with the cubic method only pays off for sizable ; in practice Strassen is switched in above a threshold and the base case falls back to the cache-friendly schoolbook multiply. Theoretically, though, it was the first sub-cubic algorithm, and a long line of refinements has pushed the exponent down toward .

strassen.pypython
from typing import Sequence

Matrix = list[list[int]]

def naive_multiply(left: Sequence[Sequence[int]], right: Sequence[Sequence[int]]) -> Matrix:
  """
    Schoolbook matrix product — each of the n^2 entries is a dot product.\n
    The Theta(n^3) baseline Strassen is measured against, and the base case\n
    real implementations fall back to below the crossover size.\n
  """
  rows: int = len(left)
  inner: int = len(right)
  columns: int = len(right[0]) if right else 0
  product: Matrix = [[0 for _ in range(columns)] for _ in range(rows)]

  for row in range(rows):
    for index in range(inner):
      # skip zero multipliers; they contribute nothing to the accumulator.
      multiplier: int = left[row][index]
      if multiplier == 0:
        continue

      # fold this term of the dot product into every column of the output row.
      left_row, right_row = product[row], right[index]
      for column in range(columns):
        left_row[column] += multiplier * right_row[column]

  return product

def _add(left: Matrix, right: Matrix) -> Matrix:
  """
    Entrywise sum of two equally shaped square blocks.\n
  """
  size: int = len(left)
  return [[left[row][col] + right[row][col] for col in range(size)] for row in range(size)]

def _subtract(left: Matrix, right: Matrix) -> Matrix:
  """
    Entrywise difference of two equally shaped square blocks.\n
  """
  size: int = len(left)
  return [[left[row][col] - right[row][col] for col in range(size)] for row in range(size)]

def _split_quadrants(matrix: Matrix) -> tuple[Matrix, Matrix, Matrix, Matrix]:
  """
    Split an even-sized square matrix into its four corner quadrants\n
    (top-left, top-right, bottom-left, bottom-right).\n
  """
  half: int = len(matrix) // 2
  top_left: Matrix = [row[:half] for row in matrix[:half]]
  top_right: Matrix = [row[half:] for row in matrix[:half]]
  bottom_left: Matrix = [row[:half] for row in matrix[half:]]
  bottom_right: Matrix = [row[half:] for row in matrix[half:]]
  return top_left, top_right, bottom_left, bottom_right

def _join_quadrants(
  top_left: Matrix, top_right: Matrix, bottom_left: Matrix, bottom_right: Matrix
) -> Matrix:
  """
    Reassemble four quadrants into one square matrix.\n
  """
  top: Matrix = [
    left_row + right_row for left_row, right_row in zip(top_left, top_right)
  ]
  bottom: Matrix = [
    left_row + right_row for left_row, right_row in zip(bottom_left, bottom_right)
  ]
  return top + bottom

def _pad_to_power_of_two(matrix: Matrix) -> Matrix:
  """
    Zero-pad a square matrix up to the next power-of-two side length so the\n
    halving recursion always divides evenly. Padding rows and columns of\n
    zeros do not change the product's top-left block.\n
  """
  # find the next power-of-two side length at or above the current size.
  size: int = len(matrix)
  target: int = 1
  while target < size:
    target *= 2

  # already a power of two: hand back a fresh copy, nothing to pad.
  if target == size:
    return [list(row) for row in matrix]

  # copy the original into the top-left block of a zero-filled square.
  padded: Matrix = [[0 for _ in range(target)] for _ in range(target)]
  for row in range(size):
    for column in range(size):
      padded[row][column] = matrix[row][column]

  return padded

def _strassen_square(left: Matrix, right: Matrix, leaf_size: int) -> Matrix:
  """
    Strassen's recursion on two equal power-of-two square matrices.\n
    Below `leaf_size` it falls back to the schoolbook multiply, matching the\n
    crossover strategy of real libraries.\n
  """
  # below the crossover, the schoolbook multiply wins on overhead.
  size: int = len(left)
  if size <= leaf_size:
    return naive_multiply(left, right)

  # break both operands into their four corner quadrants.
  a11, a12, a21, a22 = _split_quadrants(left)
  b11, b12, b21, b22 = _split_quadrants(right)

  # seven products replace the naive eight.
  m1: Matrix = _strassen_square(_add(a11, a22), _add(b11, b22), leaf_size)
  m2: Matrix = _strassen_square(_add(a21, a22), b11, leaf_size)
  m3: Matrix = _strassen_square(a11, _subtract(b12, b22), leaf_size)
  m4: Matrix = _strassen_square(a22, _subtract(b21, b11), leaf_size)
  m5: Matrix = _strassen_square(_add(a11, a12), b22, leaf_size)
  m6: Matrix = _strassen_square(_subtract(a21, a11), _add(b11, b12), leaf_size)
  m7: Matrix = _strassen_square(_subtract(a12, a22), _add(b21, b22), leaf_size)

  # recover the four output blocks by addition alone.
  c11: Matrix = _add(_subtract(_add(m1, m4), m5), m7)
  c12: Matrix = _add(m3, m5)
  c21: Matrix = _add(m2, m4)
  c22: Matrix = _add(_add(_subtract(m1, m2), m3), m6)

  return _join_quadrants(c11, c12, c21, c22)

def strassen(left: Sequence[Sequence[int]], right: Sequence[Sequence[int]],
             leaf_size: int = 1) -> Matrix:
  """
    Product of two conformable integer matrices via Strassen's algorithm.\n
    Rectangular and non-power-of-two inputs are zero-padded to a common\n
    power-of-two square, multiplied, then trimmed back to the true shape.\n
    `leaf_size` sets the crossover below which the schoolbook multiply runs.\n
  """
  # an empty dimension means an empty (or zero-filled) product.
  rows: int = len(left)
  inner: int = len(right)
  columns: int = len(right[0]) if right else 0
  if rows == 0 or inner == 0 or columns == 0:
    return [[0 for _ in range(columns)] for _ in range(rows)]

  # pad both operands into a common power-of-two square so halving divides evenly.
  side: int = max(rows, inner, columns)
  square_left: Matrix = _pad_to_power_of_two(
    [list(row) + [0 for _ in range(side - inner)] for row in left]
    + [[0 for _ in range(side)] for _ in range(side - rows)]
  )
  square_right: Matrix = _pad_to_power_of_two(
    [list(row) + [0 for _ in range(side - columns)] for row in right]
    + [[0 for _ in range(side)] for _ in range(side - inner)]
  )

  # multiply the squares, then trim back to the true output shape.
  square_product: Matrix = _strassen_square(square_left, square_right, max(1, leaf_size))
  return [square_product[row][:columns] for row in range(rows)]

Crossover: when divide and conquer actually wins

Both algorithms carry a fat constant: Karatsuba does extra additions and subtractions, Strassen does block additions per level versus the naive products' simpler bookkeeping. Below some threshold the schoolbook method's smaller constant wins outright, so every real implementation switches back to schoolbook at the base case rather than recursing down to .

The fast method has a better exponent but a fatter constant, so its cost curve starts above schoolbook and crosses it only at a threshold . Below schoolbook is cheaper (left of the dashed line); above it the fast method wins. Libraries switch over only past the crossover.

This is the recurring theme of asymptotically fast algorithms — a better exponent is worthless until is large enough to overcome the constant, so the clever method is layered on top of a simple base case, not used in place of it.

Toom–Cook: splitting into more pieces

Karatsuba splits each operand into two halves and uses three products. The same idea generalizes. View an -digit integer written in base as a degree- polynomial whose coefficients are the chunks: is the polynomial evaluated at . Multiplying the integers is multiplying the polynomials and then substituting (a shift-and-add that resolves carries).

The product has degree , so it is determined by its values at points. Toom–Cook (specifically Toom-) runs three phases:

  • evaluate and at small points (e.g. ),
  • multiply the paired values pointwise — recursive products of -size numbers,
  • interpolate the products back into the coefficients of , then evaluate at .

Karatsuba is Toom-: two chunks, products. Toom- uses three chunks and products of third-size operands, giving — a better exponent than Karatsuba's , at the price of a messier interpolation step with larger additive constants. In general Toom- achieves , and as grows the exponent creeps toward — but the evaluate/interpolate overhead grows too, so no fixed reaches near-linear.

FFT multiplication: near-linear

The Fast Fourier Transform takes the polynomial view to its conclusion. Instead of a fixed handful of evaluation points, evaluate at the complex -th roots of unity. The FFT performs that evaluation for all roots at once in (a divide-and-conquer that splits a polynomial into its even- and odd-indexed coefficients), the pointwise multiply is , and the inverse FFT interpolates back in another . Substituting and resolving carries finishes the integer product. The total is

for the polynomial multiply, and for integer multiplication once one accounts for the growing precision of the roots (Schönhage–Strassen). A 2019 result of Harvey and van der Hoeven removed the last factor, reaching — conjectured to be optimal.

The through-line is one substitution: multiplying two integers is multiplying two polynomials in the base variable, then carrying. Every method here — the schoolbook grid, Karatsuba, Toom–Cook, FFT — is a different way to multiply those polynomials, trading more evaluation-and-interpolation bookkeeping for fewer recursive products. Fast multiplication therefore leads directly to fast polynomial arithmetic, and shows up wherever big integers or high-degree polynomials do: cryptography, computer algebra, and signal processing among them.4

The exponent ladder for integer multiplication. Splitting into more chunks and reusing products drops the exponent from (schoolbook) toward (FFT); each rung trades heavier evaluate/interpolate overhead for fewer recursive products.

Frontiers of multiplication

The two stories in this lesson — integer and matrix multiplication — both remain open at their frontiers, and both illustrate a gap between what is asymptotically fastest and what any machine will ever run.

Integer multiplication: the conjectured floor was reached. For decades the best known bound was Schönhage–Strassen's (1971), the FFT method of the previous section. It was long conjectured that is the true complexity, with no room below. In 2019 Harvey and van der Hoeven proved exactly that bound: an algorithm multiplying -bit integers in . It is a landmark — it likely closes the problem — but it is a galactic algorithm: the constant hidden in the is so enormous (the method lifts the numbers into a high-dimensional FFT over a cleverly chosen ring) that the crossover where it beats Schönhage–Strassen exceeds any input that will ever be multiplied. It settles the theory while changing nothing in practice, where the tower of schoolbook, Karatsuba, Toom–Cook, and Schönhage–Strassen from earlier still governs which method a library picks by operand size.

Matrix multiplication: a moving target. Strassen's opened a race that is still running. Let be the matrix-multiplication exponent — the smallest number such that matrices can be multiplied in for every . Strassen put ; the Coppersmith–Winograd algorithm (1990) and its laser-method descendants (Williams 2012, and Alman & Williams 2021) have pushed the record to . Every one of these is galactic — the constants make them useless below astronomically large — so practical libraries still use Strassen above a tuned threshold and schoolbook below it. The lower bound is only the trivial (you must at least read the inputs), and whether is one of the central open questions of algorithm theory.

The matrix-multiplication exponent over time. Strassen's began a descent toward the trivial floor ; every improvement past Strassen is galactic (useless at practical sizes), and whether is open.

When the algebra is discovered by search. Strassen found his seven products by hand; the modern twist is to let a machine search for such identities. AlphaTensor (Fawzi et al., 2022) cast the discovery of low-rank matrix multiplication schemes as a single-player game and, with reinforcement learning, rediscovered Strassen-like decompositions and found new ones for specific small sizes (for instance a scheme over using fewer multiplications than the previously known best). It is the same idea this lesson turns on — trade multiplications for additions — but with the search for which additions automated rather than reasoned out. The frontier of spend cheap operations to save expensive ones is now partly a search problem.

Takeaways

  • Schoolbook integer multiplication is ; the naive divide split into four half-size products, , stays .
  • Karatsuba uses the identity to compute the middle coefficient from one extra product, cutting four subproblems to three: .
  • Schoolbook matrix multiplication is ; naive block divide, , stays .
  • Strassen recovers the four output blocks from seven products instead of eight: — the first sub-cubic matrix multiply.
  • The master theorem reads off every one of these by comparing the branching exponent to the work exponent ; both fast methods are leaf-heavy.
  • Toom–Cook generalizes Karatsuba by splitting into chunks and treating the digits as polynomial coefficients: Toom- uses products for ; Karatsuba is Toom-.
  • The FFT evaluates at roots of unity to multiply the underlying polynomials in , giving near-linear integer multiplication (down to since 2019).
  • Both win only asymptotically: a fat constant means the clever method is layered above a schoolbook base case, switched in only past a crossover size .
  • At the frontier both problems are open or galactic: Harvey–van der Hoeven (2019) reached the conjectured for integers, and the matrix-multiplication exponent has fallen to (Alman–Williams) with the trivial the only known floor — all useless at real sizes, and now partly discovered by machine search (AlphaTensor).

Footnotes

  1. CLRS, Ch. 4 — Divide-and-Conquer: integer and matrix multiplication as divide-and-conquer case studies, trading multiplications for additions.
  2. Erickson, Algorithms, Ch. 1 — Recursion: Karatsuba's -bit integer multiplication via the identity reducing four half-size products to three, with the recurrence.
  3. CLRS, Ch. 4 (§4.2) — Strassen's algorithm for matrix multiplication and its recurrence; the original result is V. Strassen (1969), Gaussian elimination is not optimal, Numerische Mathematik 13, 354–356.
  4. Skiena, The Algorithm Design Manual, §13 — Numerical Problems: high-precision arithmetic, the link from integer to polynomial multiplication, and FFT-based methods for large operands.
Practice

╌╌ END ╌╌