On multiplying great numbers by partition

Enuntiatio

Propositio V. Theorema. The product of two whole numbers, each of n digits, may be formed with no more than three multiplications of numbers of half the magnitude, together with a number of additions, subtractions, and shiftings proportional to n.

Let the multiplicands be split at the middle digit, so that x = a·B + b and y = c·B + d, where B = 10^m is the splitting base. The schoolbook expansion demands the four products ac, ad, bc, bd. I say that ad + bc is obtained from a single further product (a+b)(c+d) by subtracting the two products ac and bd already in hand; whence three multiplications suffice where four were thought necessary. Iterated upon the halves, this partition obeys the recurrence T(n) = 3·T(n/2) + O(n), and therefore consumes Θ(n^{log₂3}) ≈ Θ(n^{1.585}) digit-operations — for all sufficiently great n, strictly fewer than the Θ(n²) of multiplication by the schools.

Expressio

def karatsuba(x: int, y: int) -> int:
    """Multiply two integers by Karatsuba's three-product partition.

    Returns x * y exactly, using recursively three half-length
    multiplications in place of the schoolbook four.

    >>> karatsuba(1234, 5678)
    7006652
    >>> karatsuba(1234, 5678) == 1234 * 5678
    True
    >>> karatsuba(0, 999999)
    0
    >>> karatsuba(-31, 41) == -31 * 41
    True
    >>> all(karatsuba(i, j) == i * j for i in range(60) for j in range(60))
    True
    """
    # Reduce signed operands to the non-negative case.
    if x < 0 or y < 0:
        sign = -1 if (x < 0) ^ (y < 0) else 1
        return sign * karatsuba(abs(x), abs(y))

    # Base case: a single-digit factor is multiplied directly.
    if x < 10 or y < 10:
        return x * y

    # Split each operand at its middle digit: value = high * B + low.
    n = max(len(str(x)), len(str(y)))
    m = n // 2
    B = 10 ** m
    a, b = divmod(x, B)   # x = a*B + b,  0 <= b < B
    c, d = divmod(y, B)   # y = c*B + d,  0 <= d < B

    # Three recursive products, not four.
    ac = karatsuba(a, c)
    bd = karatsuba(b, d)
    ad_plus_bc = karatsuba(a + b, c + d) - ac - bd  # the Karatsuba identity

    # Recompose: x*y = ac*B^2 + (ad+bc)*B + bd.
    return ac * (10 ** (2 * m)) + ad_plus_bc * B + bd


if __name__ == "__main__":
    import doctest
    doctest.testmod()
    # Concrete witnesses of correctness.
    assert karatsuba(99, 99) == 9801
    assert karatsuba(31415926, 27182818) == 31415926 * 27182818
    assert karatsuba(2 ** 200, 3 ** 130) == (2 ** 200) * (3 ** 130)
    assert karatsuba(-123456789, 987654321) == -123456789 * 987654321
    print("Karatsuba verified.")

Demonstratio

The demonstration falls into two parts: first the correctness of the identity, then the cost of its iteration.

Part I — The identity is exact.

Let B = 10^m and write each n-digit operand in the radix-B positional form

x = a·B + b, y = c·B + d,

where a, c are the high parts and b, d the low parts, with 0 ≤ b, d < B. This is precisely what divmod(x, B) produces, so the split loses nothing. By the distributive law,

x·y = (aB + b)(cB + d) = ac·B² + (ad + bc)·B + bd. (∗)

The naive reading of (∗) requires the four products ac, ad, bc, bd. Now observe the algebraic step at the heart of the method:

(a + b)(c + d) = ac + ad + bc + bd,

and therefore

ad + bc = (a + b)(c + d) − ac − bd. (†)

The right side of (†) reuses ac and bd, which are already computed for (∗); the only new product is (a + b)(c + d). Substituting (†) into (∗) gives the recomposition coded above. Thus exactly three multiplications — ac, bd, and (a+b)(c+d) — produce all of x·y, the remaining operations being additions, subtractions, and multiplications by powers of ten (mere digit-shifts).

That the recursion terminates and is everywhere exact follows by strong induction on max(len(x), len(y)). Base: if either factor is below ten, the product is returned directly and is exact. Step: assume karatsuba multiplies exactly any pair of operands shorter than n digits. Because m = n//2 ≥ 1 whenever n ≥ 2, each of a, b, c, d is strictly shorter than n digits; and a + b, c + d each carry at most one digit beyond the half-length, hence remain strictly shorter than n digits as well. The three recursive calls therefore act on strictly smaller inputs and, by hypothesis, return exact products; equation (∗) combined with (†) then yields x·y exactly. The sign reduction at the top handles negatives by the law of signs and does not disturb the magnitude computation. Hence karatsuba(x, y) = x·y for all integers x, y. (This is what the inline assertions and the exhaustive small-input doctest confirm on concrete numbers.)

Part II — The cost is Θ(n^{log₂3}).

Let T(n) count the digit-operations to multiply two operands of at most n digits. The non-recursive work of one call — forming a, b, c, d by a split; the sums a+b and c+d; two subtractions; and the final shifts and additions — each touches O(n) digits, so it is bounded by c·n for some constant c. The three recursive calls each act on operands of at most ⌈n/2⌉ + 1 digits, which is n/2 + O(1); absorbing the additive constant into the analysis gives the recurrence

T(n) = 3·T(n/2) + O(n), T(1) = O(1). (‡)

Apply the Master Theorem with a = 3, b = 2, f(n) = Θ(n). The watershed exponent is

log_b a = log₂ 3 ≈ 1.58496.

Since f(n) = Θ(n^1) = O(n^{log₂3 − ε}) for any ε with 0 < ε < log₂3 − 1 (such ε exists because log₂3 − 1 ≈ 0.585 > 0), we are in Case 1 of the Master Theorem, the leaves of the recursion tree dominating the cost. Therefore

T(n) = Θ(n^{log₂3}) = Θ(n^{1.585}).

For the sceptic who distrusts the theorem, expand (‡) directly. At depth k there are 3^k subproblems, each of size n/2^k, contributing 3^k · c·(n/2^k) = c·n·(3/2)^k digit-operations. Summing over the k = 0,…,log₂n levels,

T(n) ≤ c·n · Σ_{k=0}^{log₂ n} (3/2)^k = c·n · ( (3/2)^{log₂n + 1} − 1 ) / ( 3/2 − 1 ) = O( n · (3/2)^{log₂ n} ).

But (3/2)^{log₂ n} = n^{log₂(3/2)} = n^{log₂3 − 1}, whence n · (3/2)^{log₂n} = n^{log₂3}, recovering T(n) = O(n^{log₂3}); the leaf-level alone contributes 3^{log₂n} = n^{log₂3} subproblems, giving the matching lower bound Ω(n^{log₂3}).

Finally, the comparison with the schoolbook method. School multiplication forms all n² digit-products, costing Θ(n²). Since log₂3 < 2, we have n^{log₂3} = o(n²); the ratio n²/n^{log₂3} = n^{2 − log₂3} = n^{0.415} → ∞. Hence for all sufficiently great n the partition method strictly dominates, as asserted.

Q.E.D.