La forma más rápida de calcular un módulo entero de 128 bits un entero de 64 bits

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 rest (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 comstackdor para enteros de 128 bits, ni la capacidad de la architecture 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?

Puedes usar la versión de división de la Multiplicación campesina rusa .

Para encontrar el rest, 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 progtwig terminado, pero los algoritmos básicos para la aritmética de precisión múltiple se pueden encontrar en el Arte de la progtwigción de computadoras de Knuth, Volumen 2. Puede encontrar el algoritmo de división que se describe en línea aquí . 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.

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 comstackdor 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: 

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 sum 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 rest 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 sumr 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 rest 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 rest de la variable contendrá lo que estás buscando. (Por favor, no me maten si me equivoqué con la endiarancia. Ni siquiera soy un progtwigdor)

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 rest directamente, sino también el rest 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.

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.

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

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 architecture de destino.

La instrucción div r/m64 hace 128b / 64b => 64b división (también produce rest 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 comstack ( Godbolt compiler explorer ) a una o dos instrucciones div (que suceden dentro de una llamada de función libgcc ). Si hubiera una manera más rápida, libgcc probablemente usaría eso en su lugar.

 #include  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 architecture de destino. udiv_qrnnd es una macro asm en línea que realiza la división 2N / N => N sin firmar para la architecture de destino.

Para x86-64 (y otras architectures con instrucciones de dividir hardware), el camino rápido (cuando high_half(A) < B ; garantizando que div no tendrá fallas) es solo dos twigs 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 las tablas internas de Agner Fog . 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 rest.


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 comstackció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 comstackda JIT o estáticamente comstackda (incluso además de la sobrecarga de precálculo).

IDK donde estaría el punto de equilibrio. Para la comstackció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.

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)

Una vez que tiene el recíproco, el cociente Q = A * 1 / B.

El rest R = A – Q * B.

Para determinar si esto sería más rápido que la fuerza bruta (ya que habrá muchas más multiplicaciones ya que utilizaremos registros de 32 bits para simular números de 64 y 128 bits, perfilarlo).

Si B es constante en su código, puede precalcular el recíproco y simplemente calcular usando las dos últimas fórmulas. Esto, estoy seguro, será más rápido que el cambio de bit.

Espero que esto ayude.

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.

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.

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?