Умножение сдвигами и деление на константу умножением

Если ты разложишь умножение на 10 на сдвиги, объяснишь, откуда берётся число 43691 при делении на три и почему его округляют вверх, а не вниз, - листай до «Теперь сам».

Как думаешь, есть ли в скомпилированной программе деление

В программе на C написано x / 10 - например, чтобы напечатать число по цифрам. Процессор умеет делить, у него есть для этого команда.

Прежде чем читать дальше, ответь себе: будет ли эта команда в машинном коде? Не будет. Компилятор заменит деление умножением на странную константу 0xCCCCCCCD и сдвигом вправо. Откуда такая константа и почему результат при этом точный - предмет урока.

На что это похоже

Разделить на пять в уме проще, если умножить на два и перенести запятую: 135 / 5 = 135 · 2 / 10 = 27. Деление заменено умножением на число, обратное делителю, - 1/5 = 0,2 - и «сдвигом» на десятичный разряд, который ничего не стоит.

С делением на три так гладко не выходит: 1/3 = 0,3333…, дробь бесконечная. Но если нужна только целая часть ответа и делимое не слишком велико, бесконечная точность не нужна: 0,3334 даст тот же целый ответ, что и точная треть, для всех делимых до нескольких тысяч. Вопрос только в том, сколько знаков оставить и в какую сторону округлить.

Механизм: обратное число, масштаб и округление

Умножение на константу. Константа в двоичной записи - сумма степеней двойки, а умножение на степень двойки - сдвиг. 10 = 8 + 2, значит x · 10 = (x << 3) + (x << 1). Блок единиц выгоднее записать разностью: 7 = 8 − 1, x · 7 = (x << 3) − x. Это умножение столбиком, в котором множитель известен заранее и нулевые строки просто не написаны.

Деление на степень двойки - сдвиг вправо. Для отрицательных чисел он округляет вниз, а деление в C - к нулю, поэтому компилятор сначала прибавляет поправку 2ᵏ − 1 к отрицательным - маской из знака, без ветвления: (x + (x >> 31 & (2ᵏ − 1))) >> k.

Деление на произвольную константу d. Заранее считают M ≈ 2ˢ / d - обратное число в масштабе 2ˢ. Тогда x / d ≈ x · M / 2ˢ, а деление на 2ˢ - сдвиг:

d = 3,  s = 17:      2¹⁷ / 3 = 43690,67        M = 43691   (округление вверх)

x · M >> 17   при x = 1000:      1000 · 43691 = 43 691 000,   >> 17  →  333

Почему вверх. M немного больше точной дроби: 43691 = (2¹⁷ + 1) / 3. Значит

x · M / 2¹⁷  =  x / 3  +  x / (3 · 2¹⁷)
                └ точный ответ    └ лишнее: меньше 1/3, пока x < 2¹⁷

Дробная часть точного ответа x / 3 бывает 0, 1/3 или 2/3. Лишнее меньше 1/3, и до следующего целого сумма не дотягивает - целая часть остаётся верной. Если округлить M вниз, ошибка станет отрицательной, и там, где x / 3 - ровно целое, результат окажется чуть меньше него и округлится на единицу вниз: так ломается уже x = 3.

Чем больше сдвиг s, тем точнее M и тем для больших x формула верна. Для 16-битных x при делении на 3 хватает s = 17, и множитель помещается в 16 битов. Для 32-битных - s = 33, множитель 0xAAAAAAAB, а для деления на 10 - s = 35 и 0xCCCCCCCD, та самая константа из начала урока. Если x и M помещаются в w битов, произведению хватает 2w битов, а сдвиг берёт его старшую часть. Некоторым делителям нужен множитель шире w битов; тогда прямой формуле нужна ещё более широкая арифметика либо преобразование со сложением и сдвигом.

Правило

