The other day it was necessary to deal with this algorithm, but a cursory search in google did not give anything worthwhile. On Habré, too, there was only one article that did not really help me. Having understood, I will try to share with the public in an accessible form:

Algorithm

The Karatsuba algorithm is a fast multiplication method with the complexity of calculating n ^{log ₂ 3} . At that time, as a naive algorithm, multiplication in a column requires n ² operations. It should be noted that when the length of the numbers is shorter than a few dozen characters (or rather determined experimentally), ordinary multiplication works faster.

Imagine that there are two numbers A and B of length n in some number system BASE:

A = a _n-1 a _n-2 ... a ₀

B = b _n-1 a _n-2 ... a ₀ , where a _? b _? - value acc. discharging numbers.

Each of them can be represented as the sum of two parts, halves of length m = n / 2 (if n is odd, then one part is shorter than the other by one digit:

A ₀ = a _m-1 a _m-2 ... a ₀

A ₁ = a _n-1 a _n-2 ... a _m

A = A ₀ + A ₁ * BASE ^m



B ₀ = b _m-1 b _m-2 ... b ₀

B ₁ = b _n-1 b _n-2 ... b _m

B = B ₀ + B ₁ * BASE ^m



Then: A * B = (A ₀ + A ₁ * BASE ^m ) * (B ₀ + B ₁ * BASE ^m ) = A ₀ * B ₀ + A ₀ * B ₁ * BASE ^m + A ₁ * B ₀ * BASE ^m + A ₁ * B ₁ * BASE ^{2 * m} = A ₀ * B ₀ + ( A ₀ * B ₁ + A ₁ * B ₀ ) * BASE ^m + A ₁ * B ₁ * BASE ^{2 * m}

Here you need 4 multiplication operations (parts of the formula * BASE ^{? * M} are not multiplication, in fact, indicating the place where the result is written, the digit). But in other way:

(A ₀ + A ₁ ) * (B ₀ + B ₁ ) = A ₀ * B ₀ + A ₀ * B ₁ + A ₁ * B ₀ + A ₁ * B ₁

Looking at the highlighted parts in both formulas. After simple transformations, the number of multiplication operations can be reduced to the 3rd by replacing two multiplications by one and several addition and subtraction operations, the execution time of which is an order of magnitude less:

A ₀ * B ₁ + A ₁ * B ₀ = (A ₀ + A ₁ ) * (B ₀ + B ₁ ) - A ₀ * B ₀ - A ₁ * B ₁

')

The final look of the expression:

A * B = A ₀ * B ₀ + ((A ₀ + A ₁ ) * (B ₀ + B ₁ ) - A ₀ * B ₀ - A ₁ * B ₁ ) * BASE ^m + A ₁ * B ₁ * BASE ^{2 * m}



Graphic representation:

Example

For example, multiply two eight-digit numbers in the decimal system 12345 and 98765:



The image clearly shows the recursive nature of the algorithm. For a number less than four digits in length, naive multiplication was applied.

C ++ implementation

Probably should start with how long numbers are stored. It is convenient to represent long numbers as arrays, where each element corresponds to a discharge, and the lower digits are stored in elements with smaller indices (that is, backwards), so it is more convenient to process them. For example:

int long_value[] = { 9, 8, 7, 6, 5, 4} // 456789

To increase performance, it is desirable to choose the maximum number within the base types for the base of the number system. But at the same time the following conditions are imposed on it:

1. The square of the maximum number in the selected number system should be placed in the selected base type. It is necessary to store the product of one digit to another in intermediate calculations.
2. The selected base type is desirable to take the sign. This will allow to get rid of several intermediate normalizations.
3. It is better that the sum of several squares of the maximum number be placed in the discharge. This will get rid of several intermediate normalizations.

Below is the working function of multiplication with comments with all the necessary auxiliary declarations and functions. For better performance, you should change the base of the number system, the type for storing the digits, and uncomment a small code snippet at the place responsible for the naive multiplication:

1. #include <cstring>
3. #define BASE 10 // number system
4. #define MIN_LENGTH_FOR_KARATSUBA 4 // numbers are shorter multiplied by a quadratic algorithm
5. typedef int digit; // taken only for digits
6. typedef unsigned long int size_length; // type for long numbers
8. using namespace std;
10. struct long_value { // type for long numbers
11. digit * values; // array with numbers written in reverse order
12. size_length length; // is long numbers
13. };
15. long_value sum (long_value a, long_value b) {
16. / * function to add two long numbers. If numbers of different length are added together.
17. * then the longer is passed as the first argument. Returns new
18. * unnormalized number.
19. * /
20. long_value s;
21. s.length = a.length + 1;
22. s.values ​​= new digit [s.length];
24. s.values ​​[a.length - 1] = a.values ​​[a.length - 1];
25. s.values ​​[a.length] = 0;
26. for (size_length i = 0; i <b.length; ++ i)
27. s.values ​​[i] = a.values ​​[i] + b.values ​​[i];
28. return s;
29. }
31. long_value & sub (long_value & a, long_value b) {
32. / * function to subtract one long number from another. Changes the contents of the first
33. * numbers. Returns a link to the first number. The result is not normalized.
34. * /
35. for (size_length i = 0; i <b.length; ++ i)
36. a.values ​​[i] - = b.values ​​[i];
37. return a;
38. }
40. void normalize (long_value l) {
41. / * Normalization of the number - bringing each digit in accordance with the number system.
42. *
43. * /
44. for (size_length i = 0; i <l.length - 1; ++ i) {
45. if (l.values ​​[i]> = BASE) { // if the number is greater than the maximum, then a transfer is organized
46. digit carryover = l.values ​​[i] / BASE;
47. l.values ​​[i + 1] + = carryover;
48. l.values ​​[i] - = carryover * BASE;
49. } else if (l.values ​​[i] <0) { // if less - loan
50. digit carryover = (l.values ​​[i] + 1) / BASE - 1;
51. l.values ​​[i + 1] + = carryover;
52. l.values ​​[i] - = carryover * BASE;
53. }
54. }
55. }
57. long_value karatsuba (long_value a, long_value b) {
58. long_value product; // resulting product
59. product.length = a.length + b.length;
60. product.values ​​= new digit [product.length];
62. if (a.length <MIN_LENGTH_FOR_KARATSUBA) { // if the number is shorter then apply a naive multiplication
63. memset (product.values, 0, sizeof (digit) * product.length);
64. for (size_length i = 0; i <a.length; ++ i)
65. for (size_length j = 0; j <b.length; ++ j) {
66. product.values ​​[i + j] + = a.values ​​[i] * b.values ​​[j];
67. / * If you change MIN_LENGTH_FOR_KARATSUBA or BASE, uncomment the following
68. * lines and pick up acc. values ​​for avoiding overflow discharges.
69. * For example, for the decimal number system, the number 100 means that it is organized
70. * transfer 1 through one digit, 200 - transfer 2 through one digit, 5000 - 5 through two.
71. * if (product.values ​​[i + j]> = 100) {
72. * product.values ​​[i + j] - = 100;
73. * product.values ​​[i + j + 2] + = 1;
74. *}
75. * /
76. }
77. } else { // multiplication by the Karatsuba method
78. long_value a_part1; // the younger part of a
79. a_part1.values ​​= a.values;
80. a_part1.length = (a.length + 1) / 2;
82. long_value a_part2; // the upper half of a
83. a_part2.values ​​= a.values ​​+ a_part1.length;
84. a_part2.length = a.length / 2;
86. long_value b_part1; // the younger part of the number b
87. b_part1.values ​​= b.values;
88. b_part1.length = (b.length + 1) / 2;
90. long_value b_part2; // the highest part of the number b
91. b_part2.values ​​= b.values ​​+ b_part1.length;
92. b_part2.length = b.length / 2;
94. long_value sum_of_a_parts = sum (a_part1, a_part2); // sum of the parts of a
95. normalize (sum_of_a_parts);
96. long_value sum_of_b_parts = sum (b_part1, b_part2); // sum of parts of number b
97. normalize (sum_of_b_parts);
98. long_value product_of_sums_of_parts = karatsuba (sum_of_a_parts, sum_of_b_parts);
99. // product of parts sums
101. long_value product_of_first_parts = karatsuba (a_part1, b_part1); // junior member
102. long_value product_of_second_parts = karatsuba (a_part2, b_part2); // senior member
103. long_value sum_of_middle_terms = sub (sub (product_of_sums_of_parts, product_of_first_parts), product_of_second_parts);
104. // find the sum of average members
106. / *
107. * Summation of a polynomial
108. * /
109. memcpy (product.values, product_of_first_parts.values,
110. product_of_first_parts.length * sizeof (digit));
111. memcpy (product.values ​​+ product_of_first_parts.length,
112. product_of_second_parts.values, product_of_second_parts.length
113. * sizeof (digit));
114. for (size_length i = 0; i <sum_of_middle_terms.length; ++ i)
115. product.values ​​[a_part1.length + i] + = sum_of_middle_terms.values ​​[i];
117. / *
118. * Stripping
119. * /
120. delete [] sum_of_a_parts.values;
121. delete [] sum_of_b_parts.values;
122. delete [] product_of_sums_of_parts.values;
123. delete [] product_of_first_parts.values;
124. delete [] product_of_second_parts.values;
125. }
127. normalize (product); // final number normalization
129. return product;
130. }

* This source code was highlighted with Source Code Highlighter .