programing

128비트 정수 모듈로 64비트 정수를 계산하는 가장 빠른 방법

javaba 2022. 8. 3. 22:39
반응형

128비트 정수 모듈로 64비트 정수를 계산하는 가장 빠른 방법

B.128이 있습니다. 계산할 수 있는 A % BA'B'는 (64)는요?

C 또는 어셈블리 언어 중 하나를 사용하고 싶은데 32비트 x86 플랫폼을 대상으로 해야 합니다.이는 128비트 정수에 대한 컴파일러 지원이나 단일 명령으로 필요한 작업을 수행하는 x64 아키텍처의 기능을 이용할 수 없음을 의미합니다.

편집:

지금까지 답변 감사합니다.단, 128비트 x 64비트 분할을 실행하는 가장 빠른 방법은 64비트 x 32비트 분할에 대한 프로세서의 네이티브 지원을 활용하는 것이 아닐까요?작은 사업부 몇 개에 대해 큰 사업부를 수행할 수 있는 방법이 있는지 아는 사람 있나요?

B는 얼마나 자주 변합니까?

저는 주로 일반적인 솔루션에 관심이 있습니다.A와 B가 매번 다를 가능성이 높다면 어떻게 계산하시겠습니까?

그러나 두 번째 가능한 상황은 B가 A만큼 자주 변화하지 않는다는 것입니다. 각 B로 나누는 A는 200개까지 있을 수 있습니다.이 경우 당신의 답변은 어떻게 다를까요?

러시아 농민 곱셈의 분할 버전을 사용할 수 있습니다.

나머지를 찾으려면 (의사 코드로) 다음을 수행합니다.

X = B;

while (X <= A/2)
{
    X <<= 1;
}

while (A >= B)
{
    if (A >= X)
        A -= X;
    X >>= 1;
}

계수는 A에 남습니다.

숫자의 , 은 매우을 64비트로서 가 있습니다).X + X를 참조해 주세요.

최대 255회(128비트 A) 루프합니다.물론 0의 제수를 미리 확인해야 합니다.

32비트 코드가 지정되어 있는 질문은 알고 있습니다만, 64비트에 대한 답변은 다른 사람에게 유용하거나 흥미로울 수 있습니다.

