Умножение Карацубы и C++ 11

Хочу в очередной раз затронуть метод реализации умножения Карацубы с использованием возможностей стандарта C++11. Данный алгоритм неоднократно рассматривался здесь («Умножение длинных чисел методом Карацубы», «Алгоритм Карацубы для умножения двух чисел»), но видимо из-за того, что я не умею их готовить, первый вариант не работал с числами разной длины, а второй делает не совсем то, что было нужно.

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

Оглавление



Введение




Всех нас учили умножать в столбик в школе. Это самый простой алгоритм, который известен уже много тысяч лет:



Даже Андрей Николаевич Колмогоров в 1956 году сформулировал гипотезу (которая заключалась в нижней оценке умножения величиной порядка ), так как если бы существовал какой-либо другой более быстрый алгоритм, то за такой огромный промежуток времени он был бы найден.

Псевдокод наивного умножения прост как и сам метод:

multiply(x[0 ... l], y[0 ... r]):
res = [0 ... r+l]
for (i = 0, i < r; ++i):
    carry = 0
    for (j = 0, j < l; ++j):
        res[i + j] += carry + x[i] * y[j]
        carry = res[i + j] / base // base -  база представления числа
        res[i + j] %= base
    res[i + l] += carry

За простоту порой приходится платить производительностью, но этот алгоритм можно оптимизировать и не вычислять остаток на каждом шаге.

Через несколько лет после формулировки гипотезы Колмогорова, Анатолий Алексеевич Карацуба нашел более быстрый метод. Его подход был обобщен до парадигмы «разделяй и властвуй». Чтобы понять как это работает, рассмотрим два числа длины , которые мы разобьем на две части длины :

Теперь заметим, что[1]:


Видно, что необходимо сделать 4 умножения и тогда сложность ничем не отличается от наивного алгоритма. Но Анатолий Алексеевич Карацуба заметил, что обойтись можно 3 умножениями чисел длины , , . Действительно:


Мы обошлись тремя умножениями вместо четырех и следовательно время работы алгоритма Карацубы удовлетворяет соотношению[2]:
,
что в итоге дает общую сложность алгоритма .
Псеводкод алгоритма умножения Карацубы:
Karatsuba_mul(X, Y):
    // X, Y - целые числа длины n
    n = max(размер X, размер Y)
    если n = 1: вернуть X * Y
    X_l = левые n/2 цифр X
    X_r = правые n/2 цифр X
    Y_l = левые n/2 цифр Y
    Y_r = правые n/2 цифр Y
    Prod1 = Karatsuba_mul(X_l, Y_l)
    Prod2 = Karatsuba_mul(X_r, Y_r)
    Prod3 = Karatsuba_mul(X_l + X_r, Y_l + Y_r)
    вернуть Prod1 * 10 ^ n + (Prod3 - Prod1 - Prod2) * 10 ^ (n / 2) + Prod2


И пример на небольших числах, чтобы закрепить механизм работы:

a = 12
b = 81

res = Karatsuba_mul(a, b):
    // размер a = размер b = 2
    n = max( размер a, размер b) // n = 2
    X_l = 1, X_r = 2 // 1 | 2
    Y_l = 8, Y_r = 1 // 8 | 1
    Prod1 = Karatsuba_mul(1, 8) // Prod1 = 8
    Prod2 = Karatsuba_mul(2, 1) // Prod2 = 2
    Prod3 = Karatsuba_mul(3, 9) // Prod3 = 27
    вернуть 8 * 10 ^ 2 + (27 - 2 - 8) * 10 + 2
-----------------------------------------------
res = 972

Реализация


