image

Деление — одна из самых дорогих операций в современных процессорах. За доказательством далеко ходить не нужно: Agner Fog[1] вещает, что на процессорах Intel / AMD мы легко можем получить Latency в 25-119 clock cycles, а reciprocal throughput — 25-120. В переводе на Русский — МЕДЛЕННО! Тем не менее, возможность избежать инструкции деления в вашем коде — есть. И в этой статье, я расскажу как это работает, в частности в современных компиляторах(они то, умеют так уже лет 20 как), а также, расскажу как полученное знание можно использовать для того чтобы сделать код лучше, быстрее, мощнее.

Собственно, я о чем: если делитель известен на этапе компиляции, есть возможность заменить целочисленное деление умножением и логическим сдвигом вправо (а иногда, можно обойтись и без него вовсе — я конечно про реализацию в Языке Программирования). Звучит весьма обнадеживающе: операция целочисленного умножения и сдвиг вправо на, например, Intel Haswell займут не более 5 clock cycles. Осталось лишь понять, как, например, выполняя целочисленное деление на 10, получить тот же результат целочисленным умножением и логическим сдвигом вправо? Ответ на этот вопрос лежит через понимание… Fixed Point Arithmetic (далее FPA). Чуть-чуть основ.

При использовании FP, экспоненту (показатель степени 2 => положение точки в двоичном представлении числа) в числе не сохраняют (в отличие от арифметики с плавающей запятой, см. IEE754), а полагают ее некой оговоренной, известной программистам величиной. Сохраняют же, только мантиссу (то, что идёт после запятой). Пример:

$0.1 = .0001 1001 1001 1001 ( 1001 )... FP, exp = 0$


0.1 — в двоичной записи имеет 'бесконечное представление', что в примере выше отмечено круглыми скобками — именно эта часть будет повторяться от раза к разу, следуя друг за другом в двоичной FP записи числа 0.1.

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

Пусть дано 16-битное целое число A и 16-битная Fraction часть числа B. Произведение A на B результатом дает число с 16 битами в целой части и 16-тью битами в дробной части. Чтобы получить только целую часть, очевидно, нужно сдвинуть результат на 16 бит вправо.

Поздравляю, вводная часть в FPA окончена.

Формируем следующую гипотезу: для выполнения целочисленного деления на 10, нам нужно выполнить умножение Числа Делимого на FP представление числа 0.1, взять целую часть и дело в шля… минуточку… А будет ли полученный результат точным, точнее его целая часть? — Ведь, как мы помним, в памяти у нас хранится лишь приближенная версия числа 0.1. Ниже я выписал три различных представления числа 0.1: бесконечно точное представление числа 0.1, обрезанное после 16-ого бита без округления представление числа 0.1 и обрезанное после 16 ого бита с округлением вверх представление числа 0.1.

$ 0001\:1001\:1001\:1001\:|\:1001\:1001 .... - infinity\:precision\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\\ 0001\:1001\:1001\:1001\:|\:0000\:0000 .... - truncating\:without\:rounding\\ 0001\:1001\:1001\:1010\:|\:0000\:0000 .... - truncating\:with\:rounding\:up\:$


Оценим погрешности truncating представлений числа 0.1:

$infinity\:precision - truncating\:without\:rounding = 0.6 * 2^{-16}\\ truncating\:with\:rounding\:up - infinity\:precision = 0.1 * 2^{-14}$


Чтобы результат умножения целого числа A, на Аппроксимацию числа 0.1 давал точную целую часть, нам нужно чтобы:

$IntegerPart( A * 0.1 ) = IntegerPart( A * ( 0.1 + 0.1 * 2 ^{-14} ) ), $

, либо

$IntegerPart( A * 0.1 ) = IntegerPart( A * ( 0.1 + 0.6 * 2 ^{-16} ) ) $


Удобнее использовать первое выражение: при $0.1 * 2 ^ {-14} * A < 0.1 $ мы всегда получим тождество (но, заметь, не все решения, что в рамках данной задачи — более чем достаточно). Решая, получаем $A < 2^{14}$. То есть, умножая любое 14-битное число A на truncating with rounding up представление числа 0.1, мы всегда получим точную целую часть, которую бы получили умножая бесконечно точно 0.1 на A. Но, по условию у нас умножается 16-битные число, а значит, в нашем случае ответ будет неточным и довериться простому умножению на truncating with rounding up 0.1 мы не можем. Вот если бы мы могли сохранить в FP представлении числа 0.1 не 16 бит, а, скажем 19, 20 — то все было бы ОК. И ведь можем же!
Внимательно смотрим на двоичное представление — truncating with rounding up 0.1: старшие три бита — нулевые, а значит, и никакого вклада в результат умножения не вносят (новых бит).
Следовательно, мы можем сдвинуть наше число влево на три бита, выполнить округление вверх и, выполнив умножение и логический сдвиг вправо сначала на 16, а затем на 3 (то есть, за один раз на 19 вообще говоря) — получим нужную, точную целую часть. Доказательство корректности такого '19' битного умножения аналогично предыдущему, с той лишь разницей, что для 16-битных чисел оно работает правильно. Аналогичные рассуждения верны и для чисел большей разрядности, да и не только для деления на 10.

Ранее я писал, что вообще говоря, можно обойтись и без какого-либо сдвига вовсе, ограничившись лишь умножением. Как? Ассемблер x86 / x64 на барабане:
В современных процессорах, есть команда MUL (есть еще аналоги IMUL, MULX — BMI2), которая принимая один, скажем 32 / 64 -битный параметр, способна выполнять 64 / 128 битное умножение, сохраняя результат частями в два регистра (старшие 32 / 64 бита и младшие, соответственно):