=> % => 의 "64b/32b => 32b division" 128b % 64b => 64b.libgccs에 도움이 .__umoddi3(는 이러한 종류의에 대한 % > 디비전 2N / N > N % = > 이 아닙니다이 경우

https://gmplib.org/manual/Integer-Division.html#Integer-Division 등 폭넓은 멀티프로세서 라이브러리를 이용할 수 있습니다.


64비트 머신의 GNU C는 타겟아키텍처 상에서 가능한 한 효율적으로 증배 및 분할하기 위한 타입과 libgcc 함수를 제공합니다.

x86-64의 명령어는 128b/64b => 64b 나눗셈(나머지도 두 번째 출력으로 생성)을 수행하지만, 지수가 오버플로우하면 실패합니다.따라서 직접 사용할 수 없습니다.A/B > 2^64-1단, gcc를 사용하여 gcc를 사용할 수 있습니다(또는 libgcc가 사용하는 코드와 같은 인라인 코드도 가능합니다.


컴파일(Godbolt 컴파일러 탐색기)은 1개 또는 2개로 컴파일됩니다.div명령(libgcc 함수 호출 내에서 발생)더 빠른 방법이 있다면 libgcc가 대신 그것을 사용할 것입니다.

#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 는 인라인 asm 매크로로 타깃아키텍처에 대해 부호 없는 2N/N => N 분할을 수행합니다.

x86-64(및 하드웨어 분할 명령이 있는 기타 아키텍처)의 경우 패스트 패스는high_half(A) < B; 증증divAgner Fog의 insn 표에 따르면, 고장난 CPU를 씹을 수 있는 두 개의 가지 가지, 약간의 솜털, 그리고번의 명령으로, 최신 x86 CPU에서는 약 50-1001 사이클이 소요됩니다.다른 작업도 병행할 수 있습니다.div 있지 않고, ,, 지, 지, 지, 지, 지, 지, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,, ,,divFP분할하다

" " " " " " " " " " " " " " " " " " 64 " " " " " " " " " " " " " " " " " 을 합니다.div「 」의 B 64비트입니다A/B 않기 때문에 64비트를 사용할 수 없습니다.A/B직접적으로 잘못될 수 있습니다.

libgcc에 해 주세요.__umodti3__udivmoddi4나머지 부분만 반환하는 포장지로 만듭니다.

1: '1: 32'divCPU의 2개입니다.AMD CPU에서는 64비트 레지스터의 작은 값이라도 실제 입력 값의 크기에 따라 성능이 달라집니다.작은 값이 일반적인 경우 64비트 또는 128비트 분할을 수행하기 전에 지점을 단순한 32비트 분할 버전으로 벤치마킹하는 것이 좋습니다.


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

x86gsmul r64명령어는 64b*64b => 128b(rdx:rax) 곱셈을 수행하며, 동일한 알고리즘을 구현하기 위해 128b * 128b => 256b 곱셈을 구성 요소로 사용할 수 있습니다.전체 256b 결과의 상위 절반만 필요하므로 몇 배수를 절약할 수 있습니다.

인텔 는 매우 을 발휘합니다.mul: 3c 지연, 1클럭당 1회 스루풋그러나 필요한 이동과 추가의 정확한 조합은 상수에 따라 다르므로 런타임에 곱셈 역수를 계산하는 일반적인 경우는 JIT 컴파일 버전 또는 정적 컴파일 버전(컴퓨팅 전 오버헤드에도 해당)만큼 효율적이지 않습니다.

IDK 。JIT의 경우 되는 J에 대해 한 횟수가 입니다.IT 컴파일의 경우 일반적으로 사용되는 생성 코드를 캐시하지 않는 한 재사용 횟수가 최대 200회 이상입니다.B만, 에서는 128비트나눗셈의 알 수 있습니다.통상적인 방법으로는 200회 재사용 범위 내일 수 있지만 IDK에서는 128비트/64비트 분할의 모듈러 곱셈 역수를 찾는 것이 얼마나 비쌀지 알 수 있습니다.

libdivide는 32비트 및 64비트 유형에서만 이 기능을 수행할 수 있습니다.그래도 좋은 출발점이 될 수 있을 것 같아요.

서명되지 않은 128비트가 서명되지 않은 63비트로 충분할 경우 최대 63사이클을 실행하는 루프에서 실행할 수 있습니다.

이것은 MSN의 오버플로 문제를 1비트로 제한하여 제안하는 해결책이라고 생각하시면 됩니다.문제를 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비트의 오버플로우(정보 손실)가 발생하여 불량 결과가 발생한다는 것입니다.

넘쳐나는 것을 깔끔하게 처리하지 못한 것이 나를 짜증나게 한다.

어셈블러 코드를 컴파일하는 방법을 모르기 때문에 컴파일 및 테스트에 도움을 주시면 감사하겠습니다.

gmplib "mpz_mod()"와 비교하여 100만 개의 루프 결과를 합계함으로써 이 문제를 해결했습니다.속도 저하(시드업 0.12)에서 1.54까지 긴 여정이었습니다.그 때문에, 이 스레드의 C코드가 늦어질 것이라고 생각합니다.

스레드의 에 대한 내용은 다음과 같습니다.
https://www.raspberrypi.org/forums/viewtopic.php?f=33&t=311893&p=1873122#p1873122httpswww..org/forums/viewtopic.php?f=33&t=311893&p=1873122#p1873122

이것은 gmplib "mpz_mod()"를 사용하여 속도를 높인 "mod_256()"이며, 더 긴 시프트에 __builtin_clzll()을 사용해야 했습니다.

typedef __uint128_t uint256_t[2];

#define min(x, y) ((x<y) ? (x) : (y))

int clz(__uint128_t u)
{
//  unsigned long long h = ((unsigned long long *)&u)[1];
  unsigned long long h = u >> 64;
  return (h!=0) ? __builtin_clzll(h) : 64 + __builtin_clzll(u);
}

__uint128_t mod_256(uint256_t x, __uint128_t n)
{
  if (x[1] == 0)  return x[0] % n;
  else
  {
    __uint128_t r = x[1] % n;
    int F = clz(n);
    int R = clz(r);
    for(int i=0; i<128; ++i)
    {
      if (R>F+1)
      {
        int h = min(R-(F+1), 128-i);
        r <<= h; R-=h; i+=(h-1); continue;
      }
      r <<= 1; if (r >= n)  { r -= n; R=clz(r); }
    }
    r += (x[0] % n); if (r >= n)  r -= n;

    return r;
  }
}

가 B보다 uint64_t +다음 중 하나:

의 「」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 「 32」 「x86」 「64」)입니다.VC\crt\src\intel\llrem.asm용맹한 자를 위해) 그래서 저는 개인적으로 그렇게 하겠습니다.

