Повторное целочисленное деление на постоянное значение времени выполнения

Ответ 1

Для этого есть библиотека libdivide:

libdivide - это библиотека с открытым исходным кодом для оптимизации целочисленного деления

libdivide позволяет вам заменить дорогостоящие целые деления на сравнительно дешевое умножение и бит-сдвиги. Компиляторы обычно делают это, но только тогда, когда делитель известен во время компиляции. libdivide позволяет вам использовать его во время выполнения. В результате целочисленное деление может стать быстрее - намного быстрее. Более того, libdivide позволяет разделить вектор SSE2 на константу времени выполнения, что особенно хорошо, потому что SSE2 не имеет целочисленного деления инструкция!

libdivide является бесплатным и открытым исходным кодом с разрешительной лицензией. Имя "libdivide" - это немного шутка, поскольку нет библиотеки как таковой: код полностью упакован как один заголовочный файл, причем как C, так и С++ API.

Вы можете прочитать об этом алгоритме на blog; например, запись.

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

Примечание. Вы можете создать еще более быструю версию libdivide. Идея состоит в том, что для каждого делителя вы всегда можете создать триплет (mul/add/shift), поэтому это выражение дает результат: (num * mul + add) → shift (здесь множится умножение). Интересно, что этот метод может превзойти версию компилятора для постоянного деления на несколько микрообъектов!

Здесь моя реализация (это не компилируется из коробки, но можно увидеть общий алгоритм):

struct Divider_u32 {
    u32 mul;
    u32 add;
    s32 shift;

    void set(u32 divider);
};

void Divider_u32::set(u32 divider) {
    s32 l = indexOfMostSignificantBit(divider);
    if (divider&(divider-1)) {
        u64 m = static_cast<u64>(1)<<(l+32);
        mul = static_cast<u32>(m/divider);

        u32 rem = static_cast<u32>(m)-mul*divider;
        u32 e = divider-rem;

        if (e<static_cast<u32>(1)<<l) {
            mul++;
            add = 0;
        } else {
            add = mul;
        }
        shift = l;
    } else {
        if (divider==1) {
            mul = 0xffffffff;
            add = 0xffffffff;
            shift = 0;
        } else {
            mul = 0x80000000;
            add = 0;
            shift = l-1;
        }
    }
}

u32 operator/(u32 v, const Divider_u32 &div) {
    u32 t = static_cast<u32>((static_cast<u64>(v)*div.mul+div.add)>>32)>>div.shift;

    return t;
}

Ответ 2

В книге "Хакерский восторг" есть "Глава 10: Целочисленное деление на константу", охватывающее 74 страницы. Вы можете найти все примеры кода бесплатно в этом каталоге: http://www.hackersdelight.org/hdcode.htm В вашем случае на рис. 10-1., 10-2 и 10-3 - это то, что вы хотите.

Проблема деления на константу d эквивалентна mutiplying на c = 1/d. Эти алгоритмы вычисляют такую ​​константу для вас. Когда у вас есть c, вы вычисляете дивиденд как таковой:

int divideByMyConstant(int dividend){
  int c = MAGIC; // Given by the algorithm

  // since 1/d < 1, c is actually (1<<k)/d to fit nicely ina 32 bit int
  int k = MAGIC_SHIFT; //Also given by the algorithm

  long long tmp = (long long)dividend * c; // use 64 bit to hold all the precision...

  tmp >>= k; // Manual floating point number =)

  return (int)tmp;
}

Ответ 3

update - в моем первоначальном ответе я заметил алгоритм, упомянутый в предыдущем потоке для генерируемого кодом компилятора для деления на константу. Код сборки был написан в соответствии с документом, связанным с предыдущим потоком. Сгенерированный код компилятора включает несколько разных последовательностей в зависимости от делителя.

