128비트 정수 모듈로 64비트 정수를 계산하는 가장 빠른 방법
128비트 부호 없는 정수 A와 64비트 부호 없는 정수 B가 있습니다. 계산하는 가장 빠른 방법은 무엇입니까 A % B? A를 B로 나눈 나머지(64비트)는 무엇입니까?
C 또는 어셈블리 언어로 이 작업을 수행하려고 하지만 32비트 x86 플랫폼을 대상으로 해야 합니다. 이것은 불행히도 128비트 정수에 대한 컴파일러 지원이나 단일 명령어에서 필요한 작업을 수행하는 x64 아키텍처의 기능을 활용할 수 없음을 의미합니다.
편집하다:
지금까지 답변에 감사드립니다. 그러나 제안된 알고리즘은 매우 느릴 것 같습니다. 128비트 x 64비트 분할을 수행하는 가장 빠른 방법은 64비트 x 32비트 분할에 대한 프로세서의 기본 지원을 활용하는 것입니다. 몇 개의 작은 부문으로 더 큰 부문을 수행하는 방법이 있는지 아는 사람이 있습니까?
Re: B는 얼마나 자주 바뀌나요?
주로 일반적인 솔루션에 관심이 있습니다. A와 B가 매번 다를 가능성이 있는 경우 어떤 계산을 수행하시겠습니까?
그러나 두 번째 가능한 상황은 B가 A만큼 자주 변하지 않는다는 것입니다. 각 B로 나누는 것은 200만큼 많을 수 있습니다. 이 경우 답은 어떻게 달라지나요?
러시아 농민 곱셈 의 나눗셈 버전을 사용할 수 있습니다 .
나머지를 찾으려면 (의사 코드에서) 다음을 실행하십시오.
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 사용). 물론 0의 제수를 미리 확인해야 합니다.
완성된 프로그램을 찾고 있을 수도 있지만 다중 정밀도 산술을 위한 기본 알고리즘은 Knuth의 Art of Computer Programming , Volume 2에서 찾을 수 있습니다 . 여기에서 온라인으로 설명된 나눗셈 알고리즘을 찾을 수 있습니다 . 알고리즘은 임의의 다중 정밀도 산술을 다루므로 필요한 것보다 더 일반적이지만 64비트 또는 32비트 숫자에서 수행되는 128비트 산술에 대해 알고리즘을 단순화할 수 있어야 합니다. (a) 알고리즘을 이해하고, (b) 그것을 C 또는 어셈블러로 변환하는 작업에 대해 합리적인 양의 준비를 하십시오.
또한 매우 영리한 어셈블러와 다중 정밀도 산술을 포함하여 기타 저수준 해커로 가득 찬 Hacker's Delight 를 확인하고 싶을 수도 있습니다 .
주어진 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비트 정수를 지원한다면 이것이 아마도 가장 쉬운 방법일 것입니다. 32비트 x86에서 MSVC의 64비트 모듈로 구현은 털이 많은 루프로 채워진 어셈블리( VC\crt\src\intel\llrem.asm용감한 경우)이므로 개인적으로 함께 갈 것입니다.
이것은 거의 테스트되지 않은 부분적으로 수정된 Mod128by64 'Russian peasant' 알고리즘 기능입니다. 불행히도 저는 델파이 사용자이므로 이 기능은 델파이에서 작동합니다. :) 그러나 어셈블러는 거의 동일하므로 ...
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;
최소 한 번 더 속도 최적화가 가능합니다! 'Huge Divisor Numbers Shift Optimization' 후 제수 상위 비트를 테스트할 수 있습니다. 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:
Mod128by64 'Russian peasant' 분할 기능의 두 버전을 모두 만들었습니다. 클래식 및 속도 최적화. 속도 최적화는 3Ghz PC에서 초당 1000.000개 이상의 임의 계산을 수행할 수 있으며 기존 기능보다 3배 이상 빠릅니다. 128 x 64 계산과 64 x 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
곱셈과 덧셈이 모두 오버플로될 수 있습니다. 나는 그것을 고려하고 여전히 약간의 수정과 함께 일반적인 개념을 사용할 수 있다고 생각하지만, 뭔가가 그것이 정말로 무섭게 될 것이라고 말합니다.
MSVC에서 64비트 모듈로 연산이 어떻게 구현되었는지 궁금해서 찾아봤습니다. 나는 어셈블리를 잘 모르고 VC\crt\src\intel\llrem.asm 소스가 없는 Express 에디션만 사용할 수 있었지만 약간의 플레이 후에 무슨 일이 일어나고 있는지 알 수 있었던 것 같습니다. 디버거 및 디스어셈블리 출력과 함께. 양의 정수와 제수가 >=2^32인 경우 나머지가 어떻게 계산되는지 알아내려고 했습니다. 물론 음수를 처리하는 코드가 있지만 자세히 살펴보지는 않았습니다.
내가 보는 방법은 다음과 같습니다.
제수 >= 2^32인 경우 제수를 32비트에 맞추는 데 필요한 만큼 피제수와 제수가 오른쪽으로 이동합니다. 다시 말해, 제수를 이진법으로 기록하는 데 n 자리가 필요하고 n > 32인 경우 제수와 피제수 둘 다의 n-32 최하위 자리는 버려집니다. 그 후, 64비트 정수를 32비트 정수로 나누는 하드웨어 지원을 사용하여 나누기가 수행됩니다. 결과가 틀릴 수도 있지만, 그 결과가 기껏해야 1만큼 빗나갈 수 있다는 것을 증명할 수 있다고 생각합니다. 나눗셈 후에 제수(원래의 것)에 결과를 곱하고 그 곱을 배당금에서 뺍니다. 그런 다음 필요한 경우 제수를 더하거나 빼서 수정합니다(나눗셈의 결과가 1만큼 벗어난 경우).
64비트에서 32비트 나누기에 대한 하드웨어 지원을 활용하여 128비트 정수를 32비트 1로 나누기가 쉽습니다. 제수가 < 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보다 크면 좋은 소식이 없습니다. 나는 MSVC가 사용하고 있다고 생각하는 앞에서 설명한 절차에서 시프트 후 결과가 1 이하로 꺼져 있다는 완전한 증거를 가지고 있지 않습니다. 그러나 나는 버려지는 부분이 제수보다 적어도 2^31배 작고, 피제수가 2^64보다 작고 제수가 2^32-1보다 크다는 사실과 관련이 있다고 생각합니다. , 따라서 결과는 2^32보다 작습니다.
배당금에 128비트가 있으면 비트를 버리는 트릭이 작동하지 않습니다. 따라서 일반적으로 최상의 솔루션은 아마도 GJ 또는 caf에서 제안한 솔루션일 것입니다. (글쎄, 비트 폐기가 작동하더라도 아마도 가장 좋을 것입니다. 128 비트 정수에 대한 나눗셈, 곱셈 뺄셈 및 수정은 느릴 수 있습니다.)
또한 부동 소수점 하드웨어를 사용할 생각도 했습니다. x87 부동 소수점 단위는 길이가 분수 64비트인 80비트 정밀도 형식을 사용합니다. 64비트로 64비트 나누기의 정확한 결과를 얻을 수 있다고 생각합니다. (나머지는 직접적으로가 아니라 "MSVC 절차"에서와 같이 곱셈과 뺄셈을 사용하는 나머지). IF 피제수 >=2^64 및 < 2^128 부동 소수점 형식으로 저장하는 것은 "MSVC 절차"에서 최하위 비트를 버리는 것과 유사합니다. 누군가가 이 경우 오류를 증명하고 유용하다고 생각할 수 있습니다. GJ의 솔루션보다 빠를 가능성이 있는지는 모르겠지만 시도해 볼 가치가 있습니다.
솔루션은 정확히 해결하려는 대상에 따라 다릅니다.
예를 들어 링 모듈로 64비트 정수에서 산술을 수행하는 경우 Montgomerys 감소 를 사용하는 것이 매우 효율적입니다. 물론 이것은 동일한 계수를 여러 번 반복하고 링의 요소를 특수 표현으로 변환하는 것이 효과적이라고 가정합니다.
이 Montgomerys 감소 속도에 대한 대략적인 추정을 하자면: 2.4Ghz Core 2에서 64비트 모듈러스와 1600ns의 지수를 사용하여 모듈식 지수를 수행하는 오래된 벤치마크가 있습니다. 이 지수는 약 96개의 모듈식 곱셈을 수행합니다( 및 모듈식 감소) 따라서 모듈식 곱셈당 약 40 사이클이 필요합니다.
@caf가 수락한 답변은 정말 훌륭하고 높은 평가를 받았지만 몇 년 동안 볼 수 없었던 버그가 포함되어 있습니다.
그 솔루션과 다른 솔루션을 테스트하는 데 도움이 되도록 테스트 도구를 게시하고 커뮤니티 위키로 만들고 있습니다.
unsigned cafMod(unsigned A, unsigned B) {
assert(B);
unsigned X = B;
// while (X < A / 2) { Original code used <
while (X <= A / 2) {
X <<= 1;
}
while (A >= B) {
if (A >= X) A -= X;
X >>= 1;
}
return A;
}
void cafMod_test(unsigned num, unsigned den) {
if (den == 0) return;
unsigned y0 = num % den;
unsigned y1 = mod(num, den);
if (y0 != y1) {
printf("FAIL num:%x den:%x %x %x\n", num, den, y0, y1);
fflush(stdout);
exit(-1);
}
}
unsigned rand_unsigned() {
unsigned x = (unsigned) rand();
return x * 2 ^ (unsigned) rand();
}
void cafMod_tests(void) {
const unsigned i[] = { 0, 1, 2, 3, 0x7FFFFFFF, 0x80000000,
UINT_MAX - 3, UINT_MAX - 2, UINT_MAX - 1, UINT_MAX };
for (unsigned den = 0; den < sizeof i / sizeof i[0]; den++) {
if (i[den] == 0) continue;
for (unsigned num = 0; num < sizeof i / sizeof i[0]; num++) {
cafMod_test(i[num], i[den]);
}
}
cafMod_test(0x8711dd11, 0x4388ee88);
cafMod_test(0xf64835a1, 0xf64835a);
time_t t;
time(&t);
srand((unsigned) t);
printf("%u\n", (unsigned) t);fflush(stdout);
for (long long n = 10000LL * 1000LL * 1000LL; n > 0; n--) {
cafMod_test(rand_unsigned(), rand_unsigned());
}
puts("Done");
}
int main(void) {
cafMod_tests();
return 0;
}
질문에 지정된 32비트 코드를 알고 있지만 64비트에 대한 답변이 다른 사람들에게 유용하거나 흥미로울 수 있습니다.
그리고 예, 64b/32b => 32b 나눗셈은 128b % 64b => 64b에 대한 유용한 빌딩 블록을 만듭니다. libgcc __umoddi3(아래 링크된 소스)는 그런 종류의 작업을 수행하는 방법에 대한 아이디어를 제공하지만 4N % 2N => 2N이 아닌 2N / N => N 나눗셈 위에 2N % 2N => 2N만 구현합니다.
더 넓은 다중 정밀도 라이브러리를 사용할 수 있습니다(예: https://gmplib.org/manual/Integer-Division.html#Integer-Division ) .
64비트 시스템의 GNU C 는 대상 아키텍처에서 가능한 한 효율적으로 곱하고 나눌 수 있는 __int128type 및 libgcc 함수를 제공합니다 .
x86-64의 div r/m64명령어는 128b/64b => 64b 나눗셈을 수행하지만(두 번째 출력으로 나머지도 생성) 몫이 오버플로되면 오류가 발생합니다. 따라서 if 를 직접 사용할 수는 없지만 A/B > 2^64-1gcc가 이를 사용하도록 할 수 있습니다(또는 libgcc가 사용하는 동일한 코드를 인라인할 수도 있습니다).
이것은 하나 또는 두 개의 명령어( libgcc 함수 호출 내에서 발생 )로 컴파일( Godbolt 컴파일러 탐색기 )합니다 . 더 빠른 방법이 있다면 libgcc가 대신 사용할 것입니다.div
#include <stdint.h>
uint64_t AmodB(unsigned __int128 A, uint64_t B) {
return A % B;
}
이 __umodti3함수가 호출 하는 함수는 전체 128b/128b 모듈로를 계산하지만 해당 함수의 구현은 libgcc 소스에서 볼 수 있듯이 제수의 상위 절반이 0인 특수한 경우를 확인합니다 . (libgcc는 대상 아키텍처에 맞게 해당 코드에서 함수의 si/di/ti 버전을 빌드합니다. 대상 아키텍처에 대해 udiv_qrnnd부호 없는 2N/N => N 분할을 수행하는 인라인 asm 매크로입니다.
x86-64 (및 하드웨어 분할 명령이 있는 다른 아키텍처)의 경우 빠른 경로 (때 high_half(A) < B; div오류가 없음을 보장하는 경우) 는 두 개의 취하지 않은 분기 , 순서가 잘못된 CPU가 씹을 수 있는 약간의 보풀 및 div r64Agner Fog의 insn 테이블 에 따르면 최신 x86 CPU에서 약 50-100 사이클이 소요 되는 단일 명령 입니다. 다른 작업이 와 병렬로 발생할 수 div있지만 정수 나누기 장치는 파이프라인이 많지 않고 div많은 uop으로 디코딩됩니다(FP 나누기와 달리).
폴백 경로 는 64비트에 불과하지만 64비트에 맞지 않아 직접 오류가 발생 div하는 경우에 대해 두 개의 64비트 명령어 만 사용합니다 .BA/BA/B
libgcc 는 나머지만 반환하는 래퍼로 __umodti3인라인 __udivmoddi4됩니다.
동일한 모듈로 반복 B
그것은 계산 고려 가치가있을 수도 고정 소수점 역수를 위해 B존재하는 경우. 예를 들어, 컴파일 시간 상수를 사용하면 gcc는 128b보다 좁은 유형에 대해 최적화를 수행합니다.
uint64_t modulo_by_constant64(uint64_t A) { return A % 0x12345678ABULL; }
movabs rdx, -2233785418547900415
mov rax, rdi
mul rdx
mov rax, rdx # wasted instruction, could have kept using RDX.
movabs rdx, 78187493547
shr rax, 36 # division result
imul rax, rdx # multiply and subtract to get the modulo
sub rdi, rax
mov rax, rdi
ret
x86의 mul r64명령어는 64b*64b => 128b (rdx:rax) 곱셈을 수행하며 동일한 알고리즘을 구현하기 위해 128b * 128b => 256b 곱셈을 구성하는 빌딩 블록으로 사용할 수 있습니다. 전체 256b 결과의 상위 절반만 필요하므로 몇 배로 절약할 수 있습니다.
최신 Intel CPU는 매우 높은 성능을 제공합니다 mul. 3c 지연 시간, 클록 처리량당 하나씩. 그러나 필요한 시프트와 덧셈의 정확한 조합은 상수에 따라 다르므로 런타임에 곱셈 역을 계산하는 일반적인 경우는 JIT 컴파일 또는 정적으로 컴파일된 버전(심지어 사전 계산 오버헤드에 추가).
손익분기점이 될 IDK. JIT 컴파일의 경우 일반적으로 사용되는 B값에 대해 생성된 코드를 캐시하지 않는 한 ~200회 재사용보다 높습니다 . "정상적인" 방법의 경우 200회 재사용 범위에 있을 수 있지만 128비트/64비트 나눗셈에 대한 모듈식 곱셈 역함수를 찾는 것이 얼마나 비용이 많이 드는지 IDK입니다.
libdivide 는 이 작업을 수행할 수 있지만 32비트 및 64비트 유형에만 해당됩니다. 그래도 아마도 좋은 출발점이 될 것입니다.
일반적으로 나눗셈은 느리고 곱셈은 빠르며 비트 시프팅은 아직 빠릅니다. 지금까지 내가 본 답변에서 대부분의 답변은 비트 시프트를 사용하는 무차별 대입 접근 방식을 사용하고 있습니다. 다른 방법이 있습니다. 더 빠른지는 두고 봐야 합니다(일명 프로파일링).
나누는 대신에 역수를 곱하십시오. 따라서 A % B를 발견하려면 먼저 B ... 1/B의 역수를 계산하십시오. 이것은 Newton-Raphson 수렴 방법을 사용하여 몇 개의 루프로 수행할 수 있습니다. 이 작업을 잘 수행하려면 테이블의 초기 값 집합에 따라 달라집니다.
Newton-Raphson의 역수 수렴 방법에 대한 자세한 내용은 http://en.wikipedia.org/wiki/Division_(digital)을 참조하십시오 .
역수를 구하면 몫 Q = A * 1/B입니다.
나머지 R = A - Q*B.
이것이 무차별 대입보다 빠른지 확인하려면(64비트 및 128비트 숫자를 시뮬레이션하기 위해 32비트 레지스터를 사용할 것이기 때문에 더 많은 곱셈이 있으므로 프로파일링합니다.
B가 코드에서 상수이면 역수를 미리 계산하고 마지막 두 공식을 사용하여 간단히 계산할 수 있습니다. 이것은 비트 시프팅보다 빠를 것이라고 확신합니다.
도움이 되었기를 바랍니다.
최신 x86 시스템이 있는 경우 SSE2+용 128비트 레지스터가 있습니다. 나는 기본 x86 이외의 다른 어셈블리를 작성하려고 시도한 적이 없지만 거기에 몇 가지 가이드가 있다고 생각합니다.
부호 없는 128비트와 부호 없는 63비트로 충분하다면 최대 63 사이클을 수행하는 루프에서 수행할 수 있습니다.
이것을 1비트로 제한하여 제안된 솔루션 MSN의 오버플로 문제를 고려하십시오. 우리는 문제를 2, 모듈식 곱셈으로 나누고 끝에 결과를 더함으로써 그렇게 합니다.
다음 예에서 upper는 최상위 64비트에 해당하고 lower는 최하위 64비트에 해당하며 div는 제수입니다.
unsigned 128_mod(uint64_t upper, uint64_t lower, uint64_t div) {
uint64_t result = 0;
uint64_t a = (~0%div)+1;
upper %= div; // the resulting bit-length determines number of cycles required
// first we work out modular multiplication of (2^64*upper)%div
while (upper != 0){
if(upper&1 == 1){
result += a;
if(result >= div){result -= div;}
}
a <<= 1;
if(a >= div){a -= div;}
upper >>= 1;
}
// add up the 2 results and return the modulus
if(lower>div){lower -= div;}
return (lower+result)%div;
}
유일한 문제는 제수가 64비트인 경우 잘못된 결과를 제공하는 1비트의 오버플로(정보 손실)가 발생한다는 것입니다.
오버플로를 처리하는 깔끔한 방법을 찾지 못했다는 사실이 저를 괴롭힙니다.
C에는 미리 정의된 128비트 정수 유형이 없기 때문에 A의 비트는 배열로 표현되어야 합니다. B(64비트 정수)는 unsigned long long int 변수에 저장할 수 있지만 A와 B에서 효율적으로 작업하려면 B의 비트를 다른 배열에 넣어야 합니다.
그 후, B는 A보다 작은 가장 큰 B가 될 때까지 Bx2, Bx3, Bx4, ...로 증가합니다. 그런 다음 밑수 2에 대한 일부 빼기 지식을 사용하여 (AB)를 계산할 수 있습니다.
이것이 당신이 찾고 있는 솔루션의 종류입니까?
ReferenceURL : https://stackoverflow.com/questions/2566010/fastest-way-to-calculate-a-128-bit-integer-modulo-a-64-bit-integer
'IT이야기' 카테고리의 다른 글
| curl과 함께 PATCH 동사를 사용하는 방법 (0) | 2021.09.23 |
|---|---|
| Windows의 전역 npm 설치 위치 (0) | 2021.09.23 |
| 템플릿 구조/클래스를 친구로 선언하는 방법 (0) | 2021.09.23 |
| runST 및 함수 구성 (0) | 2021.09.23 |
| 바이트 대신 Enum Int32의 기본 유형을 만들어야 하는 이유 (0) | 2021.09.23 |