@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;
}

Mod128by64 '러시아 농민' 분할 기능: 클래식과 속도를 최적화했습니다.최적화된 속도는 3Ghz PC에서 초당 1000.000회 이상의 랜덤 계산을 수행할 수 있으며 기존 기능보다 3배 이상 빠릅니다.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;

이것은 부분적인 속도 변조 Mod128by64 '러시아 농민' 알고리즘 함수에서 거의 테스트되지 않았습니다.유감스럽게도 저는 델파이 사용자이기 때문에 이 기능은 델파이에서도 동작합니다. :) 하지만 어셈블러는 거의 비슷하기 때문에...

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이면 65번째 비트로 추가 bh 레지스터를 사용할 필요가 없습니다.따라서 루프의 언롤 부분은 다음과 같습니다.

  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:

완성된 프로그램을 찾고 있을지도 모르지만, 다정밀 산술의 기본 알고리즘은 Knuth의 Art of Computer Programming, Volume 2에서 찾을 수 있습니다.여기서 온라인으로 설명되어 있는 나눗셈 알고리즘을 찾을 수 있습니다.알고리즘은 임의의 멀티 정밀도 산술에 대응하고 있기 때문에 필요 이상으로 일반적이지만, 64비트 또는 32비트 자릿수로 실행되는 128비트 산술에서는 단순화할 수 있습니다.(a)알고리즘을 이해하고 (b)C 또는 어셈블러로 변환할 수 있도록 준비한다.

당신은 또한 해커의 딜라이트를 체크하고 싶을지도 모른다.해커 딜라이트는 매우 영리한 어셈블러와 여러 가지 정밀도 계산을 포함한 다른 하위 수준의 해킹으로 가득 차 있다.

최신 x86 머신을 사용하는 경우 SSE2+용 128비트 레지스터가 있습니다.기본 x86 이외에는 조립품을 써 본 적이 없지만, 가이드도 있을 것 같습니다.

몇 가지 생각을 나누고 싶습니다.

유감스럽지만 그것은 MSN이 제안하는 것처럼 간단하지 않다.

식 중:

