La manière la plus rapide de calculer un nombre entier de 128 bits modulo un nombre entier de 64 bits
j'ai un entier a non signé de 128 bits et un entier B non signé de 64 bits. Quelle est la façon la plus rapide de calculer A % B -c'est le reste (64 bits) de diviser A par B?
je cherche à le faire en langage C ou assembleur, 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 un seul instruction.
Edit:
je vous Remercie pour les réponses. Cependant, il me semble que les algorithmes suggérés seraient assez lents - la manière la plus rapide d'effectuer une division de 128 bits par 64 bits ne serait-elle pas d'utiliser le support natif du processeur pour la division de 64 bits par 32 bits? Est-ce que quelqu'un sait s'il y a une façon d'effectuer la division plus grande en termes de quelques divisions plus petites?
Re: Combien de fois B le changement?
principalement je suis 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 autant que 200 à diviser par chaque B. En quoi votre réponse différerait-elle dans ce cas?
13 réponses
vous pouvez utiliser la version de division de la Multiplication paysanne russe .
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 aurez besoin de mettre en œuvre les déplacements, les comparaisons et les soustractions pour opérer sur des valeurs composées d'une paire de nombres 64 bits, mais c'est assez trivial.
ce sera boucle au plus 255 fois (avec un 128 bits). Bien sûr, vous devez faire une pré-enregistrement pour un diviseur de zéro.
peut-être que vous êtes à la recherche d'un programme fini, mais les algorithmes de base pour l'arithmétique multi-précision peut être trouvé dans Art of Computer Programming , Volume 2. Vous pouvez trouver l'algorithme de division décrit en ligne ici . Les algorithmes traitent de multi-précision arithmétique arbitraire, et sont donc plus générale que vous avez besoin, mais vous devriez être en mesure de les simplifier pour 128 bit arithmétique fait sur 64 ou 32 bits chiffres. Être préparé pour une quantité raisonnable de travail (a) comprendre l'algorithme, et (b) Le convertir en C ou assembleur.
vous pourriez également vouloir vérifier Hacker's Delight , qui est plein de très intelligent assembleur et d'autres hackery de bas niveau, y compris certains arithmétique multi-précision.
Donné 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 supporte des entiers 64 bits, alors c'est probablement la façon la plus facile de faire.
La mise en œuvre par MSVC d'un modulo 64 bits sur 32 bits x86 est un assemblage en boucle ( VC\crt\src\intel\llrem.asm pour les braves), donc je dirais personnellement.
c'est presque non testé partiellement la vitesse modicated Mod128by64 'fonction d'algorithme paysan russe. Malheureusement je suis un utilisateur Delphi donc cette fonction fonctionne sous Delphi. :) Mais l'assembleur est presque la même...
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 > "151900920"FFFFFF:FFFFFFFF
mov ecx, 8 //Load byte shift counter
mov esi, [eax + 12] //Do fast 56 bit (7 bytes) shift...
shr esi, cl //esi = "151900920"XXXXXX
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 '"énorme diviseur nombre D'Optimisation de décalage", nous pouvons tester les diviseurs high bit, s'il est 0, nous n'avons pas besoin d'utiliser le Registre supplémentaire bh en tant que 65th bit pour le stocker. Donc déroulé le cadre de boucle peut ressembler à:
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:
j'aimerais partager quelques réflexions.
ce n'est pas aussi simple que MSN Le propose j'en ai peur.
dans l'expression:
(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
la multiplication et l'addition peuvent déborder. Je pense qu'on pourrait en tenir compte et utiliser le concept général avec quelques modifications, mais quelque chose me dit que ça va devenir vraiment effrayant.
j'étais curieux de savoir comment l'opération modulo 64 bits a été implémentée en MSVC et moi avons essayé de trouver quelque chose. Je ne sais pas vraiment assembler et tout ce que j'avais disponible était L'édition Express, sans la source de VC\crt\src\intel\llrem.asm, mais je pense que j'ai réussi à avoir une idée de ce qui se passe, après un peu de jouer avec le débogueur et la sortie de démontage. J'ai essayé de comprendre comment le reste est calculé dans le cas des nombres entiers positifs et le diviseur >=2^32. Il y a du code qui traite des nombres négatifs, bien sûr, mais je n'ai pas creuser dans.
Voici comment je le vois:
Si diviseur >= 2^32 fois le dividende et le diviseur sont décalés à droite autant que nécessaire pour s'adapter le diviseur en 32 bits. En d'autres termes: s'il faut n chiffres pour écrire le diviseur en binaire et n > 32, N-32 les chiffres les moins significatifs à la fois du diviseur et du dividende sont écartés. Après cela, la division est effectuée en utilisant le support matériel pour diviser des entiers 64 bits par des entiers 32 bits. Le résultat pourrait être incorrect, mais je pense qu'il peut être prouvé que le résultat est peut-être éteint au plus à 1. Après la division, le diviseur (original) est multiplié par le résultat et le produit soustrait du dividende. Puis c'est corrigé en ajoutant ou en soustrayant le diviseur si nécessaire (si le résultat de la division a été désactivé par un).
il est facile de diviser 128 bits entier par 32 bits un en tirant parti de la prise en charge du matériel 64 bits par division 32 bits. Dans le cas où le diviseur < 2^32, on peut calculer les autres n'effectuent que 4 divisions comme suit:
supposons que le dividende soit stocké dans:
DWORD dividend[4] = ...
le reste ira dans:
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.
après ces 4 étapes le reste de la variable tiendra ce que vous recherchez. (S'il vous plaît, ne me tuez pas si je me suis trompé sur l'Enness. Je ne suis même pas programmeur)
dans le cas où le diviseur est plus grand que 2^32-1 Je n'ai pas de bonnes nouvelles. Je n'ai pas de preuve complète que le résultat après le quart de travail n'est pas éteint de plus de 1, dans la procédure que j'ai décrite plus tôt, que je crois que MSVC utilise. Je pense cependant qu'il a quelque chose à voir avec le fait, que la partie qui est écartée est au moins 2^31 fois moins que le diviseur, le dividende est moins de 2^64 et le diviseur est plus grand que 2^32-1, de sorte que le résultat est moins de 2^32.
si le dividende a 128 bits l'astuce avec les bits de rejet ne fonctionnera pas. Donc, en général, cas la meilleure solution est probablement celle proposée par GJ ou caf. (Eh bien, ce serait probablement le meilleur même si le fait de jeter des bouts fonctionnait. Division, multiplication, soustraction et correction sur 128 bits entier peut être plus lent.)
je pensais aussi utiliser le matériel à virgule flottante. l'unité à virgule flottante x87 utilise un format de précision de 80 bits avec une fraction de 64 bits de long. Je pense que l'on peut obtenir le résultat exact de 64 bits de 64 bits division. (Pas le reste directement de, mais aussi le reste utilisant la multiplication et la soustraction comme dans la "procédure MSVC"). Si le dividende >=2^64 et < 2^128 le stocker dans le format point flottant semble similaire à jeter les bits les moins significatifs dans la "procédure MSVC". Peut-être quelqu'un peut prouver l'erreur dans ce cas est lié et le trouvez utile. Je ne sais pas s'il a une chance d'être plus rapide que la solution de GJ, mais peut-être que ça vaut le coup d'essayer.
La solution dépend exactement ce que vous essayez de résoudre.
E. G. si vous faites l'arithmétique dans un anneau modulo un entier 64-bit alors en utilisant Montgomerys réduction est très efficace. Bien sûr, cela suppose que vous le même module plusieurs fois et qu'il paie pour 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 Montgomerys: j'ai un vieux benchmark qui effectue une exponentiation modulaire avec un module 64 bits et un exposant en 1600 ns sur un noyau 2,4 Ghz 2. Cette exponentiation produit environ 96 multiplications modulaires (et réductions modulaires) et nécessite donc environ 40 cycles par multiplication modulaire.
j'ai fait les deux versions de Mod128by64 'fonction de division paysanne russe: classique et vitesse optimisée. La vitesse optimisée peut faire sur mon PC 3Ghz plus de 1000.000 calculs aléatoires par seconde et est plus de trois fois plus rapide que la fonction classique. Si nous comparons le temps d'exécution de calcul 128 par 64 et de calcul 64 Par 64 bits modulo que cette fonction est seulement environ 50% plus lent.
paysan russe classique:
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;
paysan russe optimisé en vitesse:
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 > "151910920"FFFFFF:FFFFFFFF
mov ecx, 8 //Load byte shift counter
mov esi, [eax + 12] //Do fast 56 bit (7 bytes) shift...
shr esi, cl //esi = "151910920"XXXXXX
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;
la réponse acceptée par @caf était vraiment sympa et hautement cotée, mais elle contient un bug qui n'a pas été vu depuis des années.
pour aider à tester cette solution et d'autres, je poste un harnais de test et je le fais wiki communautaire.
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;
}
je sais que la question a spécifié un code 32 bits, mais la réponse pour 64 bits peut être utile ou intéressante pour d'autres.
and yes, 64b/32b => 32b division does make a useful building-block for 128b % 64b => 64b. le __umoddi3 de libgcc (source liée ci-dessous) donne une idée de comment faire ce genre de chose, mais il n'implémente que 2N % 2n => 2N en plus d'une division 2N / N => N, pas 4N % 2n => 2n.
plus Large des bibliothèques multi-précision sont disponibles, par exemple https://gmplib.org/manual/Integer-Division.html#Integer-Division .
GNU C sur les machines 64 bits fournit un __int128 type , et les fonctions libgcc pour multiplier et diviser aussi efficacement que possible sur l'architecture cible.
x86-64 div r/m64 l'instruction fait 128b / 64b = > division 64b (produisant également le reste comme une seconde sortie), mais il se trompe si le quotient déborde. Donc vous ne pouvez pas l'utiliser directement si A/B > 2^64-1 , mais vous pouvez faire en sorte que gcc l'utilise pour vous (ou même en ligne le même code que libgcc utilise).
Cette compile ( Godbolt compilateur explorer ), à un ou deux div instructions (qui se produisent à l'intérieur d'un libgcc appel de fonction). S'il y avait un moyen plus rapide, libgcc l'utiliserait probablement à la place.
#include <stdint.h>
uint64_t AmodB(unsigned __int128 A, uint64_t B) {
return A % B;
}
la fonction __umodti3 qu'elle appelle calcule un modulo 128b/128b complet, mais la mise en œuvre de cette fonction vérifie le cas spécial où la moitié haute du diviseur est 0, comme vous pouvez voir dans le Source libgcc . (libgcc construit la version si/di/ti de la fonction à partir de ce code, comme approprié pour l'architecture cible. udiv_qrnnd est une macro asm en ligne qui n'est pas signée 2N/N => n division pour l'architecture cible.
pour x86-64 (et d'autres architectures avec une instruction de division du matériel), le fast-path (quand high_half(A) < B ; en garantissant div ne sera pas fautif) est juste deux branches non-prises , un certain fluff pour les CPU hors-ordre à mâcher à travers, et un simple div r64 instruction, qui prend environ 50-100 cycles sur les CPU x86 modernes, selon Agner Fog's insn tables . Certains autres travaux peuvent se produire en parallèle avec div , mais l'Unité de division entière n'est pas très pipeliné et div décode à beaucoup de uops (contrairement à la division FP).
le chemin de repli utilise encore seulement deux 64 bits div instructions pour le cas où B est seulement 64 bits, mais A/B ne rentre pas en 64 bits donc A/B serait directement fautif.
notez que le __umodti3 de libgcc" just inlines __udivmoddi4 dans un wrapper qui ne renvoie que le reste.
répétés modulo par le même B
il pourrait être intéressant d'envisager le calcul d'un point fixe multiplicatif inverse pour B , si un existe. Par exemple, avec les constantes de temps de compilation, gcc fait l'optimisation pour les types plus étroits que 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 instruction n'64*64 => 128b (rdx:rax) la multiplication, et peut être utilisé comme un bloc de construction pour construire un 128b * 128b => 256b se multiplient pour mettre en œuvre le même algorithme. Puisque nous n'avons besoin que de la moitié supérieure du résultat 256b complet, cela permet d'économiser quelques multiples.
les processeurs Intel modernes sont très performants mul : latence 3c, un par cadran. Cependant, la combinaison exacte des déplacements et des additions varie avec la constante, de sorte que le cas général de calcul d'une inversion multiplicative à l'exécution du temps n'est pas tout à fait aussi efficace chaque fois qu'il est utilisé comme une JIT-compiled ou statically-compiled version (même en plus du pré-calcul overhead).
IDK où le point d'équilibre serait. Pour la compilation JIT, il sera supérieur à ~200 réutilisations, à moins que vous ne mettiez en cache code généré pour les valeurs B couramment utilisées. Pour la voie "normale", il pourrait être dans la gamme de 200 réutilisations, mais IDK combien il serait coûteux de trouver un inverse multiplicatif modulaire pour 128-bit / 64-bit de division.
libdivide peut le faire pour vous, mais seulement pour les types 32 et 64 bits. Pourtant, c'est probablement un bon point de départ.
en règle générale, la division est lente et la multiplication est plus rapide, et le transfert de bits est encore plus rapide. D'après ce que j'ai vu des réponses jusqu'à présent, la plupart des réponses ont utilisé une approche de la force brute en utilisant des changements de bits. Il existe une autre voie. Si il est plus rapide reste à voir (AKA profil).
au lieu de diviser, multipliez par le réciproque. Ainsi, pour découvrir un % B, calculez d'abord la réciproque de B... 1 / B. Ceci peut être fait avec quelques boucles utilisant la méthode de convergence Newton-Raphson. Pour ce faire il dépendra d'un bon jeu de valeurs initiales dans un tableau.
pour plus de détails sur la méthode Newton-Raphson de convergence sur la réciproque, veuillez vous référer à http://en.wikipedia.org/wiki/Division_ (numérique)
une fois que vous avez la réciproque, le quotient Q = A * 1 / B.
le reste R = a-Q*B.
à déterminez si cela serait plus rapide que la force brute (car il y aura beaucoup plus de multiplications puisque nous utiliserons des registres 32 bits pour simuler des nombres 64 bits et 128 bits, profilez-le.
Si B est constant dans votre code, vous pouvez pré-calculer le réciproque et simplement calculer en utilisant les deux dernières formules. Cela, je suis sûr sera plus rapide que bit-shifting.
Espérons que cette aide.
si vous avez une machine x86 récente, il y a des registres de 128 bits pour SSE2+. Je n'ai jamais essayé d'écrire assembly pour autre chose que le x86 de base, mais je soupçonne qu'il y a des guides là-bas.
si 128-bit non signé par 63-bit non signé est assez bon, alors il peut être fait dans une boucle faisant au plus 63 cycles.
une solution proposée Msn' problème de dépassement en la limitant à 1 bit. Nous n'avons donc diviser le problème en 2, multiplication modulaire et en additionnant les résultats à la fin.
dans l'exemple suivant, upper correspond au 64 bits le plus significatif, lower au 64 bits le moins significatif et div est le diviseur.
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;
}
le seul problème est que, si le diviseur est 64 bits alors nous obtenons des débordements de 1 bit (perte d'information) donnant un résultat fautif.
ça me dérange que je n'ai pas trouvé une façon soignée de gérer les débordements.
Puisqu'il n'y a pas de type entier 128 bits prédéfini dans C, les bits de A doivent être représentés dans un tableau. Bien que B (entier 64 bits) puisse être stocké dans une variable non signé long int , il est nécessaire de mettre des bits de B dans un autre tableau afin de travailler sur A et B efficacement.
après cela, B est incrémenté comme Bx2, Bx3, Bx4,... jusqu'à ce Qu'il soit le plus grand B de moins que A. et puis (A-B) peut être calculé, en utilisant une certaine connaissance de soustraction pour la base 2.
Est-ce le genre de solution que vous recherchez?