Вот мы и готовы приступить к реализации алгоритма на языке C++. В интернете я находил несколько реализаций, использующие C-стиль написания кода, что несколько затрудняет чтение его для новичков. Поэтому я решил насколько это возможно использовать улучшения доступные в стандарте C++11. Да, это замедлит код, но ведь здесь нас интересует в первую очередь простота для понимания и удобочитаемость.

  1. Хранение числа. Используем стандартный вектор целых чисел, с которым все, изучающие C++, знакомы. Длинное число будем читать в строку и с конца разбивать на разряды, соответствующие выбранной базе (в начале — 10).
    Например, на вход получили число:
    123456789000000000
    

    В нашем контейнере оно будет хранится так:

    |0|1|2|3|4|5|...|n|
     0 0 0 0 0 0 ... 1
    

    Код функции get_number()
    vector<int> get_number(istream& is) {
        string snum;
        vector<int> vnum;
        // индикатор разрядов
        unsigned int dig = 1;
        int n = 0;
        
        is >> snum;
        
        for (auto it = snum.crbegin(); it != snum.crend(); ++it) {
            n += (*it - '0') * dig;
            dig *= dig_size;
            // если разряд равен базе, то выталкиваем число в вектор
            if (dig == base) {
                vnum.push_back(n);
                n = 0;
                dig = 1;
            }
        }
        
        if (n != 0) {
            vnum.push_back(n);
        }
        
        return vnum;
    }
    


  2. Получение числа. На вход у нас могут поступать числа разной длины и нам для успешной работы алгоритма желательно привести к одной и той же длине, кратной 2 (так как мы постоянно разбиваем наши «длинные» числа пополам). Напишем функцию extend_vec(), которая брала бы наш вектор и удлиняла его как-то так:

    first = {4}; // 4; size = 1
    second = {3, 2, 1} // 123; size = 3
    auto n = max(first.size(), second.size());
    
    extend_vec(first, n); // добавить 3 нуля
    extend_vec(second, n); // добавить 1 ноль
    

    Код функции extend_vec()
    void extend_vec(vector<int>& v, size_t len) {    
        while (len & (len - 1)) {
            ++len;
        }
            
        v.resize(len);
    }
    


  3. Умножение. Здесь стоит поговорить о нескольких оптимизациях, которые стоит сделать. Мы не будем считать остатки и переносить их в старшие разряда на каждом рекурсивном вызове, а сделаем это в конце. И для перемножения двух чисел с длинной меньше, скажем, 128 будем использовать наивный алгоритм, так как он является меньшей константой, чем алгоритм Карацубы.

    Код функции naive_mul()
    vector<int> naive_mul(const vector<int>& x, const vector<int>& y) {
        auto len = x.size();
        vector<int> res(2 * len);
            
        for (auto i = 0; i < len; ++i) {
            for (auto j = 0; j < len; ++j) {
                res[i + j] += x[i] * y[j];
            }
        }
            
        return res;
    }
    


    Код функции karatsuba_mul()
    vector<int> karatsuba_mul(const vector<int>& x, const vector<int>& y) {
        auto len = x.size();    
        vector<int> res(2 * len);
        
        if (len <= len_f_naive) {
            return naive_mul(x, y);
        }
        
        auto k = len / 2;
        
        vector<int> Xr {x.begin(), x.begin() + k};
        vector<int> Xl {x.begin() + k, x.end()};
        vector<int> Yr {y.begin(), y.begin() + k};
        vector<int> Yl {y.begin() + k, y.end()};
        
        vector<int> P1 = karatsuba_mul(Xl, Yl);
        vector<int> P2 = karatsuba_mul(Xr, Yr);    
            
        vector<int> Xlr(k);
        vector<int> Ylr(k);
        
        for (int i = 0; i < k; ++i) {
            Xlr[i] = Xl[i] + Xr[i];
            Ylr[i] = Yl[i] + Yr[i];
        }
        
        vector<int> P3 = karatsuba_mul(Xlr, Ylr);
        
        for (auto i = 0; i < len; ++i) {
            P3[i] -= P2[i] + P1[i];
        }
        
        for (auto i = 0; i < len; ++i) {
            res[i] = P2[i];
        }
    
        for (auto i = len; i < 2 * len; ++i) {
            res[i] = P1[i - len];
        }
    
        for (auto i = k; i < len + k; ++i) {
            res[i] += P3[i - k];
        }
        
        return res;
    }
    


  4. Нормализация. Осталось сделать все переносы и можно выводить результат (или использовать для дальнейших вычислений).

    Код функции finalize()
    void finalize(vector<int>& res) {
        for (auto i = 0; i < res.size(); ++i) {
            res[i + 1] += res[i] / base;
            res[i] %= base;
        }
    }
    


    И выводим результат, дополняя нулями при использование базы, большей, чем 10.

    Код функции print_res()
    void print_res(const vector<int>& v, ostream& os) {
        auto it = v.crbegin();
        
        // Passing leading zeroes
        while (!*it) {
            ++it;
        }
        
        while (it != v.crend()) {
            int z = -1;
            int num = *it;
            
            if (num == 0) {
                num += 1;
            }
            
            if (num < add_zero) {
                z = 1;         
                
                while ((num *= dig_size) < add_zero) {
                    ++z;
                }
            }
            
            if (z > 0) {
                while (z--) {
                    os << '0';
                }
            }
            os << *it++;
        }
        
        os << endl;
    }
    



