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

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

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

При реализации сложения больших чисел (например, от 2⁶⁴ и выше) обычно пишут код, похожий на этот алгоритм. Любопытно здесь то, что существует простой трюк, позволяющий существенно ускорить этот процесс на современных 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-битными числами. Можно рассматривать эти части как «разряды» по основанию 2⁶⁴, потому что каждый «разряд» имеет значение от 0 до 2⁶⁴ - 1 включительно.

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

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

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

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

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

В нашем примере по основанию 10 применение 37 цифр позволяло нам складывать до четырёх нормализованных чисел без необходимости распространения переносов. В radix 2⁵¹ representation применение 2⁶⁴ цифр позволяет нам складывать до 2¹³ нормализованных чисел, не беспокоясь о переполнении старших 13 битов.

Примечание: почему 13 битов вместо 12? В нашем случае мы проигнорируем переносы в самом старшем limb, позволив числам вернуться к нулю при переполнении выше 2²⁵⁶ - 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 2⁵¹ обгоняет по скорости сложение radix 2⁶⁴ на моём CPU Haswell даже для всего трёх сложений; и это уже с учётом преобразования из radix 2⁵¹ representation. Масштабы увеличения скорости повышаются с ростом количества сложений.

Вычитание

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

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

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

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