Помните, как долго выполняется сложение на бумаге?

 ¹¹ ¹
  6876
+ 3406
------
 10282

Начиная с «единиц», мы складываем 6 + 6 = 12, записываем 2 и переносим 1. Затем пошагово двигаемся влево, пока складываемые разряды не закончатся.

При реализации сложения больших чисел (например, от 264 и выше) обычно пишут код, похожий на этот алгоритм. Любопытно здесь то, что существует простой трюк, позволяющий существенно ускорить этот процесс на современных CPU.

Но сначала я задам вопрос: почему сложение столбиком мы начинаем с самого младшего разряда? Почему бы не начать слева?

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

Представьте, что мы бы складывали слева направо:

6 + 3 = 9. Первый разряд — это 9.
8 + 4 = 12. Так, второй разряд равен 2… но нужно перенести 1, поэтому первый разряд на самом деле оказывается равным 9 + 1 = 10… а теперь перенесём обратно эту единицу…

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

Ну а что насчёт компьютеров?

Разумеется, компьютеры не работают с числами по основанию 10. У современных десктопных и серверных CPU есть интерфейс для работы с 64-битными числами (по большей мере).

; Прибавляем 64-битное значение в B к 64-битному значению в A
add A, B
; Примечание: чтобы не усложнять, я буду использовать буквы вместо имён регистров

Пока наши числа умещаются в одно 64-битное значение, всё просто. Но что, если нам нужно, например, сложить два 256-битных целых числа x и y?

Очевидным решением было бы разбиение каждого 256-битного числа на четыре 64-битных части (обычно называемых «limb»). Поместим старшие 64 бита числа x в регистр A, следующие 64 бита в регистр B и аналогично для регистров C и D. Проделаем то же самое для числа y с регистрами E, F, G, H.

Теперь можно сложить x и y  сложением соответствующих частей:

; Эквивалент x += y
add A, E
add B, F
add C, G
add D, H

Но постойте, ведь так мы можем получить неверный результат! Если хотя бы в одном из трёх последних сложений произойдёт переполнение, то нам придётся «перенести» эту дополнительную единицу вверх, в следующую 64-битную часть. Звучит знакомо, не так ли?

К счастью, в x86 есть для этого специальная команда, называющаяся «сложением с переносом» (add with carry). adc автоматически проверяет предыдущую операцию на переполнение и при необходимости прибавляет 1. Вот, как будет выглядеть правильный код:

add D, H
adc C, G ; включаем перенос из предыдущей операции
adc B, F ; включаем перенос из предыдущей операции
adc A, E ; включаем перенос из предыдущей операции

Как и в случае со сложением в столбик по основанию 10, мы начинаем с самых младших «разрядов» (D и H) и постепенно поднимаемся до самых старших (A и E), по пути перенося 1, если это необходимо.

Но в таком виде это медленно

Любопытно, что наш исправленный код медленнее, чем изначальный (неправильный) код. Гораздо медленнее. Но почему?

Во-первых, в большинстве популярных CPU x86 adc просто выполняется медленнее, чем обычная add. Так как у adc есть третье входящее значение (флаг переноса), эта команда сложнее, чем add. Кроме того, она гораздо реже используется, чем add, поэтому у проектировщиков CPU гораздо меньше мотивации тратить площадь чипов на оптимизацию скорости adc.

Вторая причина интереснее. Давайте рассмотрим в качестве примера микроархитектуру Intel Haswell.

В CPU Haswell для выполнения одной команды add требуется один такт. Однако в идеальных условиях CPU Haswell могут выполнять за такт до четырёх команд add. Как такое возможно? Благодаря параллелизму. Современные процессоры смотрят, какие команды идут дальше, и пытаются запланировать их таким образом, чтобы они по возможности выполнялись параллельно. Так как CPU Haswell имеют 8 портов выполнения, и 4 из этих портов могут выполнять целочисленную команду add, процессор Haswell может одновременно выполнять до 4 команд add.

В нашем изначальном коде сложения все 4 команды add независимы друг от друга, поэтому процессору было легко их распараллелить. В случае же с adc каждая команда зависит от результата предыдущей команды. Процессор не имеет выбора и должен выполнять команды только последовательно, одну за другой.