В этой ситуации делитель неизвестен до времени выполнения, поэтому необходим общий алгоритм. Пример в ответе geza показывает общий алгоритм, который может быть встроен в код сборки с GCC, но Visual Studio не поддерживает встроенную сборку в режиме 64 бит. В случае с Visual Studio существует компромисс между дополнительным кодом, связанным с использованием встроенных функций, и вызовом функции, записанной в сборке. В моей системе (Intel 3770k 3.5ghz) я попытался вызвать одну функцию, которая делает | mul add adc shr |, и я также попытался использовать указатель для функции, чтобы, возможно, использовать более короткие последовательности | mul shr | или | shr (1) mul shr | в зависимости от делителя, но это обеспечило мало или вообще не получило усиления, в зависимости от компилятора. Основными накладными расходами в этом случае является вызов (против | mul add adc shr |). Даже с накладными вызовами, последовательность | call mul add adc shr ret | в среднем примерно в 4 раза быстрее, чем на моей системе.

Обратите внимание, что связанный с исходным кодом для libdivide в ответе geza не имеет общей подпрограммы, которая может обрабатывать divisor == 1. Обычная последовательность libdivide является умножением, вычитанием, сдвигом (1), добавлением, сдвигом, против geza пример С++ последовательность умножения, добавления, adc, shift.

Мой первоначальный ответ: приведенный ниже пример кода использует алгоритм, описанный в предыдущем потоке.

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

Это ссылка на документ, упомянутый в другом потоке:

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

Нижеприведенный пример основан на документе pdf и предназначен для Visual Studio, используя ml64 (64-разрядный ассемблер), работающий в Windows (64-разрядная ОС). Код с метками main... и dcm... - это код для генерации преспирания (rspre, количество конечных нулевых бит в делителе), множитель и пост-сдвиг (rspost). Код с метками dct... является кодом для проверки метода.

        includelib      msvcrtd
        includelib      oldnames

sw      equ     8                       ;size of word

        .data
arrd    dq      1                       ;array of test divisors
        dq      2
        dq      3
        dq      4
        dq      5
        dq      7
        dq      256
        dq      3*256
        dq      7*256
        dq      67116375
        dq      07fffffffffffffffh      ;max divisor
        dq      0
        .data?

        .code
        PUBLIC  main

main    PROC
        push    rbp
        push    rdi
        push    rsi
        sub     rsp,64                  ;allocate stack space
        mov     rbp,rsp
        lea     rsi,arrd                ;set ptr to array of divisors
        mov     [rbp+6*sw],rsi
        jmp     main1

main0:  mov     [rbp+0*sw],rbx          ;[rbp+0*sw] = rbx = divisor = d
        cmp     rbx,1                   ;if d <= 1, q=n or overflow
        jbe     main1
        bsf     rcx,rbx                 ;rcx = rspre
        mov     [rbp+1*sw],rcx          ;[rbp+1*sw] = rspre
        shr     rbx,cl                  ;rbx = d>>rsc
        bsr     rcx,rbx                 ;rcx = floor(log2(rbx))
        mov     rsi,1                   ;rsi = 2^floor(log2(rbx))
        shl     rsi,cl
        cmp     rsi,rbx                 ;br if power of 2
        je      dcm03
        inc     rcx                     ;rcx = ceil(log2(rcx)) = L = rspost
        shl     rsi,1                   ;rsi = 2^L
;       jz      main1                   ;d > 08000000000000000h, use compare
        add     rcx,[rbp+1*sw]          ;rcx = L+rspre
        cmp     rcx,8*sw                ;if d > 08000000000000000h, use compare
        jae     main1
        mov     rax,1                   ;[rbp+3*sw] = 2^(L+rspre)
        shl     rax,cl
        mov     [rbp+3*sw],rax
        sub     rcx,[rbp+1*sw]          ;rcx = L
        xor     rdx,rdx
        mov     rax,rsi                 ;hi N bits of 2^(N+L)
        div     rbx                     ;rax == 1
        xor     rax,rax                 ;lo N bits of 2^(N+L)
        div     rbx
        mov     rdi,rax                 ;rdi = mlo % 2^N
        xor     rdx,rdx
        mov     rax,rsi                 ;hi N bits of 2^(N+L) + 2^(L+rspre)
        div     rbx                     ;rax == 1
        mov     rax,[rbp+3*sw]          ;lo N bits of 2^(N+L) + 2^(L+rspre)
        div     rbx                     ;rax = mhi % 2^N
        mov     rdx,rdi                 ;rdx = mlo % 2^N
        mov     rbx,8*sw                ;rbx = e = # bits in word
