Karatsuba multiplication

Karatsuba multiplication computes the product of two $n$ -digit numbers with about $n^{l o g_{2} 3} \approx n^{1.585}$ single-digit multiplications, fewer than the $n^{2}$ of the schoolbook method^[1]^[2]. It splits each factor into two halves, $a = a_{1} B^{m} + a_{0}$ and $b = b_{1} B^{m} + b_{0}$ , and observes that the three products $a_{1} b_{1}$ , $a_{0} b_{0}$ and $(a_{0} + a_{1}) (b_{0} + b_{1})$ determine all four products of halves that the full product needs, since the middle term $a_{1} b_{0} + a_{0} b_{1}$ is the third minus the other two. Applying the idea recursively turns the four half-size multiplications of the obvious recursion into three. Anatoly Karatsuba discovered the method in 1960 as a student in Kolmogorov's seminar, refuting Kolmogorov's conjecture that multiplication requires order $n^{2}$ operations^[3]. It was the first subquadratic multiplication algorithm and opened the study of fast arithmetic.

§ 01

The algorithm

With $z_{2} = a_{1} b_{1}$ , $z_{0} = a_{0} b_{0}$ and $z_{1} = (a_{0} + a_{1}) (b_{0} + b_{1}) - z_{2} - z_{0} = a_{1} b_{0} + a_{0} b_{1}$ ,

ab = z_{2} B^{2 m} + z_{1} B^{m} + z_{0}

The three products have half the length, up to a carry digit from the additions, and are computed recursively; the additions, subtractions and shifts cost $O (n)$ . The running time therefore satisfies $T (n) = 3 T (n /2) + O (n)$ , whose solution is $T (n) = O (n^{l o g_{2} 3})$ ^[2].

Splitting both factors in half. The high product z₂ and low product z₀ are placed at the ends of the result, and the middle term z₁ is recovered from one more product of the sums of halves.

Schoolbook and Karatsuba multiplication of natural numbers in base 10000, counting digit operations.

(* Multiplying natural numbers stored as arrays of base-10000 digits,
   least significant first, counting single-digit multiplications and
   the digits passed through additions and subtractions. *)
let base = 10_000
let mults = ref 0 and adds = ref 0
 
let normalize a = let n = ref (Array.length a) in
  while !n > 0 && a.(!n - 1) = 0 do decr n done; Array.sub a 0 !n
 
let carry a = (* propagate carries; digits may be negative on entry *)
  let a = Array.append a [| 0; 0 |] in
  for i = 0 to Array.length a - 2 do
    let q = if a.(i) >= 0 then a.(i) / base else - ((- a.(i) + base - 1) / base) in
    a.(i) <- a.(i) - q * base; a.(i + 1) <- a.(i + 1) + q
  done;
  normalize a
 
let add a b = adds := !adds + max (Array.length a) (Array.length b);
  carry (Array.init (max (Array.length a) (Array.length b)) (fun i ->
    (if i < Array.length a then a.(i) else 0) + (if i < Array.length b then b.(i) else 0)))
let sub a b = adds := !adds + max (Array.length a) (Array.length b);
  carry (Array.init (max (Array.length a) (Array.length b)) (fun i ->
    (if i < Array.length a then a.(i) else 0) - (if i < Array.length b then b.(i) else 0)))
let shift a k = if a = [||] then a else Array.append (Array.make k 0) a
 
(* Schoolbook: every digit of a times every digit of b. *)
let schoolbook a b =
  if a = [||] || b = [||] then [||]
  else begin
    let r = Array.make (Array.length a + Array.length b) 0 in
    Array.iteri (fun i x -> Array.iteri (fun j y -> incr mults; r.(i + j) <- r.(i + j) + x * y) b) a;
    carry r
  end
 
(* Karatsuba: split at half the length, a = a1 B^m + a0, b = b1 B^m + b0,
   and use three products instead of four:
   a b = z2 B^2m + ((a0 + a1)(b0 + b1) - z2 - z0) B^m + z0. *)