Разница в скорости ещё более существенна, если мы работаем с командами SIMD (Single Instruction, Multiple Data). Например, одна команда vpaddq (Vector Packed Add Quadword) одновременно выполняет четыре 64-битных сложения. Если учесть также тот факт, что процессоры Haswell могут выполнять по две команды vpaddq за такт, можно понять, для для правильной обработки нам приходится идти на серьёзное снижение скорости.

Устраняем переносы, часть первая: на бумаге

Вернёмся на минуту к основанию 10. Как нам избавиться от необходимости переносов?

Давайте внесём изменения в то, как устроена числовая система. Для начала расширим интервал доступных значений разрядов. Вместо 0-9 будем использовать 0-9, A-Z и *:

0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ*

(Да, чтобы всё работало правильно, мне понадобился этот дополнительный символ. Подождите, скоро всё станет понятно.)

Хоть у нас и есть 37 значений, мы не используем основание 37. У чисел по-прежнему есть разряды единиц, десятков и сотен, как в обычной системе по основанию 10. Число 29 по-прежнему означает 29, а 29 + 1 по-прежнему равно 30. Единственное различие заключается в том. что разрядами мы можем выполнять счёт выше 9: 29 + 1 можно записать как 2A, 1K или даже U.

В этой новой, более гибкой числовой системе мы можем выполнять сложение без необходимости в переносах!

  672415
+ 736606
--------
  DA8A1B

Этот трюк сработает не для всех чисел нашей системы (например, для 9 + W потребуется перенос), но он будет работать, если складываемые числа нормализованы, то есть все значения их разрядов не превышают 9. На самом деле, в этой записи мы можем сложить до четырёх нормализованных чисел без переносов:

  999   <-- самое большое нормализованное трёхразрядное число
  999
  999
+ 999
-----
  ***   <-- валидный трёхразрядный результат, переносов не требуется
            (помните, что * - это наибольшее значение разряда)

То есть внеся хитрые изменения в числовую систему, мы сжульничали и избавились от части переносов. Разумеется, рано или поздно нам придётся выполнить преобразование из 37-цифровой системы по основанию 10 в обычную систему по основанию 10. Мы можем это сделать, нормализовав число так, чтобы каждый из его разрядов находится в интервале от 0 до 9:

  ¹¹ ¹ ¹
   DA8A1B
= 1409021

объяснение:
D = 10 + 3
A = 10 + 0
B = 10 + 1

Мы нормализуем число, начиная справа, определяя, сколько «десяток» в каждом разряде, вычитая эти «десятки» и перенося их в следующий разряд. 672415 и 736606 действительно дают в сумме 1409021, так что система работает!

Самое важное здесь то, что можно использовать эту методику, чтобы откладывать распространение переносов до самого конца. Мы не можем полностью от него избавиться, но можем избежать его временно. Если мы сэкономим на переносах, происходящих при промежуточных сложениях, то можем распространить их все в конце за один проход.

Устраняем переносы, часть вторая: компьютеры

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

Раньше мы разделяли 256-битное число на четыре 64-битные части, потому что процессоры x86_64 работают с 64-битными числами. Можно рассматривать эти части как «разряды» по основанию 264, потому что каждый «разряд» имеет значение от 0 до 264 - 1 включительно.

При работе с основанием 10 мы сохраняли то же самое основание, но расширяли интервал значений разряда, чтобы предотвратить возникновение переносов. К сожалению, в этом случае так поступить не удастся – в 64-битное число не уместятся значения больше, а аппаратную часть мы менять не можем. Но можно добиться того же результата, уменьшив размер основания.

Вместо того, чтобы делить 256 битов на четыре разряда по основанию 264, мы разобьём 256 битов на пять разрядов по основанию 251. Каждый разряд по-прежнему может иметь значения в интервале от 0 до 264 - 1, однако меньшее основание обеспечивает нам гибкость, необходимую, чтобы избавиться от необходимости переноса. В литературе по криптографии эту методику обычно называют «radix 251 representation».

Вот, как это будет выглядеть, когда мы разобьём 256 битов на пять частей (то есть разрядов):

