Le moyen le plus rapide de calculer un nombre entier de 128 bits modulo un nombre entier de 64 bits

Question

J'ai un entier non signé de 128 bits A et un entier non signé de 64 bits B. Quel est le moyen le plus rapide pour calculer A % B - qui est le reste (64 bits) de la division de A par B ?

Je cherche à le faire en C ou en langage d'assemblage, mais je dois cibler la plate-forme x86 32 bits. Cela signifie malheureusement que je ne peux pas profiter du support du compilateur pour les entiers de 128 bits, ni de la capacité de l'architecture x64 à effectuer l'opération requise en une seule instruction.

Editar:

Merci pour les réponses apportées jusqu'à présent. Cependant, il me semble que les algorithmes suggérés seraient assez lents - le moyen le plus rapide d'effectuer une division de 128 bits par 64 bits ne serait-il pas d'exploiter le support natif du processeur pour la division de 64 bits par 32 bits ? Quelqu'un sait-il s'il existe un moyen d'effectuer la plus grande division en termes de plusieurs petites divisions ?

Re : A quelle fréquence le B change-t-il ?

Je suis avant tout intéressé par une solution générale - quel calcul effectueriez-vous si A et B sont susceptibles d'être différents à chaque fois ?

Cependant, une deuxième situation possible est que B ne varie pas aussi souvent que A - il peut y avoir jusqu'à 200 As à diviser par chaque B. En quoi votre réponse serait-elle différente dans ce cas ?

Answer 1

Vous pouvez utiliser la version divisionnaire de Multiplication des paysans russes .

Pour trouver le reste, exécutez (en pseudo-code) :

X = B;

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

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

Le module est laissé en A.

Vous devrez implémenter les décalages, les comparaisons et les soustractions pour opérer sur des valeurs composées d'une paire de nombres de 64 bits, mais c'est assez trivial (il est probable que vous deviez implémenter le décalage à gauche par 1 en tant que X + X ).

Cette boucle sera effectuée au maximum 255 fois (avec un A de 128 bits). Bien sûr, vous devez faire une vérification préalable pour un diviseur nul.

Answer 2

Vous cherchez peut-être un programme fini, mais les algorithmes de base de l'arithmétique multiprécision se trouvent dans l'ouvrage de Knuth intitulé Art de la programmation informatique Volume 2. Vous pouvez trouver l'algorithme de division décrit en ligne aquí . Les algorithmes traitent de l'arithmétique multiprécision arbitraire, et sont donc plus généraux que ce dont vous avez besoin, mais vous devriez être capable de les simplifier pour l'arithmétique 128 bits effectuée sur des chiffres de 64 ou 32 bits. Préparez-vous à une quantité raisonnable de travail (a) pour comprendre l'algorithme, et (b) pour le convertir en C ou en assembleur.

Vous pouvez également consulter Le plaisir du hacker qui est plein d'assembleur très intelligent et d'autres piratages de bas niveau, y compris de l'arithmétique multiprécision.

Answer 3

Si votre B est assez petit pour que le uint64_t + l'opération de ne pas emballer :

Dado A = AH*2^64 + AL :

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

Si votre compilateur prend en charge les entiers 64 bits, c'est probablement le moyen le plus simple de procéder. L'implémentation par MSVC d'un modulo 64 bits sur x86 32 bits est un assemblage de boucles fastidieuses ( VC\crt\src\intel\llrem.asm pour les courageux), donc je choisirais personnellement cette solution.

Answer 4

Il s'agit d'une fonction de l'algorithme du "paysan russe" Mod128by64 partiellement modifiée par la vitesse et pratiquement non testée. Malheureusement je suis un utilisateur de Delphi donc cette fonction fonctionne sous Delphi :) Mais l'assembleur est presque le même donc...

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;

Au moins une autre optimisation de la vitesse est possible ! Après l'optimisation du décalage des nombres à grand diviseur, nous pouvons tester le bit haut des diviseurs, s'il est égal à 0, nous n'avons pas besoin d'utiliser le registre supplémentaire bh comme 65ème bit pour le stocker. Ainsi, la partie déroulée de la boucle peut ressembler à ceci :

  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:

Answer 5

La solution dépend de ce que vous essayez de résoudre exactement.

Par exemple, si vous faites de l'arithmétique dans un anneau modulo un entier de 64 bits, alors l'utilisation de Réduction de Montgomery est très efficace. Bien entendu, cela suppose que vous utilisiez le même module plusieurs fois et qu'il est intéressant de convertir les éléments de l'anneau en une représentation spéciale.

---

Pour donner une estimation très approximative de la vitesse de cette réduction de Montgomery : J'ai un vieux benchmark qui effectue une exponentiation modulaire avec un module et un exposant de 64 bits en 1600 ns sur un Core 2 de 2.4Ghz. Cette exponentiation fait environ 96 multiplications modulaires (et réductions modulaires) et nécessite donc environ 40 cycles par multiplication modulaire.