(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B

곱셈과 덧셈 모두 오버플로가 발생할 수 있습니다.그 점을 감안해서 일반적인 개념을 조금 수정해서 사용할 수 있을 것 같은데, 뭔가 굉장히 무서울 것 같아요.

MSVC에서 64비트 모듈로 연산이 어떻게 구현되는지 궁금해서 알아보려고 했습니다.어셈블리에 대해서는 잘 모르기 때문에 VC\crt\src\intel\lrem.asm 소스 없이 Express 에디션밖에 사용할 수 없었습니다만, 디버거와 디스어셈블리의 출력을 조금 조작해 본 결과, 무슨 일이 일어나고 있는지 알 수 있었다고 생각합니다.나는 정수와 제수가 >=2^32인 경우에 나머지를 어떻게 계산하는지 알아보려고 했다.물론 음수를 다루는 코드도 있지만, 저는 그것을 파고들지 않았습니다.

제 견해는 다음과 같습니다.

제수가 > = 2^32이면 배당수와 제수가 모두 32비트에 맞도록 필요한 만큼 오른쪽으로 이동한다.즉, 2진법으로 제수를 적는데 n자리 숫자가 필요하고 n > 32일 경우 제수와 배당 양쪽의 n-32자리 최소 유효 자릿수는 폐기됩니다.그 후 64비트 정수를 32비트 1로 분할하는 하드웨어 지원을 이용해 분할한다.결과는 틀릴 수도 있지만, 기껏해야 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-1보다 작은 것과 관계가 있다고 생각합니다.

배당금이 128비트일 경우 폐기 비트를 사용하는 트릭은 작동하지 않습니다.따라서 일반적으로 가장 좋은 해결책은 GJ나 카페가 제안하는 해결책일 것이다.(글쎄, 비트는 폐기해도 아마 최고일 거야.128비트 정수의 나눗셈, 곱셈 감산 및 보정은 더 느릴 수 있습니다.)

부동소수점 하드웨어를 사용할 생각도 하고 있습니다.x87 부동소수점 유닛은 길이 64비트의 80비트 정밀 포맷을 사용합니다.64비트 x 64비트 나눗셈의 정확한 결과를 얻을 수 있다고 생각합니다.(나머지는 직접적이지 않고 "MSVC 절차"와 같이 곱셈과 뺄셈을 사용하는 나머지도 마찬가지입니다).이를 부동소수점 형식으로 저장하는 배당 >=2^64 및 <2^128이 "MSVC 프로시저"에서 최하위 비트를 폐기하는 것과 유사한 것으로 보이는 경우.아마 누군가는 그 경우에 오류가 있다는 것을 증명하고 유용하다고 생각할 수 있을 것이다.GJ의 솔루션보다 더 빠를 가능성이 있는지는 모르겠지만, 시도해 볼 만한 가치가 있을지도 모릅니다.

해결 방법은 해결하려는 것이 정확히 무엇인지에 따라 달라집니다.

예를 들어 링 모듈로 64비트 정수를 연산하는 경우 Montgomerys 축소를 사용하는 것이 매우 효율적입니다.물론 이것은 같은 계수를 여러 번 반복하여 링의 요소를 특별한 표현으로 변환하는 것이 효과적이라고 가정합니다.


이 Montgomerys의 감소 속도에 대한 대략적인 견적을 제시하면 다음과 같습니다.2.4Ghz Core 2에서 64비트 계수 및 1600ns 단위의 지수를 사용하여 모듈식 지수화를 수행하는 오래된 벤치마크가 있습니다.이 지수화는 약 96개의 모듈식 곱셈(및 모듈식 감소)을 수행하므로 모듈식 곱셈당 약 40개의 사이클이 필요합니다.

일반적으로 나눗셈과 곱셈은 느리고 비트 이동은 더 빠릅니다.지금까지의 답변에서 본 바로는, 대부분의 답변은 비트 시프트를 사용한 무차별적인 어프로치를 사용하고 있습니다.다른 방법이 있어요더 빠른지 아닌지는 아직 알 수 없습니다(AKA 프로파일).

나눗셈 대신 역수를 곱하세요.따라서 A % B를 발견하려면 먼저 B ... 1/B의 역수를 계산합니다.이것은 뉴턴-라프슨 수렴법을 사용하여 몇 가지 루프로 이루어질 수 있습니다.이 작업을 제대로 수행하려면 테이블 내의 적절한 초기 값 집합에 의존합니다.

역수에서 수렴하는 뉴턴-라프슨 방법에 대한 자세한 내용은 http://en.wikipedia.org/wiki/Division_(digital)를 참조하십시오.

역수를 구하면 Q = A * 1/B입니다.

나머지 R = A - Q*B.

64비트 및 128비트 숫자를 시뮬레이트하기 위해 32비트 레지스터를 사용하기 때문에 이 속도가 브루트포스보다 빠를지 여부를 판단하려면 프로파일을 작성하십시오.

코드 B가 일정할 경우 역수를 미리 계산하고 마지막 두 공식을 사용하여 간단히 계산할 수 있습니다.이게 비트 시프트보다 빠를 거예요.

이게 도움이 됐으면 좋겠다.

나는 전투 후 9년이 지났지만 여기에 언급할 가치가 있는 2의 거듭제곱에 대한 흥미로운 O(1) 엣지 케이스가 있다.

#include <stdio.h>
// example with 32 bits and 8 bits.
int main() {
    int i = 930;
    unsigned char b = (unsigned char) i;
    printf("%d", (int) b); // 162, same as 930 % 256
}
  

C에는 미리 정의된 128비트 정수형이 없기 때문에 A의 비트는 배열로 표현해야 합니다.B(64비트 정수)는 부호 없는 int 변수에 저장할 수 있지만 A 및 B에서 효율적으로 동작하려면 B의 비트를 다른 배열에 배치해야 합니다.

그 후, B는 Bx2, Bx3, Bx4, ...로서 증가해, A보다 작은 최대 B가 된다.그런 다음 베이스 2에 대한 감산 지식을 사용하여 (A-B)를 계산할 수 있습니다.

이것이 당신이 찾고 있는 해결책입니까?

언급URL : https://stackoverflow.com/questions/2566010/fastest-way-to-calculate-a-128-bit-integer-modulo-a-64-bit-integer

반응형