|            [--------------------- 52 бита --------------------]|
|             [-------------------- 51 бит  --------------------]|
|             [-------------------- 51 бит  --------------------]|
|             [-------------------- 51 бит  --------------------]|
|             [-------------------- 51 бит  --------------------]|

Каждый разряд (limb) содержит 51 (или 52) битов исходного 256-битного числа. Оставшиеся 12 или 13 битов образуют дополнительные «цифры», которые нужны нам для предотвращения переносов. По сути, самые старшие биты каждого limb резервируются под хранилище любых переносов, которые происходят во время вычислений.

В нашем примере по основанию 10 применение 37 цифр позволяло нам складывать до четырёх нормализованных чисел без необходимости распространения переносов. В radix 251 representation применение 264 цифр позволяет нам складывать до 213 нормализованных чисел, не беспокоясь о переполнении старших 13 битов.

Примечание: почему 13 битов вместо 12? В нашем случае мы проигнорируем переносы в самом старшем limb, позволив числам вернуться к нулю при переполнении выше 2256 - 1 (точно так же ведёт себя беззнаковое сложение в C с целочисленными типами обычного размера). Поэтому мы можем назначить 52 бита самому старшему limb и забыть о том, что в нём закончится место для переносов раньше, чем в других limb.

В этом новом представлении код сложения будет выглядеть так:

; Пусть число x разбито на A, B, C, D, E (A = самый старший)
; а число y разбито на F, G, H, I, J (F = самый младший)
add A, F
add B, G
add C, H
add D, I
add E, J
; Параллелизм, ура!

Несмотря на то, что теперь нам нужно пять команд add вместо четырёх, благодаря отсутствию переносов сложение происходит гораздо быстрее.

Разумеется, нам также понадобится код для нормализации числа распространением переносов.

; Пусть число x разбито на A, B, C, D, E (A = самый старший)
; Регистр T используется в качестве временного хранилища

mov T, E ; Копируем E в T
shr T, 51 ; Сдвигаем всё, кроме переносов
add D, T ; Прибавляем переносы из E в D
and E, 0x0007FFFFFFFFFFFF ; Обнуляем переносы в E

mov T, D ; Копируем D в T
shr T, 51 ; Сдвигаем всё, кроме переносов
add C, T ; Прибавляем переносы из D в C
and D, 0x0007FFFFFFFFFFFF ; Обнуляем переносы в D

mov T, C ; Копируем C в T
shr T, 51 ; Сдвигаем всё, кроме переносов
add B, T ; Прибавляем переносы из C в B
and C, 0x0007FFFFFFFFFFFF ; Обнуляем переносы в C

mov T, B ; Копируем B в T
shr T, 51 ; Сдвигаем всё, кроме переносов
add A, T ; Прибавляем переносы из B в A
and B, 0x0007FFFFFFFFFFFF ; Обнуляем переносы в B

and A, 0x000FFFFFFFFFFFFF ; Обнуляем переносы в A

Как ни удивительно, бенчмарки показывают, что сложение методикой radix 251 обгоняет по скорости сложение radix 264 на моём CPU Haswell даже для всего трёх сложений; и это уже с учётом преобразования из radix 251 representation. Масштабы увеличения скорости повышаются с ростом количества сложений.

Вычитание

Пока мы рассматривали только сложение, однако очень легко адаптировать эту методику и для вычитания. Главное отличие сложения от вычитания заключается в том, что вычитание имеет отрицательные переносы.

Выше мы обращались со всеми limb (и их переносами) как с беззнаковыми целыми числами. Для поддержки вычитания можно обращаться с limb как со знаковыми целыми числами, позволяющими разрядам быть и положительными, и отрицательными. После внесения таких изменений можно хранить и положительные, и отрицательные переносы.

Побочным результатом этого станет то, что самый старший бит каждого limb будет резервироваться под знак. Это снижает количество операций, которое мы можем выполнять между нормализациями, с 213 до 212; в большинстве случаев это небольшая жертва.

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

