AlgoScope

Karatsuba and Strassen

algorithmadvancedTime O(n^1.585) for Karatsuba, O(n^2.807) for StrassenSpace O(n) per recursion level

Both are the same idea: a multiplication costs far more than an addition once the operands are big, so spend a few extra additions to save one multiplication, and let the saving compound through the recursion. Karatsuba splits two numbers into halves and computes the cross term from one product of sums instead of two products, three half-size multiplications instead of four, which is n to the 1.585 instead of n squared. Strassen does it for 2 x 2 block matrices with seven block products instead of eight, n to the 2.807 instead of n cubed. Neither pays off on small inputs, which is why libraries switch to them only above a threshold.

highlowxyz2 = a cz0 = b d(a+b)(c+d)z1x y12345678

Multiply 1234 by 5678 with Karatsuba. Split each at 2 digits from the right: 1234 = 12 x 10^2 + 34 and 5678 = 56 x 10^2 + 78. Schoolbook expansion needs the four products a c, a d, b c and b d; Karatsuba gets away with three.

Check your understanding

The player pauses before the one decision in this run and asks what happens next. Here it is, with the answer.

  1. (12 + 34)(56 + 78) = 6164. What is the middle term a d + b c?

    • 2840
    • 6164
    • 3324

    Answer: 2840. Subtract z2 = 672 and z0 = 2652 from 6164: the cross terms are what is left.

How it runs, step by step

  1. Multiply 1234 by 5678 with Karatsuba. Split each at 2 digits from the right: 1234 = 12 x 10^2 + 34 and 5678 = 56 x 10^2 + 78. Schoolbook expansion needs the four products a c, a d, b c and b d; Karatsuba gets away with three.

    Karatsuba multiplication of 1234 and 5678.

  2. Two products are unavoidable: the highs, z2 = 12 x 56 = 672, and the lows, z0 = 34 x 78 = 2652. They land at 10^4 and at 10^0 in the answer.

    z2 is 672, z0 is 2652.

  3. The middle term needs a d + b c. Instead of two products, multiply the sums once: (12 + 34) x (56 + 78) = 46 x 134 = 6164, which expands to a c + a d + b c + b d. Subtract the two products already known, 672 and 2652, and what remains is exactly a d + b c: z1 = 2840. Three multiplications, not four.

    z1 is 6164 minus 672 minus 2652, which is 2840.

  4. 1234 x 5678 = 672 x 10^4 + 2840 x 10^2 + 2652 = 7006652, with 3 multiplications of half-size numbers instead of 4. Recursing, 3 products at each of log2 n levels gives n^(log2 3), about n^1.585, against n^2 for schoolbook. The extra additions make it slower for small numbers; libraries switch to it somewhere around a few hundred digits.

    1234 times 5678 is 7006652 with three multiplications.

Remember

  • Karatsuba: z1 = (a + b)(c + d) - ac - bd gives the cross term from one product, three instead of four.
  • Strassen: seven products of sums and differences of the blocks, then every entry of C is a sum of them.
  • Three of four per level is n^log2(3); seven of eight is n^log2(7). Only worth it for large inputs.

Where this is used

CryptographyRSA and Diffie-Hellman in OpenSSL

A 2048-bit modular exponentiation is a few thousand multiplications of 2048-bit numbers, so the cost of one big multiply is the cost of the handshake. OpenSSL's bn_mul_recursive in crypto/bn/bn_mul.c is Karatsuba over the limb array of a BIGNUM, in the subtractive form: it splits both operands in half and gets the cross term from the single product (a0 - a1)(b1 - b0), tracking the sign of each difference, which keeps both factors inside the half width instead of carrying into an extra limb. Below BN_MUL_RECURSIVE_SIZE_NORMAL, which is 16 limbs, it falls back to bn_mul_normal, and exact small sizes go to the unrolled bn_mul_comba4 and bn_mul_comba8 kernels, because at a handful of limbs the extra additions and the recursive call cost more than the multiplication they save.

Language runtimesArbitrary precision integers in CPython and Java

