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

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

File division.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 вообще? Как он генерирует это магическое число и почему все работает?

Ответы

Ответ 1

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

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

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

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

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

Ответ 2

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

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

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

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

Ответ 3

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

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

-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. Значение я составляет менее 2 64, поэтому ошибка после умножения составляет менее 1/5. После деления на 4 погрешность составляет менее (1/5) * 2 −2.

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

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

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

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

Ответ 4

Здесь приведена ссылка на документ алгоритма, который создает значения и код, которые я вижу в 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
;       ...