Dlaczego GCC używa mnożenia przez dziwną liczbę w implementacji dzielenia liczb całkowitych?

Czytałem o operacjach montażu divi mulmontażu i postanowiłem zobaczyć je w akcji, pisząc prosty program w C:

Podział pliku. C

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

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

A następnie generowanie kodu języka asemblera za pomocą:

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

Ale patrząc na wygenerowany division.splik, nie zawiera żadnych operacji div! Zamiast tego robi jakąś czarną magię z przesunięciem bitów i magicznymi liczbami. Oto fragment kodu, który oblicza 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

Co tu się dzieje? Dlaczego GCC w ogóle nie używa div? Jak generuje tę magiczną liczbę i dlaczego wszystko działa?

Podział liczb całkowitych jest jedną z najwolniejszych operacji arytmetycznych, które można wykonać na nowoczesnym procesorze, z opóźnieniem do kilkudziesięciu cykli i złą przepustowością. (Dla x86, zobacz tabele instrukcji Agner Fog i przewodnik mikroarchitektury ).

Jeśli znasz dzielnik z wyprzedzeniem, możesz uniknąć podziału, zastępując go zestawem innych operacji (zwielokrotnienia, uzupełnienia i przesunięcia), które mają równoważny efekt. Nawet jeśli potrzebnych jest kilka operacji, często jest to o wiele szybsze niż samo dzielenie liczb całkowitych.

Zaimplementowanie /operatora C w ten sposób zamiast sekwencji obejmującej wiele instrukcji divjest tylko domyślnym sposobem GCC dzielenia według stałych. Nie wymaga optymalizacji między operacjami i niczego nie zmienia nawet podczas debugowania. (Używanie w -Osprzypadku małych rozmiarów kodu powoduje jednak użycie GCC div.) Używanie odwrotności multiplikatywnej zamiast dzielenia jest jak używanie leazamiast muliadd

W wyniku tego masz tendencję do wyświetlania divlub idivwyświetlania danych wyjściowych tylko wtedy, gdy dzielnik nie jest znany w czasie kompilacji.

Aby uzyskać informacje o tym, jak kompilator generuje te sekwencje, a także kod umożliwiający ich wygenerowanie dla siebie (prawie na pewno niepotrzebne, chyba że pracujesz z kompilatorem braindead), zobacz libdivide .

Dzielenie przez 5 jest takie samo jak mnożenie 1/5, co znowu jest takie samo jak mnożenie przez 4/5 i przesunięcie w prawo o 2 bity. Dana wartość jest CCCCCCCCCCCCCCCDw postaci szesnastkowej, która jest reprezentacją binarną 4/5, jeśli jest umieszczona po punkcie szesnastkowym (tzn. 0.110011001100Powtarza się wartość binarna dla czterech piątych - dlaczego poniżej). Myślę, że możesz wziąć to stąd! Możesz chcieć sprawdzić arytmetykę stałoprzecinkową (choć pamiętaj, że na końcu jest ona zaokrąglana do liczby całkowitej).

Dlaczego mnożenie jest szybsze niż dzielenie, a gdy dzielnik jest stały, jest to szybsza trasa.

Zobacz samouczek Wzajemne mnożenie, w którym znajduje się szczegółowy opis działania, wyjaśniający w kategoriach punktu stałego. Pokazuje, jak działa algorytm znajdowania odwrotności i jak radzić sobie z podpisanym podziałem i modulo.

Zastanówmy się przez chwilę, dlaczego 0.CCCCCCCC...(hex) lub 0.110011001100...binarny to 4/5. Podziel reprezentację binarną przez 4 (przesuń w prawo o 2 miejsca), a my otrzymamy, 0.001100110011...które poprzez trywialną inspekcję można dodać oryginał 0.111111111111..., który jest oczywiście równy 1, w ten sam sposób 0.9999999...w systemie dziesiętnym jest równy jeden. Dlatego wiemy, że x + x/4 = 1tak 5x/4 = 1, x=4/5. Jest to następnie przedstawiane jako CCCCCCCCCCCCDszesnastkowe dla zaokrąglania (ponieważ cyfra binarna poza ostatnią obecną byłaby a 1).