MUL RCX	; умножить RCX на RAX, а результат (128 бит) сохранить в RDX:RAX

В регистре RCX пусть хранится некое целое 62-битное A, а в регистре RAX пускай хранится 64 битное FA представление truncating with rounding up числа 0.1 (заметь, никаких сдвигов влево нету). Выполнив 64-битное умножение получим, что в регистре RDX сохранятся старшие 64 бита результата, или, точнее говоря — целая часть, которая для 62 битных чисел, будет точной. То есть, сдвиг вправо (SHR, SHRX) не нужен. Наличие же такого сдвига нагружает Pipeline процессора, вне зависимости поддерживает ли он OOOE или нет: как минимум появляется лишняя зависимость в, скорее всего и без того длинной цепочке таких зависимостей (aka Dependency Chain). И вот тут, очень важно упомянуть о том, что современные компиляторы, видя выражение вида some_integer / 10 — автоматически генерируют ассемблерный код для всего диапазона чисел Делимого. То есть, если вам известно что числа у вас всегда 53-ех битные (в моей задаче именно так и было), то лишнюю инструкцию сдвига вы получите все равно. Но, теперь, когда вы понимаете как это работает, можете сами с легкостью заменить деление — умножением, не полагаясь на милость компилятору. К слову сказать, получение старших битов 64 битного произведения в C++ коде реализуется интринсиком (something like mulh), что по Asm коду, должно быть равносильно строчкам инструкции {I}MUL{X} выше.

Возможно, с появлением контрактов (в С++ 20 не ждем) ситуация улучшится, и в каких-то кейсах, мы сможем довериться машине! Хотя, это же С++, тут за все отвечает программист — не иначе.

Рассуждения описанные выше — применимы к любым делителям константам, ну а ниже список полезных ссылок:

[1] https://www.agner.org/optimize/instruction_tables.pdf
[2] Круче чем Агнер Фог
[3] Телеграмм канал, с полезной информацией о оптимизациях под Intel / AMD / ARM
[4] Про деление нацело, но по Английски

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


  1. Evengard
    23.09.2019 23:41

    Правильно ли я понимаю, что «магические числа» вроде 0x92492492 — это на самом деле hex представление этих самых fixed point чисел?


    1. IvanKamynin Автор
      24.09.2019 00:33

      Именно. В моей статье, достаточно просто взять binary rep. описанный и перевести в hex. Самое главное — магии нету никакой!


  1. Nagg
    24.09.2019 01:05
    +1

    Кстати, в Intel Ice Lake уменьшили латенси у idiv более чем в 5 раз до 18 циклов (не помню на счет пропускной способности) — а это уже не так страшно как ранее.


  1. assembled
    24.09.2019 08:07
    +1

    Это случайно не перевод? Уж слишком статья пестрит бессмысленными применениями английских терминов или даже фраз:

    можем получить Latency в 25-119 clock cycles, а reciprocal throughput — 25-120.

    Fraction часть числа B.

    умножению на truncating with rounding up 0.1

    Выглядит weird и unnaturaly.
    Особенно это:
    получение старших битов 64 битного произведения в C++ коде реализуется интринсиком (something like mulh)


  1. IgorPie
    24.09.2019 11:04

    Тоже открывал эту Америку, что умножение — 1 такт, деление — 14, и велосипед изобретал.
    А вот авторы gcc это знали давно, поэтому наколхоживание псевдоумножений в коде ничего не дало, только испортило читаемость.

    Если не просят оптимизировать и ускорять, лучше не усложнять. Или сразу считать в q63.


  1. antoshkka
    24.09.2019 13:50

    Спасибо, давно искал эту тему на русском!

    Но, теперь, когда вы понимаете как это работает, можете сами с легкостью заменить деление — умножением, не полагаясь на милость компилятору

    … и такие исходники, которые junior и middle разработчики не поймут и не смогут поддерживать.

    Вот так, на ровном месте, можно себе и компании создать проблемы.


    1. IvanKamynin Автор
      24.09.2019 14:14

      Ну, сложно сказать — если контора математическая, код документируется, то вроде должно быть ОК. Собственно, в Алгоритме Ruy(printf) — сделано умножением, правда, со сдвигами. То есть, сий инженер из g00gl3 по факту написал тоже самое, что делают компиляторы видя / 10, оставив все теже сдвиги :]


    1. yurmikh
      24.09.2019 20:38
      +2

      Тема давно изложена на русском языке. Уоррен Г. «Алгоритмические трюки для программистов» bookfi.net/dl/437262/99b72b

      Ну а прием этот эффективно работает далеко не всегда.


    1. yurmikh
      24.09.2019 21:18

      И да, «магические числа» — это из Уоррена. Так что Автор читал, но умолчал об источнике.


      1. IvanKamynin Автор
        24.09.2019 21:30

        Не читал, если что. А смысл статьи больше лежит в последних абзацах. Но за источник — спасибо, я себя открывателем не называл!


  1. rebuilder
    24.09.2019 14:39

    Простите, но я не понял о чём статья:

    …если делитель известен на этапе компиляции, есть возможность заменить целочисленное деление умножением и логическим сдвигом вправо (а иногда, можно обойтись и без него вовсе — я конечно про реализацию в Языке Программирования)
    Если у меня в коде есть выражение X/10, что мне мешает просто заменить его на X*0.1?


    1. Cerberuser
      24.09.2019 14:47

      То, что операция превращается из целочисленной во floating-point, вероятно.


      1. IvanKamynin Автор
        24.09.2019 15:40

        Как-то так, да. А это, точно приводит к потере точности для числе больших 2^53.


        1. qw1
          24.09.2019 17:04

          Почему, если результат умножения 128-битный?
          UPD. А, это если заменить целочисленные расчёты на FP…