Сравнение скорости работы наивного алгоритма и алгоритма Карацубы


Для сборки тестовой программы использовался Clang++ с ключом -O3. Результаты тестирования для представления чисел с базой 10 приведены на рисунке 1.

Время расчета произведения (база 10)

Рисунок 1. Время расчета произведения двух чисел, используя представление с базой 10

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

Время расчета произведения (база 100)

Рисунок 2. Время расчета произведения двух чисел, используя представление с базой 100

Выводы


Вот и все, мы разобрали с вам этот простой и эффективный способ умножения. Надеюсь, это материал будет полезен и многие новички, которые только начинают изучение алгоритмов не будут больше впадать в ступор (ну не зашел он у меня с первого раза в свое время).

Ещё есть куда оптимизировать данную реализацию:
  • увеличить базу в которой хранятся числа в векторе. Сейчас нормализация числа делается в самом конце, что вызывает переполнение стандартных типов в C++. Возможно стоит хранить числа в массиве/векторе типа unsigned long long и вычислять остатки с переносами на каждом этапе умножения. Либо использовать «длинное» представление остатка.
  • отказаться от векторов в пользу массивов и не использовать выделение левой и правой части числа с помощью итераторов.

На этом все, всем спасибо за внимание.
Исходный проект, который использовался при написании статьи, находится здесь.

Список используемой литературы


  1. С. Дасгупта Алгоритмы: Перевод с английского А. С. Куликова под редакцией А. Шеня [Текст] / С. Дасгупта, Х. Пападимитриу, У. Вазирани. -Москва: МЦНМО, 2014 — 320 с.
  2. Karatsuba algorithm [Электронный ресурс] / Wikipedia — URL: en.wikipedia.org/wiki/Karatsuba_algorithm
  3. А. С. Куликов Алгоритмы и структуры данных [Электронный ресурс] / А. С. Куликов — URL: https://stepic.org/course/Алгоритмы-и-структуры-данных-63/syllabus
