Karatsuba algorithm: Difference between revisions

Revision as of 00:27, 12 August 2006

Latest revision as of 20:43, 4 May 2025

(421 intermediate revisions by more than 100 users not shown)

Line 1:

{{short description|Algorithm for integer multiplication}}

The '''Karatsuba''' [[multiplication algorithm]], a technique for quickly multiplying large numbers, was discovered by Karatsuba and Ofman in 1962. Its [[time complexity]] is [[Big O notation|Θ]]<math>(n^{\log_23})</math>. This makes it faster than the classical [[Big O notation|Θ]](''n''²) algorithm, although for small numbers it is not as fast.

{{Infobox algorithm

|class = [[Multiplication algorithm]]

|image = <!-- filename only, no "File:" or "Image:" prefix, and no enclosing [[brackets]] -->

|caption =

|data =

|time = <!-- Worst time big-O notation -->

|best-time =

|average-time =

|space = <!-- Worst-case space complexity; auxiliary space

(excluding input) if not specified -->

}}

[[File:Karatsuba_multiplication.svg|thumb|300px|Karatsuba multiplication of az+b and cz+d (boxed), and 1234 and 567 with z=100. Magenta arrows denote multiplication, amber denotes addition, silver denotes subtraction and cyan denotes left shift. (A), (B) and (C) show recursion with z=10 to obtain intermediate values.]]

The '''Karatsuba algorithm''' is a fast [[multiplication algorithm]] for [[Integer|integers]]. It was discovered by [[Anatoly Karatsuba]] in 1960 and published in 1962.<ref name="kara1962">

{{cite journal

| author = A. Karatsuba and Yu. Ofman

| title = Multiplication of Many-Digital Numbers by Automatic Computers

| journal = Proceedings of the USSR Academy of Sciences

| volume = 145

| year = 1962

| pages = 293–294

| url = https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=dan&paperid=26729&option_lang=eng

| postscript = . Translation in the academic journal ''[[Physics-Doklady]]'', '''7''' (1963), pp. 595–596}}

</ref><ref name="kara1995">

{{cite journal

| author = A. A. Karatsuba

| title = The Complexity of Computations

| journal = Proceedings of the Steklov Institute of Mathematics

| volume = 211

| pages =169–183

| year = 1995

| url = http://www.ccas.ru/personal/karatsuba/divcen.pdf

| postscript = . Translation from Trudy Mat. Inst. Steklova, 211, 186–202 (1995)}}

</ref><ref name="knuthV2">

Knuth D.E. (1969) ''[[The Art of Computer Programming]]. v.2.'' Addison-Wesley Publ.Co., 724 pp.

</ref> It is a [[divide-and-conquer algorithm]] that reduces the multiplication of two ''n''-digit numbers to three multiplications of ''n''/2-digit numbers and, by repeating this reduction, to at most <math> n^{\log_23}\approx n^{1.58}</math> single-digit multiplications. It is therefore [[Asymptotic complexity|asymptotically faster]] than the [[long multiplication|traditional]] algorithm, which performs <math>n^2</math> single-digit products.

The Karatsuba algorithm was the first multiplication algorithm asymptotically faster than the quadratic "grade school" algorithm.

The [[Toom–Cook multiplication|Toom–Cook algorithm]] (1963) is a faster generalization of Karatsuba's method, and the [[Schönhage–Strassen algorithm]] (1971) is even faster, for sufficiently large ''n''.

==ExampleHistory==

The standard procedure for multiplication of two ''n''-digit numbers requires a number of elementary operations proportional to <math>n^2\,\!</math>, or <math>O(n^2)\,\!</math> in [[big-O notation]]. [[Andrey Kolmogorov]] conjectured that the traditional algorithm was ''[[asymptotically optimal]],'' meaning that any algorithm for that task would require <math>\Omega(n^2)\,\!</math> elementary operations.