Комментарии (13)


  1. kvd264
    05.06.2025 08:14

    А что произойдет, если мы случайно сложим больше, чем 2^13 чисел, произойдет наихудший сценарий и переносы переполнятся? Это же как-то контролируется?


    1. jin_x
      05.06.2025 08:14

      Результат будет неверным, только и всего.


  1. domix32
    05.06.2025 08:14

    Как-то недожато. Взяли б из питона тогда имплементацию. Там числа как раз представляются в СС с большим основанием, а потом ещё и алгоритмом карацубы все это складывается, а не наивным сложением.


    1. aamonster
      05.06.2025 08:14

      Алгоритм Карацубы – это про умножение. Ну и думаю, что там сложение простое: оптимизация из статьи – для суммирования кучи чисел... По сути, вынос части работы из цикла суммирования.


      1. domix32
        05.06.2025 08:14

        С карацубой это я погорячился, да. Но суть с представлением та же.


    1. Serge3leo
      05.06.2025 08:14

      Для Python, в основных вариантах использования, это всё теряется на фоне _PyLong_New (PyObject_Malloc), поэтому, для 64-бит архитектур, CPython хранит 30-бит цифры (limb) в uint32_t/int32_t.

      Ну, а насчёт, недожато, то да, тема автоматической нормализации не раскрыта, но это объективно, 256-бит числа взяты же с потолка, для примера.

      В скалярном коде, возможно, потребуется ещё одно, шестое, сложение сложение и один условный переход, а возможно лучше иметь пять условных переходов с достаточно низкой вероятностью перехода, кто знает?

      В векторном коде, неопределённость ещё выше.


  1. Mausglov
    05.06.2025 08:14

    Мы нормализуем число, начиная справа, определяя, сколько «десяток» в каждом разряде, вычитая эти «десятки» и перенося их в следующий разряд. 672415 и 736606 действительно дают в сумме 1409021, так что система работает!

    Получается, проблему переноса разряда просто перенесли в другое место ( извините за тавтологию),и 17 ассемблерных инструкций вместо adc работают быстрее?


    1. tbl
      05.06.2025 08:14

      Получается, проблему переноса разряда просто перенесли в другое место ( извините за тавтологию),и 17 ассемблерных инструкций вместо adc работают быстрее?

      во первых, это концепт: на большом количестве слагаемых это будет выгоднее, чем adc-цепочка из-за instruction-level parallelism, либо если использовать SIMD, то можно использовать векторизацию

      во вторых, команд исполнится на ALU не 17, а 13, т.к. mov-команды здесь исполнятся блоком переименования регистров за 0 тактов: https://agner.org/optimize/microarchitecture.pdf разделы "Register allocation and renaming" различных архитектур


    1. arteast
      05.06.2025 08:14

      Не просто перенесли в другое место, а собрали в одну операцию то, что в наивном варианте надо сделать 4000 раз. Теперь сериализация из-за зависимости по данным происходит не на каждом сложении, а один раз в несколько тысяч сложений.


      1. wataru
        05.06.2025 08:14

        А почему тут сериализация раз в несколько тысяч сложений? В этом коде во время перенорсов вычисление каждого следующего разряда завист от предыдущего. Вот цепочка add B,T mov T,B SHR 51 T - тут зависимость по данным.


        1. jin_x
          05.06.2025 08:14

          Если нужно выполнить 100 сложений длинных 256-битных чисел, то это операцию из 17 инструкций нужно выполнить только один раз в самом конце. Вот и ускорение.


  1. jin_x
    05.06.2025 08:14

    Если бы у автора был процессор следующего поколения (Broadwell), он бы мог ускорить сложение, используя инструкции ADCX и ADOX, с помощью которых как раз можно распределить 2 сложения. Конечно, на сотне сложений это будет работать медленнее, но на небольшом кол-ве даст больший перфоманс.

    А ещё в AVX-512 есть 2 интересные инструкции: VPMADD52LUQ и VPMADD52HUQ, которые выполняют сложение с умножением над 52-битными целыми числами :)


  1. Voldemaar
    05.06.2025 08:14

    Правильно ли я понял, что
    vpaddq %xmm5, %xmm6, %xmm5 - это AVX
    vpaddq %ymm4, %ymm0, %ymm0 - это AVX2

    ?