Почему GCC использует умножение на странное число при реализации целочисленного деления?

Я читал о div и mulсборочных операциях, и я решил , чтобы увидеть их в действии, написав простую программу в C:

Файл деление.c

#include <stdlib.h>
#include <stdio.h>

int main()
{
    size_t i = 9;
    size_t j = i / 5;
    printf("%zu\n",j);
    return 0;
}

И затем генерирование кода на ассемблере

gcc -S division.c -O0 -masm=intel

Но, глядя на сгенерированный division.sфайл, он не содержит никаких операций div! Вместо этого он выполняет какую-то черную магию со сдвигом битов и магическими числами. Вот фрагмент кода, который вычисляет i/5:

mov     rax, QWORD PTR [rbp-16]   ; Move i (=9) to RAX
movabs  rdx, -3689348814741910323 ; Move some magic number to RDX (?)
mul     rdx                       ; Multiply 9 by magic number
mov     rax, rdx                  ; Take only the upper 64 bits of the result
shr     rax, 2                    ; Shift these bits 2 places to the right (?)
mov     QWORD PTR [rbp-8], rax    ; Magically, RAX contains 9/5=1 now, 
                                  ; so we can assign it to j

Что тут происходит? Почему GCC вообще не использует div? Как он генерирует это магическое число и почему все работает?

Целочисленное деление - это одна из самых медленных арифметических операций, которые вы можете выполнять на современном процессоре, с задержкой до десятков циклов и плохой пропускной способностью. (Для x86 см . Таблицы инструкций Agner Fog и руководство microarch ).

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

Реализация /оператора C таким способом, а не с использованием последовательности из нескольких команд, divявляется просто способом GCC по умолчанию для деления на константы. Он не требует оптимизации операций и ничего не меняет даже для отладки. (Использование -Osдля небольшого размера кода действительно заставляет GCC использовать div.) Использование мультипликативного обратного вместо деления похоже на использование leaвместо mulиadd

В результате вы склонны видеть divили idivвыводить, только если делитель не известен во время компиляции.

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

Деление на 5 - это то же самое, что умножение на 1/5, что опять же, как умножение на 4/5 и сдвиг вправо на 2 бита. Соответствующее значение CCCCCCCCCCCCCCCDв шестнадцатеричном формате, которое является двоичным представлением 4/5, если ставится после шестнадцатеричной точки (т. Е. Двоичное для четырех пятых 0.110011001100повторяется - см. Ниже, почему). Я думаю, что вы можете взять это отсюда! Возможно, вы захотите проверить арифметику с фиксированной точкой (хотя обратите внимание, что в конце она округляется до целого числа).

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

Видеть Взаимное умножение, учебное пособие для подробного описания того, как оно работает, объясняя с точки зрения фиксированной точки. Он показывает, как работает алгоритм поиска обратной величины, и как обрабатывать деление со знаком и по модулю.

Давайте на минуту рассмотрим, почему 0.CCCCCCCC...(шестнадцатеричный) или 0.110011001100...двоичный 4/5. Разделите двоичное представление на 4 (сдвиньте вправо на 2 позиции), и мы получим, 0.001100110011...который путем тривиального осмотра может быть добавлен к полученному оригиналу 0.111111111111..., который, очевидно, равен 1, точно так же, как 0.9999999...десятичное число равно единице. Таким образом, мы знаем , что x + x/4 = 1, таким образом 5x/4 = 1, x=4/5. Затем это представляется CCCCCCCCCCCCDв виде шестнадцатеричного числа для округления (поскольку двоичная цифра за пределами последней присутствующей будет a 1).

В общем, умножение намного быстрее, чем деление. Так что, если нам удастся избежать умножения на обратное, мы сможем значительно ускорить деление на константу

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

-3689348814741910323 - 0xCCCCCCCCCCCCCCCD, значение чуть более 4/5, выраженное в 0,64 с фиксированной точкой.

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

Это дает нам хотя бы приблизительное значение деления на 5, но дает ли он точный ответ, правильно округленный до нуля?

Чтобы получить точный ответ, ошибка должна быть достаточно маленькой, чтобы не переместить ответ за границу округления.

Точный ответ на деление на 5 всегда будет иметь дробную часть 0, 1/5, 2/5, 3/5 или 4/5. Поэтому положительная погрешность менее 1/5 в умноженном и сдвинутом результате никогда не переместит результат за границу округления.

Ошибка в нашей константе составляет (1/5) * 2 ^-64 . Значение i составляет менее 2 ^64, поэтому ошибка после умножения составляет менее 1/5. После деления на 4 ошибка меньше (1/5) * 2 ⁻² .

(1/5) * 2 ⁻² <1/5, поэтому ответ всегда будет равен точному делению и округлению до нуля.

К сожалению, это не работает для всех делителей.

Если мы попытаемся представить 4/7 как число с фиксированной точкой 0,64 с округлением от нуля, мы получим ошибку (6/7) * 2 ^-64 . После умножения на значение i чуть менее 2 ⁶⁴ мы получим ошибку чуть меньше 6/7, а после деления на четыре мы получим ошибку чуть меньше 1,5 / 7, которая больше 1/7.