To multiply two ''n''-digit numbers ''x'' and ''y'' represented in base ''W'', where ''n'' is even and equal to 2''m'' (if ''n'' is odd instead, or the numbers are not of the same length, this can be corrected by adding zeros at the left end of ''x'' and/or ''y''), we can write

In 1960, Kolmogorov organized a seminar on mathematical problems in [[cybernetics]] at the [[Moscow State University]], where he stated the <math>\Omega(n^2)\,\!</math> conjecture and other problems in the [[Computational complexity theory|complexity of computation]]. Within a week, Karatsuba, then a 23-year-old student, found an algorithm that multiplies two ''n''-digit numbers in <math>O(n^{\log_2 3})</math> elementary steps, thus disproving the conjecture. Kolmogorov was very excited about the discovery; he communicated it at the next meeting of the seminar, which was then terminated. Kolmogorov gave some lectures on the Karatsuba result at conferences all over the world (see, for example, "Proceedings of the International Congress of Mathematicians 1962", pp. 351–356, and also "6 Lectures delivered at the International Congress of Mathematicians in Stockholm, 1962") and published the method in 1962, in the [[Proceedings of the USSR Academy of Sciences]]. The article had been written by Kolmogorov and contained two results on multiplication, Karatsuba's algorithm and a separate result by [[Yuri Petrovich Ofman|Yuri Ofman]]; it listed "A. Karatsuba and Yu. Ofman" as the authors. Karatsuba only became aware of the paper when he received the reprints from the publisher.<ref name="kara1995"/>

: ''x'' = ''x''₁ ''W''^''m'' + ''x''₂

: ''y'' = ''y''₁ ''W''^''m'' + ''y''₂

==Algorithm==

with ''m''-digit numbers ''x''₁, ''x''₂, ''y''₁ and ''y''₂. The product is given by

===Basic step===

: ''xy'' = ''x''₁''y''₁ ''W''^2''m'' + (''x''₁''y''₂ + ''x''₂''y''₁) ''W''^''m'' + ''x''₂''y''₂

The basic principle of Karatsuba's algorithm is [[Divide-and-conquer algorithm|divide-and-conquer]], using a formula that allows one to compute the product of two large numbers <math>x</math> and <math>y</math> using three multiplications of smaller numbers, each with about half as many digits as <math>x</math> or <math>y</math>, plus some additions and digit shifts. This basic step is, in fact, a generalization of [[Multiplication algorithm#Complex number multiplication|a similar complex multiplication algorithm]], where the [[imaginary unit]] {{mvar|i}} is replaced by a power of the [[radix|base]].

Let <math>x</math> and <math>y</math> be represented as <math>n</math>-digit strings in some base <math>B</math>. For any positive integer <math>m</math> less than <math>n</math>, one can write the two given numbers as

so we need to quickly determine the numbers ''x''₁''y''₁, ''x''₁''y''₂ + ''x''₂''y''₁ and ''x''₂''y''₂. The heart of Karatsuba's method lies in the observation that this can be done with only three rather than four multiplications:

:<math>x = x_1 B^m + x_0,</math>

# compute ''x''₁''y''₁, call the result ''A''

:<math>y = y_1 B^m + y_0,</math>

# compute ''x''₂''y''₂, call the result ''B''

# compute (''x''₁ + ''x''₂)(''y''₁ + ''y''₂), call the result ''C''

# compute ''C'' - ''A'' - ''B''; this number is equal to ''x''₁''y''₂ + ''x''₂''y''₁.

where <math>x_0</math> and <math>y_0</math> are less than <math>B^m</math>. The product is then

To compute these three products of ''m''-digit numbers, we can employ the same trick again, effectively using [[recursion]]. Once the numbers are computed, we need to add them together, which takes about ''n'' operations.

<math>

==Complexity==

\begin{align}

If ''T''(''n'') denotes the time it takes to multiply two ''n''-digit numbers with Karatsuba's method, then we can write

xy &= (x_1 B^m + x_0)(y_1 B^m + y_0) \\

&= x_1 y_1 B^{2m} + (x_1 y_0 + x_0 y_1) B^m + x_0 y_0 \\

