Karatsuba Multiplication, in depth: limb arrays, the subtractive middle term, cutoffs, lopsided operands and testing

Three products from three evaluation points

Think of each operand as a polynomial in a base B. If a number has 2m limbs and you split it at limb m, you get a = a1 B^m + a0, which is the linear polynomial a(t) = a1 t + a0 evaluated at t = B^m. The product of two linear polynomials is a quadratic, c(t) = z2 t² + z1 t + z0, and a quadratic is fixed by its values at three points. Karatsuba picks three cheap points: t = 0 gives z0 = a0 b0, the point at infinity (the leading coefficients) gives z2 = a1 b1, and t = 1 or t = -1 gives the third value. Evaluating at t = 1 needs (a0 + a1)(b0 + b1). Evaluating at t = -1 needs (a0 - a1)(b0 - b1). Interpolation then recovers z1 with a few additions.

This is why Karatsuba is often called Toom-2, or Toom-2.2 in GMP's naming, where the function is mpn_toom22_mul. Toom-3 splits into three pieces, evaluates at five points and makes five recursive products instead of nine, for an exponent of log₃5, about 1.465. Each step up the Toom ladder lowers the exponent but adds evaluation and interpolation work, so higher variants only pay off for larger operands. Drop the carries and the same code multiplies polynomials over a ring.

The subtractive middle term and carry growth

The textbook middle term is z1 = (a0 + a1)(b0 + b1) - z0 - z2. It is correct, but on limb arrays it has a practical problem. Adding two m-limb numbers can produce a carry out of the top limb, so the sum may need m + 1 limbs. Lengths stop being clean halves at every level, and scratch-space accounting gets awkward.

The subtractive form uses the t = -1 evaluation instead: z1 = z0 + z2 + (a0 - a1)(b1 - b0). Expanding the product gives a0 b1 - a0 b0 - a1 b1 + a1 b0, and adding z0 and z2 leaves exactly a0 b1 + a1 b0. A difference of two m-limb numbers never needs more than m limbs. The cost is a sign: compute |a0 - a1| and |b1 - b0| with a comparison first, multiply the magnitudes, and track the product's sign as the product of the two signs. When either difference is zero the product is zero and the middle term is just z0 + z2, which saves a whole recursive call.

Squaring is even simpler. With both operands equal, the middle becomes z0 + z2 - (a0 - a1)², and a square is never negative, so no sign tracking is needed. All three recursive calls are squarings, which are cheaper at the leaves because each cross product is computed once and doubled. That is why libraries keep a separate squaring path with its own threshold.

A limb-level implementation

The implementation below works on little-endian lists of 32-bit limbs. Python ints hold only single limbs and small carries, so every carry and borrow is explicit, as in C. add_shifted recombines in place and accepts a sign, so the middle term never builds a negative intermediate.

BASE_BITS = 32
BASE = 1 << BASE_BITS
MASK = BASE - 1
CUTOFF = 32                      # limbs; measure on your machine


def normalize(a):
    while len(a) > 1 and a[-1] == 0:
        a.pop()
    return a


def compare(a, b):               # both normalized
    if len(a) != len(b):
        return 1 if len(a) > len(b) else -1
    for i in reversed(range(len(a))):
        if a[i] != b[i]:
            return 1 if a[i] > b[i] else -1
    return 0


def sub(a, b):                   # requires a >= b
    out, borrow = [], 0
    for i in range(len(a)):
        d = a[i] - (b[i] if i < len(b) else 0) - borrow
        borrow = 1 if d < 0 else 0
        out.append(d & MASK)     # d + BASE when d is negative
    return normalize(out)


def absdiff(a, b):
    """Return (|a - b|, sign of a - b)."""
    s = compare(a, b)
    return (sub(a, b), s) if s >= 0 else (sub(b, a), -1)


def schoolbook(a, b):
    out = [0] * (len(a) + len(b))
    for i, ai in enumerate(a):
        carry = 0
        for j, bj in enumerate(b):
            t = out[i + j] + ai * bj + carry
            out[i + j] = t & MASK
            carry = t >> BASE_BITS
        out[i + len(b)] += carry
    return normalize(out)


def add_shifted(acc, x, shift, sign=1):
    """acc += sign * x * BASE**shift in place. The final value must be >= 0."""
    carry, i = 0, shift
    for limb in x:
        t = acc[i] + sign * limb + carry
        acc[i] = t & MASK
        carry = t >> BASE_BITS   # floor shift: a borrow arrives as -1
        i += 1
    while carry:
        t = acc[i] + carry
        acc[i] = t & MASK
        carry = t >> BASE_BITS
        i += 1


def karatsuba(a, b):
    if min(len(a), len(b)) <= CUTOFF:
        return schoolbook(a, b)
    m = max(len(a), len(b)) // 2
    a0, a1 = normalize(a[:m]), normalize(a[m:] or [0])
    b0, b1 = normalize(b[:m]), normalize(b[m:] or [0])
    z0 = karatsuba(a0, b0)
    z2 = karatsuba(a1, b1)
    da, sa = absdiff(a0, a1)
    db, sb = absdiff(b1, b0)
    out = [0] * (len(a) + len(b) + 1)
    add_shifted(out, z0, 0)
    add_shifted(out, z2, 2 * m)
    add_shifted(out, z0, m)      # middle = z0 + z2 + sa*sb*|da*db|
    add_shifted(out, z2, m)
    if sa * sb != 0:
        add_shifted(out, karatsuba(da, db), m, sa * sb)
    return normalize(out)