Таким образом, чтобы правильно реализовать деление на 7, нам нужно умножить на число с фиксированной точкой 0,65. Мы можем реализовать это путем умножения на младшие 64 бита нашего числа с фиксированной запятой, затем добавления исходного числа (это может переполниться в бит переноса) и последующего поворота через перенос.

Вот ссылка на документ алгоритма, который создает значения и код, который я вижу в Visual Studio (в большинстве случаев), и который, как я полагаю, все еще используется в GCC для деления целого числа переменной на целое число константы.

http://gmplib.org/~tege/divcnst-pldi94.pdf

В этой статье uword имеет N битов, udword имеет 2N битов, n = числитель = дивиденд, d = знаменатель = делитель, initially изначально установлен в ceil (log2 (d)), shpre является предварительным сдвигом (используется перед умножением ) = e = количество завершающих нулевых битов в d, shpost - пост-сдвиг (используется после умножения), prec - точность = N - e = N - shpre. Цель состоит в том, чтобы оптимизировать расчет н / д с использованием до сдвига, умножения и постсдвига.

Прокрутите вниз до рисунка 6.2, который определяет, как генерируется множитель udword (максимальный размер N + 1 бит), но не дает четкого объяснения процесса. Я объясню это ниже.

Рисунок 4.2 и рисунок 6.2 показывают, как множитель может быть уменьшен до множителя N бит или меньше для большинства делителей. Уравнение 4.5 объясняет, как была получена формула, используемая для работы с N + 1 битовыми умножителями на рисунках 4.1 и 4.2.

В случае современных X86 и других процессоров время умножения является фиксированным, поэтому предварительное смещение не помогает этим процессорам, но все же помогает уменьшить множитель с N + 1 бит до N бит. Я не знаю, исключили ли GCC или Visual Studio предварительный сдвиг для целей X86.

Возвращаясь к рисунку 6.2. Числитель (дивиденд) для mlow и mhigh может быть больше, чем вымышленное слово, только когда знаменатель (делитель)> 2 ^ (N-1) (когда ℓ == N => mlow = 2 ^ (2N)), в этом случае Оптимизированная замена для n / d - это сравнение (если n> = d, q = 1, иначе q = 0), поэтому множитель не генерируется. Начальные значения mlow и mhigh будут составлять N + 1 бит, и для получения каждого значения N + 1 бита (mlow или mhigh) можно использовать два деления udword / uword. Используя X86 в 64-битном режиме в качестве примера:

; upper 8 bytes of dividend = 2^(ℓ) = (upper part of 2^(N+ℓ))
; lower 8 bytes of dividend for mlow  = 0
; lower 8 bytes of dividend for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e)
dividend  dq    2 dup(?)        ;16 byte dividend
divisor   dq    1 dup(?)        ; 8 byte divisor

; ...
        mov     rcx,divisor
        mov     rdx,0
        mov     rax,dividend+8     ;upper 8 bytes of dividend
        div     rcx                ;after div, rax == 1
        mov     rax,dividend       ;lower 8 bytes of dividend
        div     rcx
        mov     rdx,1              ;rdx:rax = N+1 bit value = 65 bit value

Вы можете проверить это с GCC. Вы уже видели, как обрабатывается j = i / 5. Посмотрите, как обрабатывается j = i / 7 (это должен быть случай умножения N + 1).

На большинстве современных процессоров умножение имеет фиксированную синхронизацию, поэтому предварительная смена не требуется. Для X86 конечный результат представляет собой последовательность из двух команд для большинства делителей и последовательность из пяти команд для делителей, таких как 7 (для того, чтобы эмулировать множитель N + 1 бита, как показано в уравнении 4.5 и на рисунке 4.2 файла PDF). Пример кода X86-64:

;       rax = dividend, rbx = 64 bit (or less) multiplier, rcx = post shift count
;       two instruction sequence for most divisors:

        mul     rbx                     ;rdx = upper 64 bits of product
        shr     rdx,cl                  ;rdx = quotient
;
;       five instruction sequence for divisors like 7
;       to emulate 65 bit multiplier (rbx = lower 64 bits of multiplier)

        mul     rbx                     ;rdx = upper 64 bits of product
        sub     rbx,rdx                 ;rbx -= rdx
        shr     rbx,1                   ;rbx >>= 1
        add     rdx,rbx                 ;rdx = upper 64 bits of corrected product
        shr     rdx,cl                  ;rdx = quotient
;       ...

Я отвечу под немного другим углом: потому что это разрешено делать.

C и C ++ определены против абстрактной машины. Компилятор преобразует эту программу в терминах абстрактной машины в конкретную машину, следуя правилу « как будто» .

• Компилятору разрешено вносить ЛЮБЫЕ изменения, если он не изменяет наблюдаемое поведение, заданное абстрактной машиной. Нет никаких оснований ожидать, что компилятор преобразует ваш код самым простым способом (даже когда многие программисты на Си предполагают это). Обычно это происходит потому, что компилятор хочет оптимизировать производительность по сравнению с простым подходом (как подробно обсуждалось в других ответах).
• Если при каких-либо обстоятельствах компилятор «оптимизирует» правильную программу для чего-то, что имеет другое наблюдаемое поведение, это ошибка компилятора.
• Любое неопределенное поведение в нашем коде (классическое переполнение со знаком является классическим примером), и этот контракт является недействительным.