Умножение длинных чисел методом Карацубы

На днях нужно было разобраться с этим алгоритмом, но беглый поиск в google ничего путнего не дал. На Хабре тоже нашлась только одна статья, которая мне не особо помогла. Разобравшись, попробую поделиться с общественностью в доступной форме:

Алгоритм

Алгоритм Карацубы — метод быстрого умножения со сложностью вычисления n^log₂3. В то время, как наивный алгоритм, умножение в столбик, требует n² операций. Следует заметить, что при длине чисел короче нескольких десятков знаков (точнее определяется экспериментально), быстрее работает обычное умножение.
Представим, что есть два числа A и B длиной n в какой-то системе счисления BASE:
A = a_n-1a_n-2...a₀
B = b_n-1a_n-2...a₀, где a_?, b_? — значение в соотв. разряде числа.
Каждое из них можно представить в виде суммы их двух частей, половинок длиной m = n / 2 (если n нечетное, то одна часть короче другой на один разряд:
A₀ = a_m-1a_m-2...a₀
A₁ = a_n-1a_n-2...a_m
A = A₀ + A₁ * BASE^m

B₀ = b_m-1b_m-2...b₀
B₁ = b_n-1b_n-2...b_m
B = B₀ + B₁ * BASE^m

Тогда: 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}
Здесь нужно 4 операции умножения (части формулы * BASE^{? * m} не являются умножением, фактически указывая место записи результата, разряд). Но с другой стороны:
( A₀ + A₁ ) * ( B₀ + B₁ ) = A₀ * B₀ + A₀ * B₁ + A₁ * B₀ + A₁ * B₁
Посмотрев на выделенные части в обоих формулах. После несложных преобразований количество операций умножения можно свести к 3-м, заменив два умножения на одно и несколько операций сложения и вычитания, время выполнения которых на порядок меньше:
A₀ * B₁ + A₁ * B₀ = ( A₀ + A₁ ) * ( B₀ + B₁ ) — A₀ * B₀ — A₁ * B₁

Окончательный вид выражения:
A * B = A₀ * B₀ + (( A₀ + A₁ ) * ( B₀ + B₁ ) — A₀ * B₀ — A₁ * B₁ ) * BASE^m + A₁ * B₁ * BASE^{2 * m}

Графическое представление:

Пример

Для примера умножим два восьмизначных числа в десятичной системе 12345 и 98765:

Из изображении явно видна рекурсивная природа алгоритма. Для числа длиной менее четырех разрядов применялось наивное умножение.

Реализация на C++

Наверное, следует начать с того, как хранятся длинные числа. Длинные числа удобно представлять в виде массивов, где каждый элемент соответсвует разряду, причем младшие разряды хранятся в элементах с меньшими индексами (то есть задом-наперед) — так их удобнее обрабатывать. Например:
int long_value[] = { 9, 8, 7, 6, 5, 4} //представление числа 456789
Для увеличения быстродействия желательно выбрать за основание системы счисления максимальное число в рамках базовых типов. Но при этом на него накладываются следующие условия:

1. Квадрат максимального числа в выбранной системе счисления должен помещаться в выбранный базовый тип. Это необходимо для хранения произведения одного разряда на другой в промежуточных вычислениях.
2. Выбранный базовый тип желательно брать знаковый. Это позволить избавиться от нескольких промежуточных нормализаций.
3. Лучше, чтобы в разряд помещалось сумма из нескольких квадратов максимального числа. Это позволит избавиться от нескольких промежуточных нормализаций.

Ниже приведена рабочая функция умножения с комментариями со всеми необходимыми вспомогательными объявлениями и функциями. Для более высокой производительности следует поменять основание системы счисления, тип для хранения разрядов и раскоментировать небольшой отрывок кода в месте отвечающем за наивное умножение:

