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

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

Алгоритм


Алгоритм Карацубы — метод быстрого умножения со сложностью вычисления nlog23. В то время, как наивный алгоритм, умножение в столбик, требует n2 операций. Следует заметить, что при длине чисел короче нескольких десятков знаков (точнее определяется экспериментально), быстрее работает обычное умножение.
Представим, что есть два числа A и B длиной n в какой-то системе счисления BASE:
A = an-1an-2...a0
B = bn-1an-2...a0, где a?, b? — значение в соотв. разряде числа.
Каждое из них можно представить в виде суммы их двух частей, половинок длиной m = n / 2 (если n нечетное, то одна часть короче другой на один разряд:
A0 = am-1am-2...a0
A1 = an-1an-2...am
A = A0 + A1 * BASEm

B0 = bm-1bm-2...b0
B1 = bn-1bn-2...bm
B = B0 + B1 * BASEm

Тогда: A * B = ( A0 + A1 * BASEm ) * ( B0 + B1 * BASEm ) = A0 * B0 + A0 * B1 * BASEm + A1 * B0 * BASEm + A1 * B1 * BASE2 * m = A0 * B0 + ( A0 * B1 + A1 * B0 ) * BASEm + A1 * B1 * BASE2 * m
Здесь нужно 4 операции умножения (части формулы * BASE? * m не являются умножением, фактически указывая место записи результата, разряд). Но с другой стороны:
( A0 + A1 ) * ( B0 + B1 ) = A0 * B0 + A0 * B1 + A1 * B0 + A1 * B1
Посмотрев на выделенные части в обоих формулах. После несложных преобразований количество операций умножения можно свести к 3-м, заменив два умножения на одно и несколько операций сложения и вычитания, время выполнения которых на порядок меньше:
A0 * B1 + A1 * B0 = ( A0 + A1 ) * ( B0 + B1 ) — A0 * B0A1 * B1

Окончательный вид выражения:
A * B = A0 * B0 + (( A0 + A1 ) * ( B0 + B1 ) — A0 * B0 — A1 * B1 ) * BASEm + A1 * B1 * BASE2 * m

Графическое представление:
схема умножения


Пример


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

Реализация на 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.
Проголосовать:
+92
Сохранить: