Karatsuba multiplication

also: divide and conquer multiplication

Multiplying two n-limb integers with three half-size multiplications instead of four, by reusing (a0+a1)(b0+b1) to obtain the middle term. That turns the schoolbook n^2 into about n^1.585, and it only pays above a threshold where the extra additions cost less than the multiplication saved.

Schoolbook multiplication of two $n$ -limb integers costs $n^{2}$ limb products. Split each in half, $a = a_{1} B + a_{0}$ and $b = b_{1} B + b_{0}$ , and the product needs four of them:

ab = a_{1} b_{1} B^{2} + (a_{1} b_{0} + a_{0} b_{1}) B + a_{0} b_{0}

The middle coefficient is recoverable from the other two plus one extra product, which turns four into three:

a_{1} b_{0} + a_{0} b_{1} = (a_{1} + a_{0}) (b_{1} + b_{0}) - a_{1} b_{1} - a_{0} b_{0}

Three half-size multiplications and a handful of additions give the recurrence $T (n) = 3 T (n /2) + O (n)$ , so

T (n) = Θ (n^{l o g_{2} 3}) \approx Θ (n^{1.585})

The recursion, with the threshold that makes it worth doing.

let karatsuba_threshold = 32
 
let rec mag_mul a b =
  if Array.length a < karatsuba_threshold
     || Array.length b < karatsuba_threshold then mag_mul_school a b
  else begin
    let k = (max (Array.length a) (Array.length b) + 1) / 2 in
    let a0, a1 = split a k and b0, b1 = split b k in
    let z0 = mag_mul a0 b0 in
    let z2 = mag_mul a1 b1 in
    let z1 = mag_sub (mag_sub (mag_mul (mag_add a1 a0) (mag_add b1 b0)) z2) z0 in
    mag_add (mag_add (mag_shift_limbs z2 (2 * k)) (mag_shift_limbs z1 k)) z0
  end

The threshold is not a detail. Below it the extra additions and allocations cost more than the multiplication saved, so a Karatsuba implementation without a schoolbook base case is slower than schoolbook on everything anyone actually multiplies. Above a few thousand limbs it is in turn superseded by Toom-Cook and then by FFT-based methods.

referenced by

Knuth's Algorithm D

Building a CAS in OCaml Part 1Arbitrary-precision integers from scratch: base-2^30 limbs, schoolbook and Karatsuba multiplication, Knuth algorithm D division, and decimal conversion.