c - tipos - unsigned int para que sirve
La forma más rápida de calcular un módulo entero de 128 bits un entero de 64 bits (13)
Como regla general, la división es lenta y la multiplicación es más rápida, y el cambio de bit es más rápido todavía. Por lo que he visto de las respuestas hasta ahora, la mayoría de las respuestas han usado un enfoque de fuerza bruta usando cambios de bits. Existe otra forma. Queda por ver si es más rápido (también conocido como perfil).
En lugar de dividir, multiplica por lo recíproco. Por lo tanto, para descubrir A% B, primero calcule el recíproco de B ... 1 / B. Esto se puede hacer con algunos bucles utilizando el método de convergencia de Newton-Raphson. Hacerlo bien dependerá de un buen conjunto de valores iniciales en una tabla.
Para obtener más detalles sobre el método de convergencia de Newton-Raphson en el recíproco, consulte http://en.wikipedia.org/wiki/Division_(digital)
Once you have the reciprocal, the quotient Q = A * 1/B.
The remainder R = A - Q*B.
To determine if this would be faster than the brute force (as there will be many more multiplies since we will be using 32-bit registers to simulate 64-bit and 128-bit numbers, profile it.
If B is constant in your code, you can pre-calculate the reciprocal and simply calculate using the last two formulae. This, I am sure will be faster than bit-shifting.
Espero que esto ayude.
Tengo un entero A de 128 bits sin signo y un entero sin signo de 64 bits B. ¿Cuál es la forma más rápida de calcular A % B , es decir, el resto (de 64 bits) de dividir A entre B?
Estoy buscando hacer esto en C o en lenguaje ensamblador, pero tengo que apuntarme a la plataforma x86 de 32 bits. Desafortunadamente, esto significa que no puedo aprovechar la compatibilidad del compilador para enteros de 128 bits, ni la capacidad de la arquitectura x64 para realizar la operación requerida en una sola instrucción.
Editar:
Gracias por las respuestas hasta ahora. Sin embargo, me parece que los algoritmos sugeridos serían bastante lentos: ¿la forma más rápida de realizar una división de 128 bits por 64 bits no sería aprovechar el soporte nativo del procesador para la división de 64 bits por 32 bits? ¿Alguien sabe si hay una manera de realizar la división más grande en términos de unas pocas divisiones más pequeñas?
Re: ¿Con qué frecuencia cambia B?
Principalmente estoy interesado en una solución general: ¿qué cálculo realizarías si A y B fueran diferentes cada vez?
Sin embargo, una segunda situación posible es que B no varía con tanta frecuencia como A - puede haber hasta 200 A para dividir por cada B. ¿En qué se diferencia su respuesta en este caso?
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 su compilador admite enteros de 64 bits, esta es probablemente la forma más fácil de hacerlo. La implementación de MSVC de un módulo de 64 bits en 32-bit x86 es un conjunto relleno de bucle peludo ( VC/crt/src/intel/llrem.asm para los valientes), así que yo personalmente iría con eso.
Esto es casi no comprobado en parte modificada por la velocidad Mod128by64 ''campesino ruso'' función del algoritmo. Lamentablemente, soy un usuario de Delphi por lo que esta función funciona en Delphi. :) Pero el ensamblador es casi lo mismo, así que ...
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;
¡Al menos una optimización de velocidad más es posible! Después de ''Enorme Divisor Numbers Shift Optimization'' podemos probar los divisores de alto bit, si es 0 no necesitamos usar el registro extra bh como el 65º bit para almacenarlo. Así que la parte desenrollada del ciclo puede verse así:
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:
He hecho ambas versiones de la función de división Mod128by64 ''campesino ruso'': clásico y velocidad optimizada. La velocidad optimizada puede hacer en mi PC 3Ghz más de 1000,000 cálculos aleatorios por segundo y es más de tres veces más rápido que la función clásica. Si comparamos el tiempo de ejecución de calcular 128 por 64 y calcular el módulo de 64 por 64, esta función es solo un 50% más lenta.
Campesino ruso clásico:
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;
Velocidad de campesino ruso optimizado:
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;
La respuesta aceptada por @caf fue muy buena y altamente calificada, sin embargo, contiene un error no visto durante años.
Para ayudar a probar esa y otras soluciones, estoy publicando un arnés de prueba y convirtiéndolo en wiki de la comunidad.
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;
}
La solución depende de qué es exactamente lo que estás tratando de resolver.
Por ejemplo, si está haciendo aritmética en un módulo de anillo, un entero de 64 bits, entonces usar la reducción de Montgomerys es muy eficiente. Por supuesto, esto supone que tiene el mismo módulo muchas veces y que vale la pena convertir los elementos del anillo en una representación especial.
Para dar solo una estimación aproximada de la velocidad de esta reducción de Montgomerys: Tengo un viejo punto de referencia que realiza una exponenciación modular con módulo de 64 bits y exponente en 1600 ns en un Core 2 de 2.4Ghz. Esta exponenciación hace aproximadamente 96 multiplicaciones modulares ( y reducciones modulares) y, por lo tanto, necesita alrededor de 40 ciclos por multiplicación modular.
Me gustaría compartir algunas ideas.
No es tan simple como me propone MSN. Tengo miedo.
En la expresión:
(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
tanto la multiplicación como la suma pueden desbordarse. Creo que uno podría tomarlo en cuenta y seguir usando el concepto general con algunas modificaciones, pero algo me dice que va a ser realmente aterrador.
Tenía curiosidad de cómo se implementó la operación de módulo de 64 bits en MSVC y traté de encontrar algo. Realmente no sé ensamblar y todo lo que tenía disponible era Express Edition, sin la fuente de VC / crt / src / intel / llrem.asm, pero creo que logré tener una idea de lo que está pasando, después de jugar un poco. con el resultado de depuración y desensamblaje. Traté de averiguar cómo se calcula el resto en el caso de los enteros positivos y el divisor> = 2 ^ 32. Hay un código que trata con números negativos, por supuesto, pero no profundicé en eso.
Así es como lo veo:
Si divisor> = 2 ^ 32 tanto el dividendo como el divisor se desplazan hacia la derecha tanto como sea necesario para ajustar el divisor en 32 bits. En otras palabras: si se necesitan n dígitos para escribir el divisor en binario yn> 32, se descartan n-32 dígitos menos significativos tanto del divisor como del dividendo. Después de eso, la división se realiza utilizando soporte de hardware para dividir enteros de 64 bits por otros de 32 bits. El resultado puede ser incorrecto, pero creo que puede probarse, que el resultado puede estar apagado a lo sumo 1. Después de la división, el divisor (el original) se multiplica por el resultado y el producto restado del dividendo. Luego se corrige al sumar o restar el divisor si es necesario (si el resultado de la división se desactivó en uno).
Es fácil dividir un entero de 128 bits por uno de 32 bits, aprovechando el soporte de hardware para la división de 64 bits por 32 bits. En caso de que el divisor <2 ^ 32, uno puede calcular el resto realizando solo 4 divisiones de la siguiente manera:
Supongamos que el dividendo se almacena en:
DWORD dividend[4] = ...
el resto entrará en:
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.
Después de esos 4 pasos, el resto de la variable contendrá lo que estás buscando. (Por favor, no me maten si me equivoqué con la endiarancia. Ni siquiera soy un programador)
En caso de que el divisor tenga más de 2 ^ 32-1, no tengo buenas noticias. No tengo una prueba completa de que el resultado después del cambio esté desactivado por no más de 1, en el procedimiento que describí anteriormente, que creo que MSVC está usando. Sin embargo, creo que tiene algo que ver con el hecho de que la parte que se descarta es al menos 2 ^ 31 veces menor que el divisor, el dividendo es menor que 2 ^ 64 y el divisor es mayor que 2 ^ 32-1 , entonces el resultado es menos de 2 ^ 32.
Si el dividendo tiene 128 bits, el truco de descartar bits no funcionará. Entonces, en general, la mejor solución es probablemente la propuesta por GJ o caf. (Bueno, sería probablemente el mejor incluso si los bits descartados funcionaban. División, resta de multiplicación y corrección en un entero de 128 bits podría ser más lento.)
También estaba pensando en usar el hardware de punto flotante. La unidad de coma flotante x87 usa un formato de precisión de 80 bits con una fracción de 64 bits de longitud. Creo que uno puede obtener el resultado exacto de la división de 64 bits por 64 bits. (No el resto directamente, sino también el resto usando la multiplicación y la resta, como en el "procedimiento MSVC"). SI el dividendo> = 2 ^ 64 y <2 ^ 128 almacenarlo en el formato de coma flotante parece similar a descartar bits menos significativos en el "procedimiento MSVC". Tal vez alguien pueda probar que el error en ese caso está vinculado y lo encuentre útil. No tengo idea de si tiene la oportunidad de ser más rápido que la solución de GJ, pero quizás valga la pena intentarlo.
Puedes usar la versión de división de la Multiplicación campesina rusa .
Para encontrar el resto, ejecute (en pseudo-código):
X = B;
while (X <= A/2)
{
X <<= 1;
}
while (A >= B)
{
if (A >= X)
A -= X;
X >>= 1;
}
El módulo queda en A.
Tendrá que implementar los cambios, las comparaciones y las restas para operar en valores formados por un par de números de 64 bits, pero eso es bastante trivial.
Esto recorrerá como máximo 255 veces (con un A de 128 bits). Por supuesto, debe hacer una comprobación previa para un divisor cero.
Quizás esté buscando un programa terminado, pero los algoritmos básicos para la aritmética de precisión múltiple se pueden encontrar en el Arte de la programación de computadoras de Knuth, Volumen 2. Puede encontrar el algoritmo de división que se describe en línea here . Los algoritmos se ocupan de la aritmética arbitraria de precisión múltiple, por lo que son más generales de lo que necesita, pero debería poder simplificarlos para la aritmética de 128 bits hecha en dígitos de 64 o 32 bits. Esté preparado para una cantidad razonable de trabajo (a) entendiendo el algoritmo, y (b) convirtiéndolo en C o ensamblador.
También es posible que desee comprobar Hacker''s Delight , que está lleno de ensamblador muy inteligente y otros hackers de bajo nivel, incluida una aritmética de precisión múltiple.
Sé que la pregunta especificaba el código de 32 bits, pero la respuesta para 64 bits puede ser útil o interesante para otros.
Y sí, la división 64b / 32b => 32b hace un bloque de construcción útil para 128b% 64b => 64b. El __umoddi3 de __umoddi3 (fuente vinculada a continuación) da una idea de cómo hacer ese tipo de cosas, pero solo implementa 2N% 2N => 2N en la parte superior de una división 2N / N => N, no 4N% 2N => 2N.
Se encuentran disponibles bibliotecas más amplias de precisión múltiple, por ejemplo, https://gmplib.org/manual/Integer-Division.html#Integer-Division .
GNU C en máquinas de 64 bits proporciona un tipo __int128 , y funciones libgcc para multiplicar y dividir de la manera más eficiente posible en la arquitectura de destino.
La instrucción div r/m64 hace 128b / 64b => 64b división (también produce resto como una segunda salida), pero falla si el cociente se desborda. Por lo tanto, no puede usarlo directamente si A/B > 2^64-1 , pero puede hacer que gcc lo use (o incluso en línea el mismo código que usa libgcc).
Esto compila ( Godbolt compiler explorer ) a una o dos instrucciones div (que suceden dentro de una libgcc función libgcc ). Si hubiera una manera más rápida, libgcc probablemente usaría eso en su lugar.
#include <stdint.h>
uint64_t AmodB(unsigned __int128 A, uint64_t B) {
return A % B;
}
La función __umodti3 que llama calcula un módulo completo de 128b / 128b, pero la implementación de esa función verifica el caso especial donde la mitad alta del divisor es 0, como se puede ver en la fuente libgcc . (libgcc crea la versión si / di / ti de la función a partir de ese código, según corresponda para la arquitectura de destino. udiv_qrnnd es una macro asm en línea que realiza la división 2N / N => N sin firmar para la arquitectura de destino.
Para x86-64 (y otras arquitecturas con instrucciones de dividir hardware), el camino rápido (cuando high_half(A) < B ; garantizando que div no tendrá fallas) es solo dos ramas no tomadas , algunas fluff para out-of- ordene el procesamiento de las CPU, y una única instrucción div r64 , que toma alrededor de 50-100 ciclos en las CPU x86 modernas, de acuerdo con agner.org/optimize . Algún otro trabajo puede estar sucediendo en paralelo con div , pero la unidad de división de enteros no está muy canalizada y div decodifica a muchos uops (a diferencia de la división FP).
El camino de retorno todavía usa solo dos instrucciones div 64 bits para el caso en que B es solo de 64 bits, pero A/B no cabe en 64 bits, por lo que A/B directamente fallaría.
Tenga en cuenta que __umodti3 de __umodti3 simplemente inserta __udivmoddi4 en un contenedor que solo devuelve el resto.
Para módulo repetido por el mismo B
Puede valer la pena considerar el cálculo de un inverso multiplicativo de punto fijo para B , si existe. Por ejemplo, con las constantes de tiempo de compilación, gcc hace la optimización para tipos más estrechos 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
La instrucción mul r64 x86 hace 64b * 64b => 128b (rdx: rax) multiplicación, y se puede usar como un bloque de construcción para construir un 128b * 128b => 256b multiplicar para implementar el mismo algoritmo. Como solo necesitamos la mitad alta del resultado completo de 256b, eso ahorra algunas multiplicaciones.
Las CPUs Intel modernas tienen mul rendimiento muy alto: latencia 3c, una por rendimiento de reloj. Sin embargo, la combinación exacta de los cambios y las adiciones requeridas varía con la constante, por lo que el caso general de calcular un inverso multiplicativo en el tiempo de ejecución no es tan eficiente cada vez que se usa como una versión compilada JIT o estáticamente compilada (incluso además de la sobrecarga de precálculo).
IDK donde estaría el punto de equilibrio. Para la compilación de JIT, será superior a ~ 200 reutilizaciones, a menos que guarde en caché el código generado para B valores B comúnmente utilizados. Para la forma "normal", posiblemente esté en el rango de 200 reutilizaciones, pero IDK lo caro que sería encontrar un inverso multiplicativo modular para la división de 128 bits / 64 bits.
libdivide puede hacer esto por usted, pero solo para tipos de 32 y 64 bits. Aún así, es probablemente un buen punto de partida.
If 128-bit unsigned by 63-bit unsigned is good enough, then it can be done in a loop doing at most 63 cycles.
Consider this a proposed solution MSNs'' overflow problem by limiting it to 1-bit. We do so by splitting the problem in 2, modular multiplication and adding the results at the end.
In the following example upper corresponds to the most significant 64-bits, lower to the least significant 64-bits and div is the divisor.
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;
}
The only problem is that, if the divisor is 64-bits then we get overflows of 1-bit (loss of information) giving a faulty result.
It bugs me that I haven''t figured out a neat way to handle the overflows.
If you have a recent x86 machine, there are 128-bit registers for SSE2+. I''ve never tried to write assembly for anything other than basic x86, but I suspect there are some guides out there.
Since there is no predefined 128-bit integer type in C, bits of A have to be represented in an array. Although B (64-bit integer) can be stored in an unsigned long long int variable, it is needed to put bits of B into another array in order to work on A and B efficiently.
After that, B is incremented as Bx2, Bx3, Bx4, ... until it is the greatest B less than A. And then (AB) can be calculated, using some subtraction knowledge for base 2.
Is this the kind of solution that you are looking for?