Поделиться публикацией
Комментарии 12
    +1
    Помимо оптимизации постоянных копирований памяти:

    — в extend_vec можно использовать метод vector::resize
    — нужно оптимизировать заглавную картинку (сейчас она размером 1920 на 1920)
      0
      И ещё: сейчас функции предполагают, что аргументы имею одинаковые и правильные длины. Если это не так, они молча возвращаюи ерунду. Это опасно. Лучше либо сделать автоматическое приведение длин в начале функции, либо хотя бы проверку на корректность длин.
      +2
      Раз уж вы написали статью о больших числах, то делать, например, вот так:

      int n = max(first.size(), second.size());

      как-то некорректно. У вас будет переполнение при длине чисел более чем в 2^31 (на x86 и x86-64 в современных компиляторах).

      Правильно:

      size_t n = max(first.size(), second.size()

      Или, если вы уж упомянули C++11, то:

      auto n = max(first.size(), second.size());

      Ну и во многих других местах тоже.
        0
        Да, спасибо за замечание, исправил.
      • НЛО прилетело и опубликовало эту надпись здесь
          +1
          1. Статья содержит реализацию, значит нужно было вместо основания выбрать не 10 или 100, а 2^k.
          И к выбирается обычно таким, чтобы 2^k умещалось в точности в машинное слово.

          2. Нужно было привести другие формулы разложения чисел для перемножения для A=A_1 + 2^k A_2, B=B_1 + 2^k B_2
          Ваше
          A B= A_1 B_1 + 2^k ((A_1 + A_2)(B_1 + B_2) — (A_1 B_1 + A_2 B_2)) + 2^(2k) A_2 B_2

          Используется в gmplib.org
          A B= A_1 B_1 — 2^k ((A_1 — A_2)(B_1 — B_2) — (A_1 B_1 + A_2 B_2)) + 2^(2k) A_2 B_2

          Еще умножение можно заменить возведением в квадрат с той же сложностью.

          3. Метод Карацубы начинает выигрывать у обычного умножения на современных компьютерах при размере чисел равных 8 машинным словам. Но это зависит от реализации и конкретного компа.

          4. Если разбивать число сразу на r — частей (метод умножения Toom–Cook), то сложность можно довести до O(n^(1+eps)).

          5. Самый быстрый метод умножения Шёнхаге — Штрассена основанный на быстром преобразование Фурье имеет сложность O(n log(n) log(log(n))). Фактически его асимптотическая сложность совпадает со сложением и вычитанием.

          6. Деление можно заменить умножением на обратный (его для точной арифметики получают алгоритмом основанным на методе Ньютона для приближенного решения нелинейных уравнений) и тогда асимптотическая сложность деления равна сложности умножения.
            0
            Шёнхаге — Штрассена продержался в лидерах до 2007-го. Алгоритм Фюрера быстрее. Есть еще один, не помню его названия, аналогичный АФ по быстродействию.
              0
              1. Идея такой реализации, которая приведена здесь, чтобы все переносы делались в конце. И нам необходимо, чтобы все промежуточные результаты умещались в стандартный тип, который может хранить числа до (причем тип должен быть знаковым, так как в алгоритме Карацубы используется вычитание. Если использовать систему , умещая в машинное слово, то где Вы будете хранить результат произведения? Поправьте меня, если я не прав.
              2. Можно добавить, но идея от этого не изменится
              3. Безусловно. На моей машине до определенной длины числа оба алгоритма работали мгновенно и только после перехода через определенный порог наивный алгоритм стал существенно замедляться.
              4. Можно, но алгоритм Тома-Кука и делает больше промежуточной работы по разбиению числа на r–частей, поэтому и существенный выигрыш будет при ещё бóльших длинах чисел.
                +1
                Когда-то, примерно в 1998 году у меня была реализация длинной арифметики на C++. По скорости она в 2-3 раза делала на тот момент GMP. Сейчас я пользуюсь GMP и не заморачиваюсь с ассемблером под конкретный процессор.

                И так, идея реализации с минимальным использованием ассемблера. Цифра везде машинное слово, а типы C++ зависят от конкретной архитектуры. Кстати для организации всяких сдвигов нужно знать где расположен значащий бит в начале или в конце машинного слова

                struct LongInt {
                int mAlloc; \\ размер выделенной памяти
                int mSize; \\ abs(mSize) текущий размер числа и sign(mSize) его знак
                unsigned int* mLimb; \\ массив цифр
                };

                unsgned int a, b;
                если a+b < a, то произошло переполнение и единичку должны перенести, a+b при этом правильный результат для младшей цифры
                если a-b > a, то произошло переполнение и единичку должны занять, a-b при этом правильный результат для младшей цифры

                Ассемблер нужен для двузначного умножения и деления. Поскольку с 60-х архитектура большинства процессоров поддерживала создание длинной арифметики (вспомните хотя бы lisp)

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

                при делении задается две цифры делимого и цифра делителя, в результате цифра частного и цифра остатка

                Как выглядит все на разных ассемблерах можно посмотреть в gmplib.org в исходниках в папке mpn. В ней mpn\generic без ассемблера.

              0
              Через FT умножение описывалось на хабре?
              А, вижу, выше упомянули. FFT алгоритм хорош.
                0
                Да, вот здесь.
                  0
                  Да, спасибо. Лет 10 назад, когда впервые встала потребность перемножения чисел по нескольку миллионов знаков, эта идея поразила своей красотой.

              Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

              Самое читаемое