Алгоритм Карацубы для умножения двух чисел

    Как-то раз пришлось реализовывать умножение длинных чисел, через половинки их представления на C++. 128-битные числа были представлены как пара 64-битных. Оказалось что перемножить два 128-битных числа и получить все 256 бит результата по сложности сравнимо с 4-мя произведениями 64-битных половинок. Как же это работает…



    Для умножения использовался метод отечественного математика Анатолия Алексеевича Карацубы. Для начала приведу свой комментарий к функции. Он очень много проясняет по алгоритму.

    //
    // Karatsuba multiplication algorithm
    //
    //            +------+------+
    //            |  x1  |  x0  |
    //            +------+------+
    //                   *
    //            +------+------+
    //            |  y1  |  y0  |
    //            +------+------+
    //                   =
    //     +-------------+-------------+
    //  +  |    x1*y1    |    x0*y0    |
    //     +----+-+------+------+------+
    //          . .             .
    //          . .             .
    //          +-+------+------+
    //       +  | x0*y1 + x1*y0 |
    //          +-+------+------+
    //
    

    В принципе по комментарию становится понятно как вычисляется результат.
    Суть метода можно понять если вывести следующую формулу:



    T символизирует отступ, соответственно T^2 — двойной отступ.
    Вот код, выполняющий эти операции:

    //
    // Params:
    //     T  - is type of x0, x1, y0 and y1 halves
    //     T2 - is type of x, y and half of res
    //
    template<typename T, typename T2>
    inline void Karatsuba_multiply(T * const x, T * const y, T2 * res)
    {
        // Define vars (depending from byte order)
    
        #define ptrRes ((T*)res)
        T2 &  lowWord = (T2&)(ptrRes[0]);
        T2 &  midWord = (T2&)(ptrRes[1]);
        T2 & highWord = (T2&)(ptrRes[2]);
        T  & highByte = (T &)(ptrRes[3]);
        #undef ptrRes
    
        const T & x0 = x[0];
        const T & x1 = x[1];
        const T & y0 = y[0];
        const T & y1 = y[1];
    
        // Multiply action
    
        lowWord  = x0 * y0;
        highWord = x1 * y1;
    
        T2 m1 = x0 * y1;
        T2 m2 = x1 * y0;
    
        midWord += m1;
        if (midWord < m1) highByte++;
        midWord += m2;
        if (midWord < m2) highByte++;
    }
    

    Тип T является типом половинки, типом x0, x1, y0, y1.
    Тип N2 является типом половины результата, в 2 раза больше типа T.

    Пример использования функции:

    int main()
    {
        typedef unsigned char   u8;
        typedef unsigned short u16;
        typedef unsigned int   u32;
    
        u16 a = 1000;
        u16 b = 2000;
        u32 r = 0;
    
        u8  * a_ptr = (u8*)&a;
        u8  * b_ptr = (u8*)&b;
        u16 * r_ptr = (u16*)(void*)&r;
    
        Karatsuba_multiply(a_ptr, b_ptr, r_ptr);
    
        cout << r;
    }
    

    Код с результатом выполнения можно посмотреть тут: codepad.org/1OTeGqhA
    Код с алгоритмом для разных порядков байтов (LE + BE) тут: codepad.org/f5Pxtiq1

    UPDATE1:
    Пользователь mark_ablov заметил, что в коде нехватает #undef.
    Исправленный код: codepad.org/13U4fuTp
    Исправленный полный код (LE + BE): codepad.org/kBazqo8f

    UPDATE2:
    Пользователь ttim заметил, что приведённый алгоритм, не совсем является методом Карацубы и указал на возможность использования 3-х умножений вместо четырёх.

    Формула, вносящая ясность:



    Таким образом, потребуется вычислить лишь 3 произведения:
    1. a*x
    2. b*y
    3. (Ta+b)*(Tx+y)

    К сожалению, мне не удастся воспользоваться этой формулой, т.к. я не имею возможности перемножать числа разрядности 128 бит. Собственно моей задачей и являлось реализовать перемножение 128-битных чисел. Дело в том что числа (Ta+b) и (Tx+y) разрядности 128 бит...

    UPDATE3:
    Пользователь susl продолжил обсуждение алгоритма и показал что вовсе не нужно перемножать 128-битные числа.

    Ещё одна формула:



    Функцию можно переписать следующим образом:

    //
    // Params:
    //     T  - is type of x0, x1, y0 and y1 halves
    //     T2 - is type of x, y and half of res
    //
    template<typename T, typename T2>
    inline void Karatsuba_multiply(T * const x, T * const y, T2 * res)
    {
        // Define vars (depending from byte order)
    
        #define ptrRes ((T*)res)
        T2 &  lowWord = (T2&)(ptrRes[0]);
        T2 &  midWord = (T2&)(ptrRes[1]);
        T2 & highWord = (T2&)(ptrRes[2]);
        T  & highByte = (T &)(ptrRes[3]);
        #undef ptrRes
    
        const T & x0 = x[0];
        const T & x1 = x[1];
        const T & y0 = y[0];
        const T & y1 = y[1];
    
        // Multiply action
    
        lowWord  = x0 * y0;
        highWord = x1 * y1;
    
        //T2 m1 = x0 * y1;
        //T2 m2 = x1 * y0;
        T2 m = (x0+x1)*(y0+y1) - (lowWord + highWord);
    
        //midWord += m1;
        //if (midWord < m1) highByte++;
        //midWord += m2;
        //if (midWord < m2) highByte++;
        midWord += m;
        if (midWord < m) highByte++;
    }
    

    Исправленный пример: http://codepad.org/w0INBD77
    Исправленный пример для LE и BE: http://codepad.org/nB9HqWt1
    Поделиться публикацией
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама
    Комментарии 53
    • +7
      Как то не привычно когда исходный код в статье в виде картинок, по привычке попытался выделить кусок кода, когда изучал его и не получилось. А чем собственно плох тот же Source Code Highlighter?
      • –6
        Ну скажем так, с картинкой проще. Она лежит на официальном хостинге хабра — значит никуда не пропадёт. Форматирование и подсветка в коде именно такие как «я задумал» — и так будет на любом устройстве. Картинки оптимизированы по объёму (все порядка 5-15Кб).
        • +1
          Видимо, Вы решили это продвигать в массы)
          Линк
          • +6
            Но их нельзя скопировать, на моём Andoid Wildfire они — нечитабельные. При этом есть отличный тег
            <source lang="cpp">
                <!-- Исходник -->
            </source>

            Код в этом теге подсветится, будет в стиле Хабра, везде легко читаться и легко копироваться.
            • +1
              С утра переделаю, спасибо за терпение.
            • +2
              Ну текст то тоже на официальном хостинге)
          • +22
            Если что, то Карацуба — это не суровый японский самурай, как может показаться, а известный советский математик.
            Звали его Анатолий Алексеевич Карацуба.
            • 0
              Я думаю это стоит добавить к статье. Спасибо, плюсую)
              • 0
                > Если что, то Карацуба — это не суровый японский самурай, как может показаться

                А почему самурай? Фамилия же не похожа на японскую.
                Я почему-то при прочтении этой фамилии представил колоритного хохла лысого и с чубом.
                • +2
                  А я, почему-то, прочитал фамилию как Курацаба, и подумал что японец)
                  • +5
                    Суровый хохол からつば :)
                    • +6
                      По-моему, хорошая такая японская фамилия. Так как половина, если не больше, японских фамилий — географического происхождения (некие аналоги «имярек с высокой горы», «имярек из далекой деревни, что в топких северных болотах»), «цуба» (場) — одно из наиболее часто встречающихся окончаний, означающих «место». «Кара-цуба» — в зависимости от прочтения, это будет либо «ореховое место», либо «пустое место», либо «китайское место» или «сложное, напряженное место» :)

                      Впрочем, это все фантазии, ни в одном справочнике японских имен-фамилий «Карацуба» нет :) Есть «Камацуба» и «Кавацуба».
                      • +1
                        Ну, в такие дебри я не лез. Просто каной записывается, поэтому вполне себе японская фамилия :)
                    • +1
                      Почему-то пограничника Карацупа все знают и никто не сомневается, а вот одной буквой отличается — и всё…
                  • +2
                    Видимо, Вы решили это продвигать в массы)
                    Линк
                    • +2
                      Не в ту ветку
                      • +1
                        Прямо как тут :)
                        • +1
                          Ох, вот это промазал так промазал =) есть к чему стремиться
                      • 0
                        Следующую статью обязательно оформлю текстом. Просто я не могу смириться, когда код выглядит непотребно. Надеюсь будет не очень жутко.
                        • +3
                          У меня есть статьи с огромным количеством кода. Как Вам кажется, смотрится ли он непотребно?
                          • 0
                            Вот снимок экрана с айфона: http://img830.imageshack.us/i/imagexqh.jpg/. Скажете это нормальные отступы? Или у вас так и есть в коде? Не могу сейчас с компа посмотреть.
                            • 0
                              Что-то imageshack меня обманул с качеством скриншота. Но на скрине видно выравнивание кода…
                              • +1
                                Попробуйте сейчас.
                                • 0
                                  Отступы вы я погляжу уменьшили и всё-равно есть значительная разница между компом и телефоном(было, стало).
                                • +1
                                  Честно говоря, вообще не понимаю, как исходники можно читать с телефона, будь они написаны текстом или в скриншоты запихнуты…
                                  • 0
                                    Их можно даже писать с телефона)) Пробовал на codepad.org, долго но проверить идейку на плюсах можно по дороге домой.
                          • +3
                            код не смотрел внимательно, но бросилось в глаза то, что нет #undef'a после того как ваш утилитарный макрос стал не нужен.
                            • 0
                              Спасибо, сделаем)
                            • +5
                              Я возможно не прав, но по-моему это не алгоритм Карацубы. Приведенный вами алгоритм не имеет никакой ценности казалось бы — для умножения двух чисел длины N он использует 4 умножения чисел длины N/2. То есть по сути ваш алгоритм работает за N^2. Алгоритм перемножения «влоб» тоже будет иметь такую ассимптотику.

                              Алгоритм Карацубы имеет такой же первый шаг, но использует следующую оптимизацию: рассмотрим (a+b)*(x+y) = (ay+bx) + ax + by.
                              Вычислим 3 (!) произведения чисел длины N/2 — (a+b)*(x+y), ax, by. Вычтем из первого второе и третье и получим ay+bx. Подсчитаем T^2(ax)+T(ay+bx)+by.

                              Как-то так.
                              • 0
                                Похоже вы правы. Сейчас перепишу статью. Спасибо!
                                • 0
                                  Посмотрел ещё раз. Да Вы правы, но и я в некотором смысле прав)) Я не могу умножать 128-битные числа, только 64-хбитные. Допишу в статье о возможности сокращения объёма вычислений. Ещё раз Спасибо!
                                  • +2
                                    Так и не надо умножать 128-битные.
                                    a = ah * 2^64 + al
                                    b = bh * 2^64 + bl
                                    тогда
                                    a * b = (ah * 2^64 + al) * (bh * 2^64 + bl) = 2^128 * ah * bh + 2^64 * (ah * bl + bh * al) + al * bl
                                    но ah * bl + bh * al = (ah + al) * (bh + bl) — (ah *bh + al * bl)

                                    поэтому вычисляем 3 64-битные умножения: (ah + al) * (bh + bl), ah *bh и al *bl
                                    и радуемся :)

                                    Кстати, алгоритм из статьи (оригинальный) можно найти в «Алгоритмических Трюках для Программиста», там целая глава про умножения была, если мне память не изменяет.
                                    • 0
                                      Спасибо, поправил статью!
                                      • 0
                                        Если бы кто для деления чисел 128-битных через 64-битные алгоритм подсказал.
                                        • +1
                                          Быстрое деление делается через быстрое нахождение обратного и быстрое умножение. Быстрое нахождение обратного делается через метод Ньютона решения уравнения 1/x=a, для чего используется быстрое умножение (точнее возведение в квадрат). Практически вы никакого выигрыша для 128-битных чисел не получите.

                                          Медленное деление вполне через 64-битные числа делается — нужно выделять старшие биты, делить их друг на друга, вычитать и повторять снова с уменьшенным делимым.
                                          • 0
                                            У меня дело не в выигрыше)) Я просто не могу пока разделить... Не могу)) Как-то запутанно с битами возиться. Хотелось бы опираться на 64-битные числа.
                                          • +1
                                            я помню только тупой «в столбик» :) реализуется через сдвиги и вычитания в цикле.
                                            делить 256 бит на 128 (ну или 128 на 128) с помощью деления на 128 на 64 бита я не уверен что можно… надо думать :)
                                            • 0
                                              Выйдет чуть ли не О(N^2), где N — число бит...
                                              Это весьма печально. Хочется пользоваться 64-мя битами и их делением…
                                          • +1
                                            Отличная книга. По английки зовётся Hacker's Delight, тут можно скачать исходники алгоритмов из книги.
                                      • +1
                                        Это ведь нужно для умножение очень преочень больших чисел, как я понимаю? А как можно попробовать умножить реально большие числа, а не 1000 на 2000? я что-то не очень понял. Пробовал так:
                                            ...
                                            typedef unsigned __int64   u8;
                                            typedef unsigned __int64   u16;
                                            typedef unsigned __int64   u32;
                                            u16 a = 17446744073709551615;
                                            ...
                                        


                                        Но Dev-C++, например, говорит, что константа слишком велика. Не подскажите каким образом можно использовать такое большое число?
                                        • +1
                                          Попробуйте так: 17446744073709551615L
                                          И потом непонятно как вы будете печатать 128-битное число на экран…
                                          • +1
                                            Вы с типами что-то не то сделали напишите 1 строку:
                                            typedef unsigned __int64 u64;
                                            • +1
                                              ммм, да. Спасибо большое за помощь
                                              • +1
                                                А зачем вообще все эти typedef? В стандартном заголовочном файле inttypes.h есть типы данных intX_t/uintX_t, включая int64_t/uint64_t.
                                                • –1
                                                  Это если код только под Linux.
                                                  С моим typedef в Visual Studio тоже работает.
                                                  • +1
                                                    Вообще-то, по стандарту C99 положено иметь этот заголовочный файл на любой платформе, равно как и stdint.h. Тем более, по приведённой вами ссылке даже написано, откуда взять файлы для Visual Studio: stdint.h и inttypes.h.
                                                    Не берите пример с Microsoft, она как всегда — на стандарты положила болт, а в результате уже понаписаны тонны непортабельного говнокода, килограммы которого я разгребал на предыдущей работе.
                                                    • –1
                                                      И как-то сложилось у меня именовать типы коротко через i8, i16, i32, i64, u8, u16, u32, u64.
                                            • +1
                                              Спасибо. К сожалению, не работает с L. Как выводить это другой вопрос. Просто хотел попробовать использовать реализованный вами алгоритм.
                                              • 0
                                                У последнего параметра тип должен быть — указатель на половину результата.
                                                Если хотите перемножить 2 64-битных числа, типы дожны быть (u32*,u32*,u64*). Тогда перемножатся два 64-битных числа и результат будет в 128-битном.
                                              • +1
                                                Следует учесть, что когда берётся разность, то промежуточно надо на 1 бит больше для хранения. То есть получается что-то типа 128*128 => 255 бит, а не 256.
                                                Проверьте на предельные значения.
                                                • +1
                                                  Код проверялся полным перебором u16*u16=>u32
                                                  Там этот ньюанс учтён через так сказать if/increment
                                                • 0
                                                  Также для данной задачи можно использовать FFT. А используя для хранения чисел массивы, можно получить числа любой длины (ограниченной только объемом свободной памяти). Причем, операции выполняются довольно быстро — перемножение двух 100К значных чисел без оптимизации у меня выполнялось примерно 7с ( данное время можно уменьшить до 1.5 — 2 с).

                                                  Интересный факт: при использовании BigDecimal и BigInteger платформы Java, потребовалось около 3 минут (скорее всего, там реализовано просто вычитание/сложение «в столбик»).
                                                  • +1
                                                    При первом же взгляде на статью бросается в глаза картинка с большой буквой «Х».
                                                    Мысль — что статья про проект XNeur.

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