dcm00:  mov     rsi,rdx                 ;rsi = mlo/2
        shr     rsi,1
        mov     rdi,rax                 ;rdi = mhi/2
        shr     rdi,1
        cmp     rsi,rdi                 ;break if mlo >= mhi
        jae     short dcm01
        mov     rdx,rsi                 ;rdx = mlo/2
        mov     rax,rdi                 ;rax = mhi/2
        dec     rbx                     ;e -= 1
        loop    dcm00                   ;loop if --shpost != 0
dcm01:  mov     [rbp+2*sw],rcx          ;[rbp+2*sw] = shpost
        cmp     rbx,8*sw                ;br if N+1 bit multiplier
        je      short dcm02
        xor     rdx,rdx
        mov     rdi,1                   ;rax = m = mhi + 2^e
        mov     rcx,rbx
        shl     rdi,cl
        or      rax,rdi
        jmp     short dct00

dcm02:  mov     rdx,1                   ;rdx = 2^N
        dec     qword ptr [rbp+2*sw]    ;dec rspost
        jmp     short dct00

dcm03:  mov     rcx,[rbp+1*sw]          ;rcx = rsc
        jmp     short dct10

;       test    rbx = n, rdx = m bit N, rax = m%(2^N)
;               [rbp+1*sw] = rspre, [rbp+2*sw] = rspost

dct00:  mov     rdi,rdx                 ;rdi:rsi = m
        mov     rsi,rax
        mov     rbx,0fffffffff0000000h  ;[rbp+5*sw] = rbx = n
dct01:  mov     [rbp+5*sw],rbx
        mov     rdx,rdi                 ;rdx:rax = m
        mov     rax,rsi
        mov     rcx,[rbp+1*sw]          ;rbx = n>>rspre
        shr     rbx,cl
        or      rdx,rdx                 ;br if 65 bit m
        jnz     short dct02
        mul     rbx                     ;rdx = (n*m)>>N
        jmp     short dct03

dct02:  mul     rbx
        sub     rbx,rdx
        shr     rbx,1
        add     rdx,rbx
dct03:  mov     rcx,[rbp+2*sw]          ;rcx = rspost
        shr     rdx,cl                  ;rdx = q = quotient
        mov     [rbp+4*sw],rdx          ;[rbp+4*sw] = q
        xor     rdx,rdx                 ;rdx:rax = n
        mov     rax,[rbp+5*sw]
        mov     rbx,[rbp+0*sw]          ;rbx = d
        div     rbx                     ;rax = n/d
        mov     rdx,[rbp+4*sw]          ;br if ok
        cmp     rax,rdx                 ;br if ok
        je      short dct04
        nop                             ;debug check
dct04:  mov     rbx,[rbp+5*sw]
        inc     rbx
        jnz     short dct01
        jmp     short main1

;       test    rbx = n, rcx = rsc
dct10:  mov     rbx,0fffffffff0000000h  ;rbx = n
dct11:  mov     rsi,rbx                 ;rsi = n
        shr     rsi,cl                  ;rsi = n>>rsc
        xor     edx,edx
        mov     rax,rbx
        mov     rdi,[rbp+0*sw]          ;rdi = d
        div     rdi
        cmp     rax,rsi                 ;br if ok
        je      short dct12
        nop
dct12:  inc     rbx
        jnz     short dct11

main1:  mov     rsi,[rbp+6*sw]          ;rsi ptr to divisor
        mov     rbx,[rsi]               ;rbx = divisor = d
        add     rsi,1*sw                ;advance ptr
        mov     [rbp+6*sw],rsi
        or      rbx,rbx
        jnz     main0                   ;br if not end table

        add     rsp,64                  ;restore regs
        pop     rsi
        pop     rdi
        pop     rbp
        xor     rax,rax
        ret     0

main    ENDP
        END