1. #include <cstring>
2.  
3. #define BASE 10 //система счисления
4. #define MIN_LENGTH_FOR_KARATSUBA 4 //числа короче умножаются квадратичным алгоритмом
5. typedef int digit; //взят только для разрядов числа
6. typedef unsigned long int size_length; //тип для длинны числа
7.  
8. using namespace std;
9.  
10. struct long_value { //тип для длинных чисел
11.   digit *values; //массив с цифрами числа записанными в обратном порядке
12.   size_length length; //длинна числа
13. };
14.  
15. long_value sum(long_value a, long_value b) {
16.   /* функция для суммирования двух длинных чисел. Если суммируются числа разной длинны
17.   * то более длинное передется в качестве первого аргумента. Возвращает новое
18.   * ненормализованное число.
19.   */
20.   long_value s;
21.   s.length = a.length + 1;
22.   s.values = new digit[s.length];
23.  
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. }
30.  
31. long_value &sub(long_value &a, long_value b) {
32.   /*функция для вычитания одного длинного числа из другого. Изменяет содержимое первого
33.   * числа. Возвращает ссылку на первое число. Результат не нормализован.
34.   */
35.   for (size_length i = 0; i < b.length; ++i)
36.     a.values[i] -= b.values[i];
37.   return a;
38. }
39.  
40. void normalize(long_value l) {
41.   /*Нормализация числа - приведение каждого разряда в соответствие с системой счисления.
42.   *
43.   */
44.   for (size_length i = 0; i < l.length - 1; ++i) {
45.     if (l.values[i] >= BASE) { //если число больше максимального, то организовавается перенос
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) { //если меньше - заем
50.       digit carryover = (l.values[i] + 1) / BASE - 1;
51.       l.values[i + 1] += carryover;
52.       l.values[i] -= carryover * BASE;
53.     }
54.   }
55. }
56.  
57. long_value karatsuba(long_value a, long_value b) {
58.   long_value product; //результирующее произведение
59.   product.length = a.length + b.length;
60.   product.values = new digit[product.length];
61.  
62.   if (a.length < MIN_LENGTH_FOR_KARATSUBA) { //если число короче то применять наивное умножение
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.         /*В случае изменения MIN_LENGTH_FOR_KARATSUBA или BASE расскоментировать следующие
68.         * строки и подобрать соотв. значения для исключения переполнения разрядов.
69.         * Например для десятичной системы счисления число 100, означает, что организовавается
70.         * перенос 1 через один разряд, 200 - перенос 2 через один разрряд, 5000 - 5 через два.
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 { //умножение методом Карацубы
78.     long_value a_part1; //младшая часть числа a
79.     a_part1.values = a.values;
80.     a_part1.length = (a.length + 1) / 2;
81.  
82.     long_value a_part2; //старшая часть числа a
83.     a_part2.values = a.values + a_part1.length;
84.     a_part2.length = a.length / 2;
85.  
86.     long_value b_part1; //младшая часть числа b
87.     b_part1.values = b.values;
88.     b_part1.length = (b.length + 1) / 2;
89.  
90.     long_value b_part2; //старшая часть числа b
91.     b_part2.values = b.values + b_part1.length;
92.     b_part2.length = b.length / 2;
93.  
94.     long_value sum_of_a_parts = sum(a_part1, a_part2); //cумма частей числа a
95.     normalize(sum_of_a_parts);
96.     long_value sum_of_b_parts = sum(b_part1, b_part2); //cумма частей числа 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.     // произведение сумм частей
100.  
101.     long_value product_of_first_parts = karatsuba(a_part1, b_part1); //младший член
102.     long_value product_of_second_parts = karatsuba(a_part2, b_part2); //старший член
103.     long_value sum_of_middle_terms = sub(sub(product_of_sums_of_parts, product_of_first_parts), product_of_second_parts);
104.     //нахождение суммы средних членов
105.  
106.     /*
107.     * Суммирование многочлена
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];
116.  
117.     /*
118.     * Зачистка
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.   }
126.  
127.   normalize(product); //конечная нормализация числа
128.  
129.   return product;
130. }

* This source code was highlighted with Source Code Highlighter.