&= z_2 B^{2m} + z_1 B^m + z_0, \\

\end{align}

</math>

where

:''T''(''n'') = 3 ''T''(''n''/2) + ''cn'' + ''d''

:<math>z_2 = x_1 y_1,</math>

for some constants ''c'' and ''d'', and this [[recurrence relation]] can be solved, giving a time complexity of Θ(''n''^ln(3)/ln(2)). The number ln(3)/ln(2) is approximately 1.585, so this method is significantly faster than long multiplication. Because of the overhead of recursion, Karatsuba's multiplication is slower than long multiplication for small values of ''n''; typical implementations therefore switch to long multiplication if ''n'' is below some threshold.

:<math>z_1 = x_1 y_0 + x_0 y_1,</math>

:<math>z_0 = x_0 y_0.</math>

These formulae require four multiplications and were known to [[Charles Babbage]].<ref>Charles Babbage, Chapter VIII – Of the Analytical Engine, Larger Numbers Treated, [https://archive.org/details/bub_gb_Fa1JAAAAMAAJ/page/n142 <!-- pg=125 --> Passages from the Life of a Philosopher], Longman Green, London, 1864; page 125.</ref> Karatsuba observed that <math>xy</math> can be computed in only three multiplications, at the cost of a few extra additions. With <math>z_0</math> and <math>z_2</math> as before and <math>z_3=(x_1 + x_0) (y_1 + y_0),</math> one can observe that

Implementations differ greatly on what the crossover point is when Karatsuba becomes faster than the classical algorithm, but generally numbers over 2³²⁰ are faster with Karatsuba. [http://www.swox.com/gmp/manual/Karatsuba-Multiplication.html][http://ozark.hendrix.edu/~burch/proj/karat/comment1.html] Karatsuba is surpassed by the asymptotically faster [[Schönhage-Strassen algorithm]] around 2^33,000.

:<math>

\begin{align}

z_1 &= x_1 y_0 + x_0 y_1 \\

&= (x_1 + x_0) (y_0 + y_1) - x_1 y_1 - x_0 y_0 \\

&= z_3 - z_2 - z_0. \\

\end{align}

</math>

Thus only three multiplications are required for computing <math>z_0, z_1</math> and <math>z_2.</math>

==An Objective Caml implementation==

===Example===

(* Decomposition of the numbers into arrays *)

To compute the product of 12345 and 6789, where ''B'' = 10, choose ''m'' = 3. We use ''m'' right shifts for decomposing the input operands using the resulting base (''B''^''m'' = ''1000''), as:

let decomposition10 a=

: 12345 = '''12''' · ''1000'' + '''345'''

let rec log10 a=

: 6789 = '''6''' · ''1000'' + '''789'''

if a<10 then 1 else 1+(log10 (a/10))

in

Only three multiplications, which operate on smaller integers, are used to compute three partial results:

let lga=log10 a in

: ''z''₂ = '''12''' '''×''' '''6''' = 72

let resultat=Array.create (lga) 0 in

: ''z''₀ = '''345''' '''×''' '''789''' = 272205

let rec decomp_aux a i=

: ''z''₁ = ('''12''' + '''345''') '''×''' ('''6''' + '''789''') − ''z''₂ − ''z''₀ = 357 '''×''' 795 − 72 − 272205 = 283815 − 72 − 272205 = 11538

if a=0 then () else

begin

We get the result by just adding these three partial results, shifted accordingly (and then taking carries into account by decomposing these three inputs in base ''1000'' as for the input operands):

let b=a/10 in

: result = ''z''₂ · (''B''^''m'')^''2'' + ''z''₁ · (''B''^''m'')^''1'' + ''z''₀ · (''B''^''m'')^''0'', i.e.

resultat.(i)<-(a-10*b);

: result = 72 · ''1000''² + 11538 · ''1000'' + 272205 = '''83810205'''.

decomp_aux b (i+1)

end

Note that the intermediate third multiplication operates on an input ___domain which is less than two times larger than for the two first multiplications, its output ___domain is less than four times larger, and base-''1000'' carries computed from the first two multiplications must be taken into account when computing these two subtractions.

in

decomp_aux a 0;

===Recursive application===

resultat;;

If ''n'' is four or more, the three multiplications in Karatsuba's basic step involve operands with fewer than ''n'' digits. Therefore, those products can be computed by [[recursion|recursive]] calls of the Karatsuba algorithm. The recursion can be applied until the numbers are so small that they can (or must) be computed directly.

(* Power function, defined recursively *)

In a computer with a full 32-bit by 32-bit [[Binary multiplier|multiplier]], for example, one could choose ''B'' = 2³¹ and store each digit as a separate 32-bit binary word. Then the sums ''x''₁ + ''x''₀ and ''y''₁ + ''y''₀ will not need an extra binary word for storing the carry-over digit (as in [[carry-save adder]]), and the Karatsuba recursion can be applied until the numbers to multiply are only one digit long.

let rec ( ** ) nb pb=

if pb=0 then 1 else nb*(nb**(pb-1));;

===[[Time complexity]] analysis===

Karatsuba's basic step works for any base ''B'' and any ''m'', but the recursive algorithm is most efficient when ''m'' is equal to ''n''/2, rounded up. In particular, if ''n'' is 2^''k'', for some integer ''k'', and the recursion stops only when ''n'' is 1, then the number of single-digit multiplications is 3^''k'', which is ''n''^''c'' where ''c'' = log₂3.

(* Gives both arrays the same length *)

let equilibre vecta vectb=

Since one can extend any inputs with zero digits until their length is a power of two, it follows that the number of elementary multiplications, for any ''n'', is at most <math>3^{ \lceil\log_2 n \rceil} \leq 3 n^{\log_2 3}\,\!</math>.

let pluslong=max (Array.length vecta) (Array.length vectb) in

let m=Array.make_matrix 2 pluslong 0 in

Since the additions, subtractions, and digit shifts (multiplications by powers of ''B'') in Karatsuba's basic step take time proportional to ''n'', their cost becomes negligible as ''n'' increases. More precisely, if ''T''(''n'') denotes the total number of elementary operations that the algorithm performs when multiplying two ''n''-digit numbers, then

for i=0 to pluslong-1 do

begin

:<math>T(n) = 3 T(\lceil n/2\rceil) + cn + d</math>

try m.(0).(i)<-vecta.(i) with _->m.(0).(i)<-0;

end;

for some constants ''c'' and ''d''. For this [[recurrence relation]], the [[master theorem (analysis of algorithms)|master theorem for divide-and-conquer recurrences]] gives the [[big O notation|asymptotic]] bound <math>T(n) = \Theta(n^{\log_2 3})\,\!</math>.

begin

try m.(1).(i)<-vectb.(i) with _->m.(1).(i)<-0;

It follows that, for sufficiently large ''n'', Karatsuba's algorithm will perform fewer shifts and single-digit additions than longhand multiplication, even though its basic step uses more additions and shifts than the straightforward formula. For small values of ''n'', however, the extra shift and add operations may make it run slower than the longhand method.

end

done;

==Implementation==

m.(0),m.(1);;

Here is the pseudocode for this algorithm, using numbers represented in base ten. For the binary representation of integers, it suffices to replace everywhere 10 by 2.<ref>{{cite book |last= Weiss |first= Mark A. |date= 2005 |title= Data Structures and Algorithm Analysis in C++ |publisher= Addison-Wesley|page= 480|isbn= 0321375319}}</ref>

(* Karatsuba multiplication function *)

let karatsuba a b=

The second argument of the split_at function specifies the number of digits to extract from the ''right'': for example, split_at("12345", 3) will extract the 3 final digits, giving: high="12", low="345".

let decompa0=decomposition10 a

and decompb0=decomposition10 b in

<syntaxhighlight lang="c">

let decompa,decompb=equilibre decompa0 decompb0 in

function karatsuba(num1, num2)

let rec karatsuba_aux xi xj yi yj=

if xi=xj(num1 &&< yi=yj10 thenor decompa.(xinum2 < 10)*decompb.(yi) else

return num1 × num2 /* fall back to traditional multiplication */

begin

let x1y1=karatsuba_aux xi ((xi+xj)/2+1) yi ((yi+yj)/2+1)

/* Calculates the size of the numbers. */

and x2y2=karatsuba_aux ((xi+xj)/2) xj ((yi+yj)/2) yj

m = max(size_base10(num1), size_base10(num2))

and x1y2=karatsuba_aux xi ((xi+xj)/2+1) ((yi+yj)/2) yj

m2 = floor(m / 2)

and x2y1=karatsuba_aux ((xi+xj)/2) xj yi ((yi+yj)/2+1)

/* m2 = ceil (m / 2) will also work */

and m2=(xi-xj)/2+1 in

x1y1*(10**(2*m2))+(x1y2+x2y1)*(10**m2)+x2y2

/* Split the digit sequences in the middle. */

end

high1, low1 = split_at(num1, m2)

in

high2, low2 = split_at(num2, m2)

karatsuba_aux (Array.length decompa -1) 0 (Array.length decompa-1) 0;;

/* 3 recursive calls made to numbers approximately half the size. */

z0 = karatsuba(low1, low2)

z1 = karatsuba(low1 + high1, low2 + high2)

z2 = karatsuba(high1, high2)

return (z2 × 10 ^ (m2 × 2)) + ((z1 - z2 - z0) × 10 ^ m2) + z0

</syntaxhighlight>

An issue that occurs when implementation is that the above computation of <math>(x_1 + x_0)</math> and <math>(y_1 + y_0)</math> for <math>z_1</math> may result in overflow (will produce a result in the range <math>B^m \leq \text{result} < 2 B^m</math>), which require a multiplier having one extra bit. This can be avoided by noting that

:<math>z_1 = (x_0 - x_1)(y_1 - y_0) + z_2 + z_0.</math>

This computation of <math>(x_0 - x_1)</math> and <math>(y_1 - y_0)</math> will produce a result in the range of <math>-B^m < \text{result} < B^m</math>. This method may produce negative numbers, which require one extra bit to encode signedness, and would still require one extra bit for the multiplier. However, one way to avoid this is to record the sign and then use the absolute value of <math>(x_0 - x_1)</math> and <math>(y_1 - y_0)</math> to perform an unsigned multiplication, after which the result may be negated when both signs originally differed. Another advantage is that even though <math>(x_0 - x_1)(y_1 - y_0)</math> may be negative, the final computation of <math>z_1</math> only involves additions.

==References==

{{Reflist}}

* A. Karatsuba and Yu Ofman, "Multiplication of Many-Digital Numbers by Automatic Computers." ''Doklady Akad. Nauk SSSR'' Vol. 145 (1962), pp. 293&ndash;294. Translation in ''Physics-Doklady'' 7 (1963), pp. 595&ndash;596.

*[http://mathworld.wolfram.com/KaratsubaMultiplication.html Karatsuba Multiplication on MathWorld]

==External links==

*Bernstein, D. J., "[http://cr.yp.to/papers/m3.pdf Multidigit multiplication for mathematicians]". Covers Karatsuba and many other multiplication algorithms.

* [http://www.cs.pitt.edu/~kirk/cs1501/animations/Karatsuba.html Karatsuba's Algorithm for Polynomial Multiplication]

*{{MathWorld|urlname=KaratsubaMultiplication|title=Karatsuba Multiplication}}

* Bernstein, D. J., "[http://cr.yp.to/papers/m3.pdf Multidigit multiplication for mathematicians]". Covers Karatsuba and many other multiplication algorithms.

{{Number-theoretic algorithms}}

[[Category:Computer arithmetic algorithms]]

[[Category:Multiplication]]

[[Category:Articles with example pseudocode]]

[[Category:Divide-and-conquer algorithms]]