Trading Expensive Multiplications for Cheap Additions
When we divide a problem of size $n$ into subproblems of size $n/2$, the Master Theorem tells us that for $T(n) = a T(n/2) + \Theta(n^d)$ with $a > 2^d$, the asymptotic complexity is governed purely by the number of leaves in the recursion tree: $\Theta(n^{\log_2 a})$. Even reducing $a$ by one single recursive call changes the exponent of the polynomial!
1. Karatsuba Fast Integer & Polynomial Multiplication (1960)
Given two $n$-limb big integers (or degree-$n$ polynomials) $X$ and $Y$, split each at $m = \lfloor n/2 \rfloor$ into high and low halves in base $B^m$:
$$X = x_1 B^m + x_0, \qquad Y = y_1 B^m + y_0$$Their product is:
$$X \cdot Y = \underbrace{x_1 y_1}_{z_2} B^{2m} + \underbrace{(x_1 y_0 + x_0 y_1)}_{z_1} B^m + \underbrace{x_0 y_0}_{z_0}$$Naively computing $x_1 y_1, x_1 y_0, x_0 y_1, x_0 y_0$ requires $a = 4$ multiplications of size $n/2$, giving $T(n) = 4T(n/2) + \Theta(n) = \Theta(n^{\log_2 4}) = \Theta(n^2)$. Anatolii Karatsuba observed (following Gauss's complex multiplication trick) that $z_1$ can be recovered using only one additional multiplication plus two subtractions:
$$z_2 = x_1 y_1, \qquad z_0 = x_0 y_0, \qquad z_1 = (x_1 + x_0)(y_1 + y_0) - z_2 - z_0$$Now $a = 3$, so $T(n) = 3T(n/2) + \Theta(n) = \Theta\!\left(n^{\log_2 3}\right) \approx \Theta\!\left(n^{1.58496}\right)$.
2. Strassen's Sub-Cubic Matrix Multiplication (1969)
Given two $n \times n$ matrices $A$ and $B$, partition each into four $\frac{n}{2} \times \frac{n}{2}$ quadrant blocks:
$$\begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix} = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix} \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}$$Standard block multiplication requires $8$ recursive $\frac{n}{2} \times \frac{n}{2}$ matrix products, giving $T(n) = 8T(n/2) + \Theta(n^2) = \Theta(n^3)$. Volker Strassen shocked the mathematical community by defining 7 intermediate block products $M_1, \dots, M_7$ using $18$ matrix additions/subtractions:
- $M_1 = (A_{11} + A_{22})(B_{11} + B_{22})$
- $M_2 = (A_{21} + A_{22}) B_{11}$
- $M_3 = A_{11} (B_{12} - B_{22})$
- $M_4 = A_{22} (B_{21} - B_{11})$
- $M_5 = (A_{11} + A_{12}) B_{22}$
- $M_6 = (A_{21} - A_{11})(B_{11} + B_{12})$
- $M_7 = (A_{12} - A_{22})(B_{21} + B_{22})$
From which the four output quadrants are reconstructed in $\Theta(n^2)$ additions:
$$C_{11} = M_1 + M_4 - M_5 + M_7, \quad C_{12} = M_3 + M_5, \quad C_{21} = M_2 + M_4, \quad C_{22} = M_1 - M_2 + M_3 + M_6$$This yields $T(n) = 7T(n/2) + \Theta(n^2) = \Theta\!\left(n^{\log_2 7}\right) \approx \Theta\!\left(n^{2.80735}\right)$.
Inspect the Radix-2 Cooley-Tukey FFT butterfly y_k = E_k + ω_n^k O_k, Karatsuba 3-multiplication algebraic identity, and CDQ Divide-and-Conquer left-to-right contribution transfer.
Invariant & State
| Algorithm | Recurrence | Asymptotic Time | Practical Crossover Threshold |
|---|---|---|---|
| Schoolbook Big-Int / Poly Multiply | $4T(n/2) + O(n)$ | $\Theta(n^2)$ | $n < 32$ limbs ($64$-bit words) |
| Karatsuba Multiplication | $3T(n/2) + O(n)$ | $\Theta(n^{1.585})$ | $32 \le n \le 2000$ limbs |
| Toom-Cook (Toom-3) / FFT Multiplication | $5T(n/3) + O(n)$ / FFT | $\Theta(n^{1.465})$ / $\Theta(n \log n)$ | $n > 2000$ limbs (GMP library) |
| Cache-Blocked GEMM vs. Strassen | $8T(n/2)$ vs. $7T(n/2) + 18(n/2)^2$ | $\Theta(n^3)$ vs. $\Theta(n^{2.807})$ | Strassen wins for $n \ge 512\text{–}1024$ |
In both Karatsuba ($a=3, b=2, d=1$) and Strassen ($a=7, b=2, d=2$), the branching factor $a$ strictly exceeds $b^d$ ($3 > 2^1$ and $7 > 2^2$). Consequently, work grows geometrically down the recursion tree by a factor of $3/2$ or $7/4$ per level. Stopping the recursion at a tuned base-case leaf size $n_0$ (e.g., $n_0 = 32$ for Karatsuba, $n_0 = 64$ for Strassen) replaces the bottom $\log_2 n_0$ levels of recursive allocation overhead with vectorizable SIMD loops, speeding up total runtime by $5\times\text{–}15\times$!
Two practical factors govern algebraic D&C in production:
- Heap Allocation Thrashing: If each Karatsuba call allocates three
std::vectors on the heap,malloc/freelock contention destroys performance. Passing a single pre-allocated contiguous scratch arena bump-pointer reduces runtime by over $4\times$. - Floating-Point Catastrophic Cancellation in Strassen: Because Strassen computes $C_{11} = M_1 + M_4 - M_5 + M_7$ by adding and subtracting quadrants of vastly different magnitudes, IEEE-754 floating-point error is bounded only in the normwise sense ($\|\Delta C\| \le c_n u \|A\| \|B\|$), not the componentwise sense ($|\Delta C_{ij}| \le c_n u (|A||B|)_{ij}$). If $A$ and $B$ have poorly scaled rows/columns, small entries in $C$ can lose all significant digits!
Production C++20 Implementation: Arena-Backed Karatsuba Polynomial Multiplier
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <vector>
// Multiplies two polynomials a[0..n-1] and b[0..n-1] (n = power of 2) into res[0..2n-1]
// Uses a pre-allocated scratch buffer to perform ZERO heap allocations during recursion!
void karatsuba_rec(const std::int64_t* a, const std::int64_t* b,
std::int64_t* res, std::int64_t* scratch, std::size_t n) noexcept {
std::fill_n(res, 2 * n, 0);
// Base case: Schoolbook O(n^2) loop auto-vectorizes cleanly for n <= 32
if (n <= 32) {
for (std::size_t i = 0; i < n; ++i) {
for (std::size_t j = 0; j < n; ++j) {
res[i + j] += a[i] * b[j];
}
}
return;
}
const std::size_t m = n / 2;
const std::int64_t* a0 = a;
const std::int64_t* a1 = a + m;
const std::int64_t* b0 = b;
const std::int64_t* b1 = b + m;
// Carve out workspace from linear scratch arena:
// sum_a[m], sum_b[m], z_m[2m], next_scratch[...]
std::int64_t* sum_a = scratch;
std::int64_t* sum_b = scratch + m;
std::int64_t* z_m = scratch + 2 * m;
std::int64_t* next_scratch = scratch + 4 * m;
for (std::size_t i = 0; i < m; ++i) {
sum_a[i] = a0[i] + a1[i];
sum_b[i] = b0[i] + b1[i];
}
// 1. z0 = a0 * b0 stored directly in res[0 .. 2m-1]
karatsuba_rec(a0, b0, res, next_scratch, m);
// 2. z2 = a1 * b1 stored directly in res[2m .. 4m-1]
karatsuba_rec(a1, b1, res + 2 * m, next_scratch, m);
// 3. z_m = (a0 + a1) * (b0 + b1) stored in z_m[0 .. 2m-1]
karatsuba_rec(sum_a, sum_b, z_m, next_scratch, m);
// Combine: z1 = z_m - z0 - z2 added into res[m .. 3m-1]
for (std::size_t i = 0; i < 2 * m; ++i) {
z_m[i] -= (res[i] + res[2 * m + i]);
}
for (std::size_t i = 0; i < 2 * m; ++i) {
res[m + i] += z_m[i];
}
}
std::vector<std::int64_t> karatsuba_multiply(const std::vector<std::int64_t>& A,
const std::vector<std::int64_t>& B) {
std::size_t n = 1;
const std::size_t max_len = std::max(A.size(), B.size());
while (n < max_len) n <<= 1;
std::vector<std::int64_t> a_pad(n, 0), b_pad(n, 0);
std::copy(A.begin(), A.end(), a_pad.begin());
std::copy(B.begin(), B.end(), b_pad.begin());
std::vector<std::int64_t> res(2 * n, 0);
std::vector<std::int64_t> scratch(6 * n, 0); // 4m + 4(m/2) + ... <= 4n scratch space
karatsuba_rec(a_pad.data(), b_pad.data(), res.data(), scratch.data(), n);
res.resize(A.size() + B.size() - 1);
return res;
}When multiplying polynomials, $a_0 + a_1$ has the exact same degree $m-1$ as $a_0$ and $a_1$. When multiplying base-$2^{64}$ big integers, however, the sum $x_0 + x_1$ can produce a $1$-bit carry into limb $m+1$, making the recursive operands $(m+1)$ limbs wide! Moreover, subtracting $z_2$ and $z_0$ requires signed borrow handling. Production big-int libraries (like GMP) either work with redundant limb representations or handle the top carry bit explicitly outside the $m$-limb recursion.