let rec karatsuba ?(threshold = 1) a b =
  let n = max (Array.length a) (Array.length b) in
  if n <= threshold || Array.length a = 0 || Array.length b = 0 then schoolbook a b
  else begin
    let m = n / 2 in
    let lo x = normalize (Array.sub x 0 (min m (Array.length x)))
    and hi x = if Array.length x > m then Array.sub x m (Array.length x - m) else [||] in
    let a0 = lo a and a1 = hi a and b0 = lo b and b1 = hi b in
    let z0 = karatsuba ~threshold a0 b0 and z2 = karatsuba ~threshold a1 b1 in
    let z1 = sub (sub (karatsuba ~threshold (add a0 a1) (add b0 b1)) z2) z0 in
    add (add (shift z2 (2 * m)) (shift z1 m)) z0
  end
 
let of_string s = (* decimal string to base-10000 digits *)
  let n = String.length s in
  Array.init ((n + 3) / 4) (fun i -> let hi = n - 4 * i in let lo = max 0 (hi - 4) in int_of_string (String.sub s lo (hi - lo)))
  |> normalize
let to_string a = if a = [||] then "0" else
    let parts = Array.to_list (Array.mapi (fun i d -> if i = Array.length a - 1 then string_of_int d else Printf.sprintf "%04d" d) a) in
    String.concat "" (List.rev parts)

Checking against Zarith.

open Kara

Running it.

85397342226735670654635508695465744950348885357650690581833014653136349617576232159238592126769012450
schoolbook 169 digit products, Karatsuba 86, both equal to Zarith: true
300 random products up to 300 digits agree with Zarith: true

The same identity multiplies polynomials, with the halves taken by degree; there are no carries, and the method is used in computer algebra and in cryptographic implementations of polynomial and binary-field arithmetic.

§ 02

Cost in practice

Digit multiplications and additions for operands of increasing length: schoolbook, pure Karatsuba, and Karatsuba that switches to schoolbook below 16 digits.

open Kara
 
(* Digit multiplications for n-digit operands (n base-10000 digits). *)

Running it.

digits   schoolbook    Karatsuba (mul + add)   cutoff 16 (mul + add)
    16          256        133 + 1392             256 + 0         
    64         4096       1208 + 13535           2235 + 1356      
   256        65536      11099 + 128224         19719 + 18141     
  1024      1048576     100400 + 1176312       176079 + 185871    
  4096     16777216     905724 + 10674789     1581589 + 1756051

Digit operations on a logarithmic scale. Pure Karatsuba saves multiplications but performs many more additions; switching to the schoolbook method for small operands keeps most of the savings and removes most of the additions.

Recursing all the way down saves multiplications but spends them on additions, and at small sizes the schoolbook method is faster. Real implementations therefore recurse only above a threshold, tuned to the machine; in GMP the Karatsuba range begins at a few dozen machine words^[7]. On 4,096-digit operands the version with a cutoff does about a fifth of the schoolbook work in total.

§ 03

Generalizations

Karatsuba's method evaluates the linear polynomials $a_{1} t + a_{0}$ and $b_{1} t + b_{0}$ at $t = 0, 1, \infty$ , multiplies the values, and interpolates the quadratic product. Toom observed in 1963 that splitting into $k$ parts and evaluating at $2 k - 1$ points reduces $k^{2}$ products to $2 k - 1$ , giving exponents $lo g_{k} (2 k - 1)$ that approach 1; Cook made the construction practical in his 1966 thesis^[4]^[9]. Toom-3, with three parts and five products, costs $O (n^{1.465})$ . Evaluating at roots of unity instead, by the fast Fourier transform, gives Schönhage and Strassen's $O (n lo g n lo g lo g n)$ algorithm of 1971^[5], and Harvey and van der Hoeven reached $O (n lo g n)$ in 2019^[6]. Montgomery found Karatsuba-like formulae for five, six and seven parts that use fewer products than Toom–Cook for the same number of parts^[8]. Libraries combine all of these, each above its own threshold^[7].

§ 04

History

In 1960 Andrey Kolmogorov conjectured in his seminar at Moscow State University that multiplying two $n$ -digit numbers requires on the order of $n^{2}$ elementary operations. Within a week, Anatoly Karatsuba, then 23, found the three-product method and disproved it^[3]. The result was published in 1962 under the names of Karatsuba and Yuri Ofman, although Ofman's contribution was a separate result^[1]^[3]. Toom's generalization followed in 1963^[4], Cook's thesis in 1966^[9], and FFT-based multiplication in 1971^[5].