Ogólnie rzecz biorąc, mnożenie jest znacznie szybsze niż dzielenie. Jeśli więc uda nam się uniknąć mnożenia przez odwrotność, możemy znacznie przyspieszyć dzielenie o stałą

Zmarszczka polega na tym, że nie możemy dokładnie przedstawić odwrotności (chyba że podział był potęgą dwóch, ale w takim przypadku zwykle możemy po prostu przekształcić podział na odrobinę przesunięcia). Dlatego, aby zapewnić prawidłowe odpowiedzi, musimy uważać, aby błąd w naszej wzajemności nie powodował błędów w naszym wyniku końcowym.

-3689348814741910323 to 0xCCCCCCCCCCCCCCCCCD, która jest wartością nieco ponad 4/5 wyrażoną w stałym punkcie 0,64.

Kiedy pomnożymy 64-bitową liczbę całkowitą przez stałą liczbę 0,64, otrzymamy wynik 64,64. Obcinamy wartość do 64-bitowej liczby całkowitej (skutecznie zaokrąglając ją do zera), a następnie wykonujemy dalsze przesunięcie, które dzieli się przez cztery i ponownie obcinamy. Patrząc na poziom bitów, jasne jest, że możemy traktować oba skróty jako pojedyncze obcięcie.

To wyraźnie daje nam przynajmniej przybliżenie podziału przez 5, ale czy daje nam dokładną odpowiedź poprawnie zaokrągloną do zera?

Aby uzyskać dokładną odpowiedź, błąd musi być wystarczająco mały, aby nie przesuwać odpowiedzi poza zaokrągloną granicę.

Dokładna odpowiedź na podział przez 5 zawsze będzie miała ułamkową część 0, 1/5, 2/5, 3/5 lub 4/5. Dlatego dodatni błąd mniejszy niż 1/5 w pomnożonym i przesuniętym wyniku nigdy nie przesunie wyniku poza zaokrągloną granicę.

Błąd w naszej stałej wynosi (1/5) * ^2-64 . Wartość i jest mniejsza niż 2 ^64, więc błąd po pomnożeniu jest mniejszy niż 1/5. Po podzieleniu przez 4 błąd jest mniejszy niż (1/5) * 2 ⁻² .

(1/5) * 2 ⁻² <1/5, więc odpowiedź zawsze będzie równa dokładnemu podziałowi i zaokrągleniu do zera.

Niestety nie działa to dla wszystkich dzielników.

Jeśli spróbujemy przedstawić 4/7 jako stałą liczbę 0,64 z zaokrągleniem od zera, otrzymamy błąd (6/7) * ^2-64 . Po pomnożeniu przez wartość i nieco poniżej 2 ⁶⁴ otrzymujemy błąd poniżej 6/7, a po podzieleniu przez cztery otrzymujemy błąd nieco poniżej 1,5 / 7, który jest większy niż 1/7.

Tak więc, aby poprawnie wdrożyć podział przez 7, musimy pomnożyć przez stałą liczbę 0,65. Możemy to zaimplementować, mnożąc przez dolne 64 bity naszego stałego numeru punktu, a następnie dodając pierwotną liczbę (może to przelać się do bitu przenoszenia), a następnie wykonując obrót przez przeniesienie.

Oto link do dokumentu algorytmu, który generuje wartości i kod, które widzę w Visual Studio (w większości przypadków) i który, jak zakładam, jest nadal używany w GCC do dzielenia zmiennej całkowitej przez stałą liczbę całkowitą.

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