Умножение на константу - сумма или разность сдвигов по двоичной записи константы. Деление беззнакового x на константу d: x * M >> s, где M = ⌈2ˢ / d⌉, а s выбрано так, чтобы лишнее x · (M · d − 2ˢ) / (d · 2ˢ) оставалось меньше 1/d для всех x нужной ширины. Умножение должно сохранить все биты, нужные после сдвига: если M шире слова, двойной ширины может не хватить. Остаток - x − (x // d) · d.

Разбор: деление на 3 и на 10, проверенное на всех числах

x = 1234
print(x * 10, (x << 3) + (x << 1), x * 7, (x << 3) - x)
print(1000 * 43691 >> 17, 65535 * 43691 >> 17, 65535 // 3)
print(3 * 43690 >> 17, 3 * 43691 >> 17)
print(1234 * 52429 >> 19, 1234 - (1234 * 52429 >> 19) * 10)

12340 12340 8638 8638
333 21845 21845
0 1
123 4

Последняя строка - деление на 10: 52429 = ⌈2¹⁹ / 10⌉, это 0xCCCD. Частное 123 и остаток 4: так функция печати отщепляет от числа очередную цифру.

Один пример ничего не доказывает, поэтому - проверка на всех 16-битных числах. Множитель M, сдвиг S и делитель d задаются перед запуском; программа печатает число ошибок и первое число, на котором формула ошиблась, или −1. Четыре прогона сведены в таблицу:

oshibok = 0
pervaya = -1
for x in range(65536):
    if x * M >> S != x // d:
        oshibok = oshibok + 1
        if pervaya < 0:
            pervaya = x
print(oshibok, pervaya)

делитель   множитель M      сдвиг S    ошибок    первая ошибка
   3       43691 (вверх)      17           0          -
   3       43690 (вниз)       17      21 845          3
   3       21846 (вверх)      16      10 923     32 768
  10       52429 (вверх)      19           0          -

Множитель, округлённый вниз, ошибается на каждом числе, кратном трём, кроме нуля. Множитель с масштабом на бит меньше верен для первой половины диапазона и ломается ровно на 2¹⁵: лишнее доросло до трети.

Разбор: делитель 7, обратный по модулю и граница приёма

С семёркой гладко не выходит.

делитель   множитель M      сдвиг S    ошибок    первая ошибка
   7       37450 (вверх)      18       3 121     43 693
   7       74899 (вверх)      19           0          -

Верный множитель - 74 899, а это 0x12493: 17 битов. В 16-битное слово он не помещается, и произведение x · M не помещается в двойное слово. То же в 32 битах: для семёрки множитель занимает 33 бита. Компилятор выходит из положения лишним сложением и сдвигом - берёт множитель без старшего бита, а недостающее x · 2¹⁶ добавляет отдельно. Делители 3, 5 и 10 «удобные», 7 - первый «неудобный».

У числа 43691 есть второе имя - из раздела про числа. 3 · 43691 = 131 073 = 2¹⁷ + 1, значит по модулю 2¹⁶ это единица: 43691 - обратный элемент к тройке по модулю 65 536. Отсюда третий приём - точное деление. Если заранее известно, что x делится на 3 нацело, частное - это x · 43691 по модулю 2¹⁶, без сдвига и без старшей половины произведения:

print(43691 * 3 & 0xFFFF, 369 * 43691 & 0xFFFF, 124 * 43691 & 0xFFFF, 0xFFFF // 3)

1 123 43732 21845

369 = 3 · 123, и умножение на обратный вернуло 123. А 124 на три не делится - и результат 43 732 больше наибольшего возможного частного 21 845. Это готовая проверка делимости без деления: умножить на обратный и сравнить с границей.

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

Что ломается без этого

Первое: младшая половина вместо старшей. x * M для 32-битных чисел - 64-битное произведение, и нужны его старшие биты. В языке, где произведение двух 32-битных чисел обрезается до 32 битов, формула молча возвращает мусор. В тренажёре по той же причине пример взят 16-битный: 0xAAAAAAAB на 32-битное число - это до 2⁶⁴, за пределом точного счёта тренажёра 2⁵³.

Второе: формула для беззнаковых, число со знаком. Для отрицательного x нужна поправка на единицу - иначе округление пойдёт вниз, а не к нулю. Компилятор добавляет её сам; при ручном переносе формулы о ней забывают.

Третье: остаток через % после деления. Написать x / 10 и x % 10 рядом - нормально: компилятор сделает одно умножение и вычислит остаток как x − q · 10. А вот заменять % на вычитание в цикле («отнимай 10, пока больше») - превращать одну команду в тысячу.

Четвёртое: сдвиг вместо деления отрицательного. Из главы про дополнительный код: -7 >> 1 равно −4, а -7 / 2 в C равно −3. С поправкой (x + (x >> 31 & 1)) >> 1 получается −3.

Зачем это в работе

Это происходит в каждой программе на компилируемом языке, просто без твоего участия. Печать числа - деление на 10 в цикле; перевод времени из секунд в часы - деление на 3600; индекс в хеш-таблице простого размера - остаток по константе. Во всех трёх местах в машинном коде стоит умножение.

Из этого следует практический совет: если делитель известен, пиши его константой, а не читай из конфигурации в переменную. Константу компилятор превратит в умножение, переменную - нет. По той же причине размер таблицы, равный степени двойки, быстрее простого: маска против умножения со сдвигом, а оба они быстрее настоящего деления.

Генераторы случайных чисел и хеш-функции используют умножение на подобранную константу по модулю 2³² или 2⁶⁴ как дешёвое перемешивание; обратимо оно ровно потому, что у нечётной константы есть обратный элемент.

Теперь сам

1. Разложи x * 12 и x * 31 на сдвиги так, чтобы слагаемых было два.
2. Найди множитель для деления на 5 при s = 18 и проверь его на x = 1000.
3. Почему при делении на 3 множитель 43690 ломается именно на x = 3?
4. M = 21846, s = 16. Чему равно лишнее для x = 32 768 и почему на этом числе формула ошибается?
5. Число 369 умножили на 43691 и взяли 16 младших битов - вышло 123. Что получится для 370 и о чём это говорит?