Учитывая, что целочисленное деление является относительно медленным процессом по сравнению с другим типом целочисленной арифметики, я бы хотел его оптимизировать. Есть ли какой-то альтернативный формат, в котором я мог бы хранить d, чтобы процесс разделения выполнялся быстрее? Может быть, обратная какая-то форма?
Ответ 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;
}
Ответ 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