Three details matter. The split point comes from the longer operand, and the halves are normalized because they may have leading zero limbs. The output buffer has a spare limb, because z0 + z2 at offset m can carry past the final length before the subtraction brings it back. And the middle product is skipped when either sign is zero.

Worked example: 1234 x 5678 in base 100

Run one level by hand in base B = 100, so each limb is two decimal digits. Multiply a = 1234 by b = 5678. The split is a1 = 12, a0 = 34, b1 = 56, b0 = 78, with m = 1 limb.

The direct product 1234 x 5678 is 7,006,652, so the recombination is right. Now compare the additive form: (a0 + a1)(b0 + b1) = 46 x 134 = 6164, and 6164 - 2652 - 672 = 2840 as well. Notice that 134 does not fit in one base-100 limb: that is carry growth in a toy example. The middle value 2840 also exceeds the base, so recombination carries 28 into the next limb, which is the job of add_shifted.

Cost, constants and the cutoff

The recurrence is T(n) = 3T(n/2) + cn, where the linear term covers the splits, the two differences and the shifted additions. The recursion tree has 3^k nodes at depth k, each doing work proportional to n/2^k, so the per-level cost grows by a factor of 1.5 at each level and the leaves dominate. The total is Theta(n^{log2 3}), about n^1.585. The Big-O guide explains why constants still matter at real sizes.

Here they are large. Each level does several linear passes with branches and allocations, while schoolbook is a tight multiply-accumulate loop. So every implementation runs schoolbook below a cutoff. Real libraries publish their choices:

CPython's cutoff applies to both operands: it uses the school algorithm unless both contain more than 70 digits, which is about 2,100 bits. Java switches to Karatsuba only when both magnitude arrays exceed 80 ints, roughly 2,560 bits, so RSA-2048 multiplications in Java never reach it. For your own code, measure: the cutoff depends on limb width, allocation cost and the CPU's multiplier.

Lopsided operands and scratch memory

Splitting at half the longer operand works badly when sizes differ a lot. If a has 1,000 limbs and b has 60, the high half of b is empty and z2 is zero. CPython handles this in k_lopsided_mul: cut the long operand into slices the size of the short one, multiply each slice with balanced Karatsuba, and add the partial products at their offsets:

def lopsided(a, b):              # assumes len(a) >= 2 * len(b)
    k = len(b)
    out = [0] * (len(a) + k + 1)
    for start in range(0, len(a), k):
        chunk = normalize(a[start:start + k])
        add_shifted(out, karatsuba(chunk, b), start)
    return normalize(out)

The cost is about len(a) / k balanced k-by-k products. A dispatcher in karatsuba should route here when one length is at least twice the other.

Memory is the other practical issue. The version above allocates fresh lists at every level. C implementations allocate one scratch buffer of about 2n limbs up front and pass slices down the recursion, since each level's needs shrink geometrically. Depth is only log₂(n / cutoff), so allocation and cache behaviour matter far more than stack depth.

Failure modes and how to test for them

Most Karatsuba bugs give correct answers for most inputs, which is what makes them dangerous.

The last row deserves a warning. The sign comparison and the zero shortcut take different paths for different secret values, which is a timing side channel in code that handles keys. Constant-time libraries fold the sign in with masks or skip Karatsuba at key sizes where the gain is small anyway.

Testing is cheap because Python's own multiplication is a trusted oracle:

import random

def to_limbs(x):
    out = []
    while True:
        out.append(x & MASK)
        x >>= BASE_BITS
        if not x:
            return out

def from_limbs(a):
    return sum(limb << (BASE_BITS * i) for i, limb in enumerate(a))

edge = [0, 1, MASK, (1 << 4096) - 1, 1 << 4095]
for _ in range(3000):
    x = random.choice(edge + [random.getrandbits(random.randint(1, 9000))])
    y = random.choice(edge + [random.getrandbits(random.randint(1, 9000))])
    assert from_limbs(karatsuba(to_limbs(x), to_limbs(y))) == x * y

The edge values matter most: all-ones operands force the longest carry chains and powers of two make a difference zero.

Beyond Karatsuba

Karatsuba suits a middle band of sizes, from a few thousand bits until Toom-3 takes over. Past that, FFT-based methods multiply long polynomials pointwise in a transformed domain and carry once. The FFT article covers the transform itself. Schönhage-Strassen, published in 1971, runs in O(n log n log log n). In 2019 Harvey and van der Hoeven gave an O(n log n) algorithm, which is of theoretical interest: its crossover point is far beyond any practical size.

In modular exponentiation every multiplication has the same size, so tuning for that one size beats asymptotics. In fast Fibonacci by doubling, operands grow to millions of digits and the multiplication algorithm dominates run time.

What to do next

• Type in the limb implementation and run the property test, including the edge values, until it passes 3,000 random cases.
• Benchmark schoolbook against Karatsuba for 8 to 512 limbs and set CUTOFF at the crossover you measure.
• Add the lopsided dispatcher and a squaring path that uses (a0 - a1) squared.
• Replace per-call allocations with one preallocated scratch buffer.
• Read CPython's k_mul and k_lopsided_mul in Objects/longobject.c and map each block to the steps on this page.
• If you write cryptographic code, decide explicitly whether data-dependent branches are acceptable, and use a constant-time library if they are not.
• Move on to Toom-3 to see how five evaluation points generalise the same idea.