Both runtimes carry two multiplication paths and a measured switch point between them. CPython's k_mul in Objects/longobject.c takes over once the shorter operand exceeds KARATSUBA_CUTOFF = 70 internal digits, 30 bits each on a 64-bit build, so roughly 600 decimal digits, and it doubles that cutoff for squaring. java.math.BigInteger switches to Karatsuba once both magnitudes exceed 80 ints and to Toom-Cook 3 once one exceeds 240. The threshold exists because schoolbook multiplication is a tight pass over contiguous memory, and the asymptotically better method only wins once the saved multiplications outweigh its allocations and extra passes.

Numerical librariesGMP under computer algebra systems and GCC

GMP does not have one multiplication algorithm but a ladder - basecase, Karatsuba, Toom-3, Toom-4, Toom-6.5, Toom-8.5, then FFT - with each rung taking over above its own threshold, MUL_TOOM22_THRESHOLD through MUL_FFT_THRESHOLD. Those thresholds are compile-time constants taken from a gmp-mparam.h already tuned for the detected CPU; they are not measured on your machine, though GMP ships a tune/tuneup program you can run by hand and paste the results back in. Sage and PARI/GP rest on GMP, and so does GCC, which requires GMP and MPFR so it can fold constant arithmetic exactly at compile time instead of on the target's floating point.

Scientific computingStrassen in dense linear algebra

Tuned BLAS libraries multiply matrices with the cubic algorithm because it streams cache-friendly blocks, and LAPACK is built on that kernel. The BLIS-based implementation by Huang, Smith, Henry and van de Geijn shows the usual objections are softer than they look: by folding the sums and differences of blocks into the packing step a high-performance GEMM already performs, one-level Strassen needs no workspace beyond the buffers GEMM allocates anyway, and it beats conventional GEMM well below the sizes Strassen is normally said to require. It stays an opt-in path rather than the default because the error bound is weaker and odd dimensions must be handled at every level. The same seven-product scheme is what DeepMind's AlphaTensor searched for, treating it as a rank decomposition of the 2 x 2 multiplication tensor and turning up new ones for larger blocks in modular arithmetic.

Why it works this way

Why a few extra additions change the exponent at all

The additions live in the combine step, which is linear, while the multiplications decide how many subproblems the recursion spawns. Karatsuba is T(n) = 3T(n/2) + O(n), and the master theorem gives n^log2(3) because the 3^d subproblems at depth d outgrow the linear work at every level; four subproblems would give n^2 instead. So the win is not 25 percent off the work, it is a change in how fast the work multiplies as you descend.

Why not split into three or four pieces?

You can, and that is Toom-Cook: cut each number into k parts, evaluate at 2k - 1 points, multiply there, then interpolate back, for an exponent of log_k(2k - 1). That is 1.585 for k = 2, 1.465 for k = 3, 1.404 for k = 4, and it approaches 1 as k grows. The catch is the interpolation, whose additions, shifts and small divisions grow quickly with k, so each k only wins over a window of operand sizes. That is why a serious bignum library does not pick one k, it switches between them as the operands grow.

The halves are not quite halves: a + b can carry

b and d fit in m digits, but a + b and c + d can need m + 1, so the middle recursive call runs on slightly wider operands than the other two and implementations allocate an extra limb for it rather than assuming the split size holds. The split point needs the same care: x and y must be cut at the same m, and if one operand is far shorter than the other, splitting by the longer one's length can leave a = 0 so the recursion stops shrinking. CPython keeps a separate lopsided path for exactly that case, when one operand has at least twice the digits of the other.

What Strassen costs you: accuracy and memory

The seven products are built from sums and differences of blocks, so intermediate values can be far larger than the entries of the answer and then cancel. Strassen satisfies only a norm-wise error bound rather than the entry-by-entry one the cubic algorithm gives, and the bound degrades with each level of recursion, which is why LAPACK and the standard BLAS keep the classical kernel. A direct implementation also holds temporaries for the seven products, and each level wants even dimensions, so odd sizes have to be padded or peeled off before recursing.

Read more

Next up