W artykule uword ma N bitów, słowo ma 2N bitów, n = licznik = dywidenda, d = mianownik = dzielnik, ℓ jest początkowo ustawiony na pułap (log2 (d)), shpre jest przed przesunięciem (używane przed pomnożeniem ) = e = liczba ostatnich zerowych bitów w d, shpost jest przesunięciem (używane po pomnożeniu), prec jest precyzją = N - e = N - shpre. Celem jest optymalizacja obliczeń n / d przy użyciu zmiany wstępnej, pomnożenia i zmiany następczej.

Przewiń w dół do rysunku 6.2, który definiuje sposób generowania mnożnika udword (maksymalny rozmiar to N + 1 bit), ale nie wyjaśnia dokładnie procesu. Wyjaśnię to poniżej.

Rysunek 4.2 i rysunek 6.2 pokazują, w jaki sposób mnożnik można zredukować do mnożnika N bitowego lub mniejszego dla większości dzielników. Równanie 4.5 wyjaśnia, w jaki sposób wyprowadzono wzór używany do radzenia sobie z mnożnikami bitów N + 1 na rysunkach 4.1 i 4.2.

W przypadku współczesnych X86 i innych procesorów czas mnożenia jest stały, więc wstępne przesunięcie nie pomaga w tych procesorach, ale nadal pomaga zmniejszyć mnożnik z N + 1 bitów do N bitów. Nie wiem, czy GCC czy Visual Studio wyeliminowały wstępne zmiany dla celów X86.

Wracając do rysunku 6.2. Licznik (dywidenda) dla mlow i mhigh może być większy niż słowo ud, tylko gdy mianownik (dzielnik)> 2 ^ (N-1) (gdy ℓ == N => mlow = 2 ^ (2N)), w tym przypadku zoptymalizowane zastąpienie dla n / d jest porównaniem (jeśli n> = d, q = 1, w innym przypadku q = 0), więc nie jest generowany mnożnik. Początkowe wartości mlow i mhigh będą wynosić N + 1 bitów, a do podziału każdej wartości bitowej N + 1 (mlow lub mhigh) można użyć dwóch podziałów udword / uword. Używanie X86 w trybie 64-bitowym jako przykład:

; 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

Możesz to przetestować za pomocą GCC. Już wiesz, jak obsługiwane jest j = i / 5. Zobacz, jak obsługiwane jest j = i / 7 (co powinno być przypadkiem mnożnika N + 1 bit).

W większości obecnych procesorów mnożenie ma ustalony czas, więc zmiana wstępna nie jest potrzebna. W przypadku X86 wynikiem końcowym jest sekwencja dwóch instrukcji dla większości dzielników i pięć sekwencji instrukcji dla dzielników takich jak 7 (w celu emulacji mnożnika bitów N + 1, jak pokazano w równaniu 4.5 i rysunku 4.2 pliku pdf). Przykładowy kod 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
;       ...

Odpowiem z nieco innej strony: ponieważ wolno to zrobić.

C i C ++ są zdefiniowane na maszynie abstrakcyjnej. Kompilator przekształca ten program w kategoriach abstrakcyjnej maszyny do betonu maszyny następujących po co-jeśli reguły.

• Kompilator może dokonywać DOWOLNYCH zmian, o ile nie zmienia obserwowalnego zachowania określonego przez maszynę abstrakcyjną. Nie ma uzasadnionych oczekiwań, że kompilator przekształci Twój kod w najprostszy możliwy możliwy sposób (nawet jeśli wielu programistów C zakłada to). Zwykle robi to, ponieważ kompilator chce zoptymalizować wydajność w porównaniu z prostym podejściem (jak omówiono szczegółowo w innych odpowiedziach).
• Jeśli w jakichkolwiek okolicznościach kompilator „zoptymalizuje” poprawny program do czegoś, co ma inne obserwowalne zachowanie, jest to błąd kompilatora.
• Wszelkie nieokreślone zachowanie w naszym kodzie (klasyczny przepełnienie z podpisem jest klasycznym przykładem), a niniejsza umowa jest nieważna.