Самый быстрый способ вычислить 128-разрядное целое число по модулю 64-разрядное целое число
У меня есть 128-разрядное целое число без знака A и 64-разрядное целое число без знака B. Что состоит в том, чтобы вычислить самый быстрый путь A % B - это - (64-разрядный) остаток от деления A на B?
Я надеюсь делать это или в C или в ассемблере, но я должен быть нацелен на 32-разрядную x86 платформу. Это, к сожалению, означает, что я не могу использовать в своих интересах поддержку компилятора 128-разрядных целых чисел, ни способности x64 архитектуры выполнить необходимую операцию в единственной инструкции.
Править:
Спасибо за ответы до сих пор. Однако кажется мне, что предложенные алгоритмы были бы довольно медленными - не будет самый быстрый способ выполнить 128-разрядное 64-разрядным подразделением для усиления собственной поддержки процессора 64-разрядного 32-разрядным подразделением? Кто-либо знает, существует ли способ выполнить более крупное подразделение с точки зрения нескольких меньших подразделений?
Ре: Как часто B изменяется?
Прежде всего, я интересуюсь общим решением - какое вычисление Вы выполнили бы, если A и B, вероятно, будут отличаться каждый раз?
Однако вторая возможная ситуация состоит в том, что B не варьируется так же часто как - могут быть целых 200 для деления на каждый B. Как Ваш ответ отличался бы в этом случае?
10 ответов
Вы можете использовать версию деления Русского крестьянского умножения.
Чтобы найти остаток, выполните (в псевдокоде):
X = B;
while (X <= A/2)
{
X <<= 1;
}
while (A >= B)
{
if (A >= X)
A -= X;
X >>= 1;
}
Модуль остался в A.
Вам потребуется реализовать сдвиги, сравнения и вычитания для работы со значениями, состоящими из пары 64-битных чисел, но это довольно тривиально (скорее всего, вы должны реализовать сдвиг влево на 1 как X + X).
Это зациклится не более 255 раз (при 128-битном A). Конечно, нужно сделать предварительную проверку на нулевой делитель.
Учитывая A = AH * 2 ^ 64 + AL :
A % B == (((AH % B) * (2^64 % B)) + (AL % B)) % B
== (((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
Если ваш компилятор поддерживает 64-битные целые числа, это, вероятно, самый простой способ.
Реализация MSVC 64-битного модуля на 32-битном x86 представляет собой сборку с запутанным циклом ( VC \ crt \ src \ intel \ llrem.asm для смелых), так что я ' Я лично согласен с этим.
Если у вас последняя машина x86, для SSE2 + есть 128-битные регистры. Я никогда не пробовал писать сборку для чего-либо, кроме базового x86, но подозреваю, что есть некоторые руководства.
Это почти непроверенная, частично модифицированная по скорости функция алгоритма Mod128by64 "Русский крестьянин". К сожалению, я пользователь Delphi, поэтому эта функция работает под Delphi. :) Но ассемблер почти такой же, так что ...
function Mod128by64(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
// : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//Divisor = edx:ebp
//Dividend = bh:ebx:edx //We need 64 bits + 1 bit in bh
//Result = esi:edi
//ecx = Loop counter and Dividend index
push ebx //Store registers to stack
push esi
push edi
push ebp
mov ebp, [edx] //Divisor = edx:ebp
mov edx, [edx + 4]
mov ecx, ebp //Div by 0 test
or ecx, edx
jz @DivByZero
xor edi, edi //Clear result
xor esi, esi
//Start of 64 bit division Loop
mov ecx, 15 //Load byte loop shift counter and Dividend index
@SkipShift8Bits: //Small Dividend numbers shift optimisation
cmp [eax + ecx], ch //Zero test
jnz @EndSkipShiftDividend
loop @SkipShift8Bits //Skip 8 bit loop
@EndSkipShiftDividend:
test edx, $FF000000 //Huge Divisor Numbers Shift Optimisation
jz @Shift8Bits //This Divisor is > $00FFFFFF:FFFFFFFF
mov ecx, 8 //Load byte shift counter
mov esi, [eax + 12] //Do fast 56 bit (7 bytes) shift...
shr esi, cl //esi = $00XXXXXX
mov edi, [eax + 9] //Load for one byte right shifted 32 bit value
@Shift8Bits:
mov bl, [eax + ecx] //Load 8 bits of Dividend
//Here we can unrole partial loop 8 bit division to increase execution speed...
mov ch, 8 //Set partial byte counter value
@Do65BitsShift:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
setc bh //Save 65th bit
sub edi, ebp //Compare dividend and divisor
sbb esi, edx //Subtract the divisor
sbb bh, 0 //Use 65th bit in bh
jnc @NoCarryAtCmp //Test...
add edi, ebp //Return privius dividend state
adc esi, edx
@NoCarryAtCmp:
dec ch //Decrement counter
jnz @Do65BitsShift
//End of 8 bit (byte) partial division loop
dec cl //Decrement byte loop shift counter
jns @Shift8Bits //Last jump at cl = 0!!!
//End of 64 bit division loop
mov eax, edi //Load result to eax:edx
mov edx, esi
@RestoreRegisters:
pop ebp //Restore Registers
pop edi
pop esi
pop ebx
ret
@DivByZero:
xor eax, eax //Here you can raise Div by 0 exception, now function only return 0.
xor edx, edx
jmp @RestoreRegisters
end;
Возможна как минимум еще одна оптимизация скорости! После «Оптимизации сдвига огромных чисел делителя» мы можем проверить старший бит делителя, если он равен 0, нам не нужно использовать дополнительный регистр bh в качестве 65-го бита для хранения в нем. Итак, развернутая часть цикла может выглядеть так:
shl bl,1 //Shift dividend left for one bit
rcl edi,1
rcl esi,1
sub edi, ebp //Compare dividend and divisor
sbb esi, edx //Subtract the divisor
jnc @NoCarryAtCmpX
add edi, ebp //Return privius dividend state
adc esi, edx
@NoCarryAtCmpX:
Возможно, вы ищете законченную программу, но основные алгоритмы арифметики с высокой точностью можно найти в книге Кнута Искусство компьютерного программирования , том 2. Вы можете найти алгоритм деления, описанный в Интернете здесь . Алгоритмы имеют дело с произвольной арифметикой с высокой точностью и поэтому являются более общими, чем вам нужно, но вы должны иметь возможность упростить их для 128-битной арифметики, выполняемой с 64- или 32-битными цифрами. Будьте готовы к разумному объему работы (а) по пониманию алгоритма и (б) по его преобразованию на C или ассемблер.
Вы также можете попробовать Hacker's Delight , который полон очень умного ассемблера и другого низкоуровневого хакерства, включая некоторую арифметику с высокой точностью.
Решение зависит от того, что именно вы пытаетесь решить.
Например. если вы выполняете арифметические операции в кольце по модулю 64-битного целого числа, то использование редукции Монтгомери очень эффективно. Конечно, это предполагает, что у вас много раз один и тот же модуль упругости, и что стоит преобразовать элементы кольца в специальное представление.
Чтобы дать очень приблизительную оценку скорости этого сокращения Монтгомери: у меня есть старый тест, который выполняет модульное возведение в степень с 64-битным модулем и экспонентой за 1600 нс на ядре 2 с частотой 2,4 ГГц. Это возведение в степень дает примерно 96 модульное умножение (и модульное сокращение) и, следовательно, требует около 40 циклов на модульное умножение.
Поскольку в C нет предопределенного 128-битного целочисленного типа, биты A должны быть представлены в виде массива. Хотя B (64-битное целое число) может быть сохранено в переменной unsigned long long int , необходимо поместить биты B в другой массив, чтобы эффективно работать с A и B.
После этого B увеличивается на Bx2, Bx3, Bx4, ... до тех пор, пока не станет наибольшим B меньше A. И затем (AB) может быть вычислено, используя некоторые знания вычитания для базы 2.
Is это решение, которое вы ищете?
Как правило, деление происходит медленно, умножение - быстрее, а битовый сдвиг - еще быстрее. Из того, что я видел до сих пор, в большинстве ответов использовался подход грубой силы с использованием битовых сдвигов. Есть другой способ. Будет ли это быстрее, еще неизвестно (AKA профилирует это).
Вместо деления умножьте на обратное. Таким образом, чтобы обнаружить A% B, сначала вычислите обратную величину B ... 1 / B. Это можно сделать с помощью нескольких циклов, используя метод сходимости Ньютона-Рафсона. Это будет зависеть от хорошего набора начальных значений в таблице.
Для получения более подробной информации о методе схождения обратного числа по методу Ньютона-Рафсона, пожалуйста, обратитесь к http://en.wikipedia.org/wiki/Division_ (digital)
. частное Q = A * 1 / B.
Остаток R = A - Q * B.
Чтобы определить, будет ли это быстрее, чем грубая сила (так как будет намного больше умножений, поскольку мы будем использовать 32-битные регистры для имитации 64-битных и 128-битных чисел, профилируйте это.
Если B является константой в вашем коде, вы можете предварительно вычислить обратную величину и просто вычислить, используя последние две формулы. Я уверен, что это будет быстрее, чем битовый сдвиг.
Надеюсь, это поможет.
Я сделал обе версии функции разделения Mod128by64 'Русский крестьянин': классическую и оптимизированную по скорости. Оптимизированная по скорости может выполнять на моем 3Ghz PC более 1000.000 случайных вычислений в секунду и более чем в три раза быстрее, чем классическая функция. Если сравнить время выполнения вычислений 128 на 64 и вычислений 64 на 64 по модулю, то эта функция всего лишь на 50% медленнее.
Классический Русский Крестьянин:
function Mod128by64Clasic(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
// : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//edx:ebp = Divisor
//ecx = Loop counter
//Result = esi:edi
push ebx //Store registers to stack
push esi
push edi
push ebp
mov ebp, [edx] //Load divisor to edx:ebp
mov edx, [edx + 4]
mov ecx, ebp //Div by 0 test
or ecx, edx
jz @DivByZero
push [eax] //Store Divisor to the stack
push [eax + 4]
push [eax + 8]
push [eax + 12]
xor edi, edi //Clear result
xor esi, esi
mov ecx, 128 //Load shift counter
@Do128BitsShift:
shl [esp + 12], 1 //Shift dividend from stack left for one bit
rcl [esp + 8], 1
rcl [esp + 4], 1
rcl [esp], 1
rcl edi, 1
rcl esi, 1
setc bh //Save 65th bit
sub edi, ebp //Compare dividend and divisor
sbb esi, edx //Subtract the divisor
sbb bh, 0 //Use 65th bit in bh
jnc @NoCarryAtCmp //Test...
add edi, ebp //Return privius dividend state
adc esi, edx
@NoCarryAtCmp:
loop @Do128BitsShift
//End of 128 bit division loop
mov eax, edi //Load result to eax:edx
mov edx, esi
@RestoreRegisters:
lea esp, esp + 16 //Restore Divisors space on stack
pop ebp //Restore Registers
pop edi
pop esi
pop ebx
ret
@DivByZero:
xor eax, eax //Here you can raise Div by 0 exception, now function only return 0.
xor edx, edx
jmp @RestoreRegisters
end;
Оптимизированный по скорости Русский Крестьянин:
function Mod128by64Oprimized(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
// : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//Divisor = edx:ebp
//Dividend = ebx:edx //We need 64 bits
//Result = esi:edi
//ecx = Loop counter and Dividend index
push ebx //Store registers to stack
push esi
push edi
push ebp
mov ebp, [edx] //Divisor = edx:ebp
mov edx, [edx + 4]
mov ecx, ebp //Div by 0 test
or ecx, edx
jz @DivByZero
xor edi, edi //Clear result
xor esi, esi
//Start of 64 bit division Loop
mov ecx, 15 //Load byte loop shift counter and Dividend index
@SkipShift8Bits: //Small Dividend numbers shift optimisation
cmp [eax + ecx], ch //Zero test
jnz @EndSkipShiftDividend
loop @SkipShift8Bits //Skip Compute 8 Bits unroled loop ?
@EndSkipShiftDividend:
test edx, $FF000000 //Huge Divisor Numbers Shift Optimisation
jz @Shift8Bits //This Divisor is > $00FFFFFF:FFFFFFFF
mov ecx, 8 //Load byte shift counter
mov esi, [eax + 12] //Do fast 56 bit (7 bytes) shift...
shr esi, cl //esi = $00XXXXXX
mov edi, [eax + 9] //Load for one byte right shifted 32 bit value
@Shift8Bits:
mov bl, [eax + ecx] //Load 8 bit part of Dividend
//Compute 8 Bits unroled loop
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove0 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow0
ja @DividentAbove0
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow0
@DividentAbove0:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow0:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove1 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow1
ja @DividentAbove1
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow1
@DividentAbove1:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow1:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove2 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow2
ja @DividentAbove2
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow2
@DividentAbove2:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow2:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove3 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow3
ja @DividentAbove3
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow3
@DividentAbove3:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow3:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove4 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow4
ja @DividentAbove4
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow4
@DividentAbove4:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow4:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove5 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow5
ja @DividentAbove5
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow5
@DividentAbove5:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow5:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove6 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow6
ja @DividentAbove6
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow6
@DividentAbove6:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow6:
shl bl, 1 //Shift dividend left for one bit
rcl edi, 1
rcl esi, 1
jc @DividentAbove7 //dividend hi bit set?
cmp esi, edx //dividend hi part larger?
jb @DividentBelow7
ja @DividentAbove7
cmp edi, ebp //dividend lo part larger?
jb @DividentBelow7
@DividentAbove7:
sub edi, ebp //Return privius dividend state
sbb esi, edx
@DividentBelow7:
//End of Compute 8 Bits (unroled loop)
dec cl //Decrement byte loop shift counter
jns @Shift8Bits //Last jump at cl = 0!!!
//End of division loop
mov eax, edi //Load result to eax:edx
mov edx, esi
@RestoreRegisters:
pop ebp //Restore Registers
pop edi
pop esi
pop ebx
ret
@DivByZero:
xor eax, eax //Here you can raise Div by 0 exception, now function only return 0.
xor edx, edx
jmp @RestoreRegisters
end;
Я хотел бы поделиться несколькими мыслями.
Боюсь, что все не так просто, как предлагает MSN.
В выражении:
(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
и умножение, и сложение могут переполняться. Я думаю, что это можно учесть и использовать общую концепцию с некоторыми модификациями, но что-то мне подсказывает, что это будет очень страшно.
Мне было интересно, как реализована 64-битная операция модуляции в MSVC, и я попытался что-нибудь выяснить. Я не очень хорошо знаю ассемблер, и все, что у меня было доступно, это Express edition, без исходников VC\crt\src\intel\llrem.asm, но я думаю, что мне удалось получить некоторое представление о том, что происходит, после небольшой игры с отладчиком и выводом дизассемблера. Я пытался понять, как вычисляется остаток в случае положительных целых чисел и делителя >=2^32. Конечно, есть код, который работает с отрицательными числами, но я не стал копаться в этом.
Вот как я это вижу:
Если делитель >= 2^32, то и делимое, и делитель сдвигаются вправо настолько, насколько это необходимо, чтобы вместить делитель в 32 бита. Другими словами: если для записи делителя в двоичном виде требуется n цифр и n > 32, то отбрасываются n-32 младших разряда как делителя, так и делимого. После этого выполняется деление с использованием аппаратной поддержки деления 64-битных целых чисел на 32-битные. Результат может быть неправильным, но я думаю, что можно доказать, что результат может отличаться не более чем на 1. После деления делитель (исходный) умножается на результат, а произведение вычитается из дивиденда. Затем его корректируют, добавляя или вычитая делитель, если это необходимо (если результат деления отклонился на единицу).
Легко разделить 128-битное целое число на 32-битное, используя аппаратную поддержку 64-битного деления на 32-битное. В случае, если делитель < 2^32, можно вычислить остаток, выполнив всего 4 деления следующим образом:
Предположим, что делитель хранится в:
DWORD dividend[4] = ...
остаток будет в:
DWORD remainder;
1) Divide dividend[3] by divisor. Store the remainder in remainder.
2) Divide QWORD (remainder:dividend[2]) by divisor. Store the remainder in remainder.
3) Divide QWORD (remainder:dividend[1]) by divisor. Store the remainder in remainder.
4) Divide QWORD (remainder:dividend[0]) by divisor. Store the remainder in remainder.
После этих 4 шагов переменная остатка будет содержать то, что вы ищете. (Пожалуйста, не убивайте меня, если я неправильно указал порядок следования. Я даже не программист)
В случае, если делитель больше 2^32-1, у меня нет хороших новостей. У меня нет полного доказательства того, что результат после сдвига отклоняется не более чем на 1, в процедуре, которую я описал ранее и которую, как я полагаю, использует MSVC. Однако я думаю, что это связано с тем, что отбрасываемая часть по крайней мере в 2^31 раз меньше делителя, делимое меньше 2^64, а делитель больше 2^32-1, поэтому результат меньше 2^32.
Если делимое имеет 128 бит, то трюк с отбрасыванием битов не сработает. Так что в общем случае лучшим решением, вероятно, является решение, предложенное GJ или caf. (Ну, оно было бы лучшим, даже если бы отбрасывание битов работало. Деление, умножение, вычитание и коррекция на 128-битных целых числах могут быть медленнее.)
Я также думал об использовании аппаратного обеспечения с плавающей запятой. x87 блок с плавающей запятой использует формат 80-битной точности с дробью длиной 64 бита. Я думаю, что можно получить точный результат деления 64 бита на 64 бита. (Не непосредственно остаток, а остаток с использованием умножения и вычитания, как в "процедуре MSVC"). Если делимое >=2^64 и < 2^128, то хранение его в формате с плавающей точкой похоже на отбрасывание младших значащих бит в "процедуре MSVC". Может быть, кто-то сможет доказать, что ошибка в этом случае является пограничной, и сочтет это полезным. Я понятия не имею, есть ли у него шанс оказаться быстрее, чем решение GJ, но, возможно, стоит попробовать.