Maneira mais rápida de calcular um módulo inteiro de 128 bits um inteiro de 64 bits

Eu tenho um inteiro não assinado de 128 bits A e um inteiro sem sinal de 64 bits B. Qual é o caminho mais rápido para calcular A % B – que é o (64 bits) restante da divisão de A por B?

Eu estou olhando para fazer isso em linguagem C ou assembly, mas eu preciso direcionar a plataforma x86 de 32 bits. Isso infelizmente significa que não posso aproveitar o suporte do compilador para números inteiros de 128 bits nem a capacidade da arquitetura x64 de executar a operação necessária em uma única instrução.

Editar:

Obrigado pelas respostas até agora. No entanto, parece-me que os algoritmos sugeridos seriam bastante lentos – a maneira mais rápida de executar uma divisão de 128 bits por 64 bits não seria aproveitar o suporte nativo do processador para a divisão de 32 bits por 64 bits? Alguém sabe se existe uma maneira de realizar a divisão maior em termos de algumas divisões menores?

Re: Quantas vezes B muda?

Primeiramente estou interessado em uma solução geral – qual cálculo você executaria se A e B fossem diferentes a cada vez?

No entanto, uma segunda situação possível é que B não varia tão freqüentemente quanto A – pode haver até 200 Quanto a dividir por B. Como sua resposta diferiria nesse caso?

Você pode usar a versão de divisão do Russian Peasant Multiplication .

Para encontrar o restante, execute (no pseudo-código):

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

O módulo é deixado em A.

Você precisará implementar os turnos, comparações e subtrações para operar em valores compostos de um par de números de 64 bits, mas isso é bastante trivial.

Isso fará um loop de no máximo 255 vezes (com um A de 128 bits). Claro que você precisa fazer uma pré-verificação para um divisor zero.

Talvez você esteja procurando por um programa acabado, mas os algoritmos básicos para a aritmética multi-precisão podem ser encontrados em Art of Computer Programming , Volume 2 de Knuth. Você pode encontrar o algoritmo de divisão descrito on-line aqui . Os algoritmos lidam com aritmética arbitrária de multi-precisão e, portanto, são mais gerais do que você precisa, mas você deve ser capaz de simplificá-los para aritmética de 128 bits feita em dígitos de 64 ou 32 bits. Esteja preparado para uma quantidade razoável de trabalho (a) compreendendo o algoritmo e (b) convertendo-o para C ou assembler.

Você também pode querer dar uma olhada no Hacker’s Delight , que é cheio de muito inteligente assembler e outras hackers de baixo nível, incluindo alguma aritmética multi-precisão.

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 

Se o seu compilador suporta números inteiros de 64 bits, então este é provavelmente o caminho mais fácil. A implementação do MSVC de um módulo de 64 bits em x86 de 32 bits é um assembly cheio de loops peludos ( VC\crt\src\intel\llrem.asm para os bravos), então eu pessoalmente vou com isso.

Isso é quase não testado em parte modificada velocidade mod128by64 function de algoritmo ‘camponês russo’. Infelizmente eu sou um usuário Delphi, então esta function funciona sob o Delphi. 🙂 Mas o montador é quase o mesmo assim …

 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; 

Pelo menos mais uma otimização de velocidade é possível! Depois de ‘Huge Divisor Numbers Shift Optimization’ podemos testar divisores de bit alto, se for 0, não precisamos usar o registro bh extra como o 65º bit para armazenar nele. Então parte desenrolada do loop pode se parecer com:

  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: 

Eu gostaria de compartilhar alguns pensamentos.

Não é tão simples quanto o MSN propõe, eu tenho medo.

Na expressão:

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

tanto multiplicação quanto adição podem transbordar. Eu acho que alguém poderia levar isso em conta e ainda usar o conceito geral com algumas modificações, mas algo me diz que vai ficar realmente assustador.

Eu estava curioso para saber como a operação do módulo de 64 bits foi implementada no MSVC e tentei descobrir algo. Eu realmente não sei assembly e tudo o que eu tinha disponível era edição Express, sem a fonte de VC \ crt \ src \ intel \ llrem.asm, mas eu acho que consegui ter uma idéia do que está acontecendo, depois de um pouco de jogar com a saída do depurador e da desassembly. Eu tentei descobrir como o restante é calculado no caso de inteiros positivos e o divisor> = 2 ^ 32. Existe algum código que lida com números negativos, é claro, mas eu não investiguei isso.

Aqui está como eu vejo isso:

Se divisor> = 2 ^ 32, tanto o dividendo quanto o divisor são deslocados para a direita, tanto quanto necessário, para ajustar o divisor em 32 bits. Em outras palavras: se forem necessários n dígitos para escrever o divisor em binário e n> 32, n-32 serão descartados dígitos menos significativos do divisor e do dividendo. Depois disso, a divisão é realizada usando o suporte de hardware para dividir inteiros de 64 bits por 32 bits. O resultado pode estar incorreto, mas acho que pode ser provado que o resultado pode estar no máximo 1. Após a divisão, o divisor (original) é multiplicado pelo resultado e o produto é subtraído do dividendo. Em seguida, é corrigido adicionando ou subtraindo o divisor, se necessário (se o resultado da divisão estiver desligado por um).

É fácil dividir o inteiro de 128 bits pelo suporte de hardware alavancador de 32 bits para a divisão de 64 bits por 32 bits. No caso do divisor <2 ^ 32, pode-se calcular o restante executando apenas 4 divisões da seguinte forma:

Vamos supor que o dividendo seja armazenado em:

 DWORD dividend[4] = ... 

o restante irá para:

 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. 

Após esses 4 passos, o restante variável irá conter o que você está procurando. (Por favor, não me mate se eu tiver a endianess errada. Eu nem sou um programador)

Caso o divisor seja maior que 2 ^ 32-1 não tenho boas notícias. Não tenho uma prova completa de que o resultado após o turno esteja desativado em no máximo 1, no procedimento que descrevi anteriormente, que acredito que o MSVC esteja usando. Eu penso, entretanto, que tem algo a ver com o fato de que a parte que é descartada é pelo menos 2 ^ 31 vezes menor que o divisor, o dividendo é menor que 2 ^ 64 e o divisor é maior que 2 ^ 32-1 , então o resultado é menor que 2 ^ 32.

Se o dividendo tiver 128 bits, o truque para descartar os bits não funcionará. Então, em geral, a melhor solução é provavelmente aquela proposta por GJ ou caf. (Bem, seria provavelmente o melhor mesmo se os bits de descarte funcionassem. Divisão, subtração de multiplicação e correção no inteiro de 128 bits podem ser mais lentos.)

Eu também estava pensando em usar o hardware de ponto flutuante. A unidade de ponto flutuante x87 usa um formato de precisão de 80 bits com fração de 64 bits. Acho que podemos obter o resultado exato da divisão de 64 bits por 64 bits. (Não o restante diretamente, mas também o restante usando multiplicação e subtração como no “procedimento MSVC”). SE o dividendo> = 2 ^ 64 e <2 ^ 128 armazenando-o no formato de ponto flutuante parecer semelhante a descartar bits menos significativos no "procedimento MSVC". Talvez alguém possa provar que o erro nesse caso está vinculado e achar útil. Eu não tenho ideia se tem uma chance de ser mais rápido que a solução do GJ, mas talvez valha a pena tentar.

A solução depende do que exatamente você está tentando resolver.

Por exemplo, se você está fazendo aritmética em um módulo de anel de um inteiro de 64 bits, em seguida, usando a redução de Montgomerys é muito eficiente. É claro que isso pressupõe que você tenha o mesmo módulo várias vezes e que compensa converter os elementos do anel em uma representação especial.


Para dar apenas uma estimativa muito aproximada sobre a velocidade desta redução de Montgomerys: Eu tenho um benchmark antigo que executa uma exponenciação modular com módulo de 64 bits e expoente em 1600 ns em um Núcleo de 2.4Ghz 2. Esta exponenciação faz cerca de 96 multiplicações modulares ( e reduções modulares) e, portanto, precisa de cerca de 40 ciclos por multiplicação modular.

Eu fiz ambas as versões da function de divisão ‘Camponês Russo’ Mod128by64: clássico e velocidade otimizada. A velocidade otimizada pode fazer no meu PC 3Ghz mais de 1000.000 cálculos randoms por segundo e é mais de três vezes mais rápida do que a function clássica. Se compararmos o tempo de execução do cálculo de 128 por 64 e calcular o módulo de 64 por 64 bits, essa function é apenas cerca de 50% mais lenta.

Camponês clássico russo:

 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; 

Velocidade otimizada camponesa russa:

 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; 

A resposta aceita pelo @caf foi muito boa e altamente cotada, mas contém um bug que não é visto há anos.

Para ajudar a testar essa e outras soluções, estou postando um equipamento de teste e tornando-o wiki da comunidade.

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

Eu sei que a pergunta especificada código de 32 bits, mas a resposta para 64 bits pode ser útil ou interessante para os outros.

E sim, a divisão 64b / 32b => 32b faz um bloco de construção útil para 128b% 64b => 64b. O __umoddi3 da libgcc (fonte ligada abaixo) dá uma idéia de como fazer esse tipo de coisa, mas implementa apenas 2N% 2N => 2N em cima de uma divisão 2N / N => N, não 4N% 2N => 2N.

Bibliotecas de multiprecisão mais abrangentes estão disponíveis, por exemplo, https://gmplib.org/manual/Integer-Division.html#Integer-Division .


O GNU C em máquinas de 64 bits fornece um tipo __int128 e libgcc para multiplicar e dividir da maneira mais eficiente possível na arquitetura de destino.

A instrução div r/m64 x86-64 faz a divisão 128b / 64b => 64b (também produzindo o restante como uma segunda saída), mas falha se o quociente estourar. Então você não pode usá-lo diretamente se A/B > 2^64-1 , mas você pode obter o gcc para usá-lo para você (ou mesmo inline o mesmo código que libgcc usa).


Isto compila ( explorador do compilador Godbolt ) para uma ou duas instruções div (que acontecem dentro de uma chamada de function libgcc ). Se houvesse uma maneira mais rápida, a libgcc provavelmente usaria isso.

 #include  uint64_t AmodB(unsigned __int128 A, uint64_t B) { return A % B; } 

A function __umodti3 ele chama calcula um módulo 128b / 128b completo, mas a implementação dessa function verifica o caso especial em que a metade alta do divisor é 0, como você pode ver na fonte libgcc . (O libgcc constrói a versão si / di / ti da function a partir desse código, conforme apropriado para a arquitetura de destino. udiv_qrnnd é uma macro asm inline que faz uma divisão não-programada 2N / N => N para a arquitetura de destino.

Para x86-64 (e outras arquiteturas com uma instrução de divisão de hardware), o caminho rápido (quando high_half(A) < B ; garantindo div não high_half(A) < B ) é apenas duas ramificações não tomadas , alguma penugem para fora-de- pedidos de CPU para mastigar, e uma única instrução div r64 , que leva cerca de 50-100 ciclos em CPUs modernas x86, de acordo com as tabelas de insn de Agner Fog . Algum outro trabalho pode estar acontecendo em paralelo com div , mas a unidade de divisão inteira não é muito pipeline e div decodifica para muitos uops (ao contrário da divisão FP).

O caminho de fallback ainda usa apenas duas instruções div 64 bits para o caso em que B é de apenas 64 bits, mas A/B não se encheckbox em 64 bits, portanto, A/B falharia diretamente.

Note que __umodti3 da libgcc apenas incorpora __udivmoddi4 em um wrapper que apenas retorna o resto.


Para módulo repetido pelo mesmo B

Pode valer a pena considerar o cálculo de um inverso multiplicativo de ponto fixo para B , se existir um. Por exemplo, com constantes de tempo de compilation, o gcc faz a otimização para tipos menores 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 

A instrução mul r64 do x86 mul r64 64b * 64b => 128b (rdx: rax), e pode ser usada como um bloco de construção para construir um 128b * 128b => 256b para implementar o mesmo algoritmo. Como só precisamos da metade alta do resultado total de 256b, isso economiza algumas multiplicações.

Os modernos processadores da Intel têm um desempenho muito alto: 3c latência, um por clock. No entanto, a combinação exata de turnos e acréscimos necessários varia com a constante, portanto, o caso geral de calcular um inverso multiplicativo em tempo de execução não é tão eficiente cada vez que é usado como uma versão JIT-compilada ou estaticamente compilada (mesmo na parte superior da sobrecarga de pré-cálculo).

IDK onde o ponto de equilíbrio seria. Para compilation JIT, ele será maior que ~ 200 reutilizações, a menos que você armazene em cache código gerado para valores B comumente usados. Para o modo "normal", pode ser que esteja no intervalo de 200 reutilizações, mas IDK como seria caro encontrar um inverso multiplicativo modular para a divisão de 128 bits / 64 bits.

A libdivide pode fazer isso para você, mas apenas para os tipos 32 e 64 bits. Ainda assim, é provavelmente um bom ponto de partida.

Como regra geral, a divisão é lenta e a multiplicação é mais rápida, e a mudança de bits é mais rápida ainda. Pelo que vi das respostas até agora, a maioria das respostas tem usado uma abordagem de força bruta usando turnos de bits. Existe outro caminho. Se é mais rápido continua a ser visto (perfil AKA).

Em vez de dividir, multiplique pelo recíproco. Assim, para descobrir A% B, calcule primeiro o recíproco de B … 1 / B. Isso pode ser feito com alguns loops usando o método de convergência de Newton-Raphson. Para fazer isso bem, dependerá de um bom conjunto de valores iniciais em uma tabela.

Para mais detalhes sobre o método de Newton-Raphson de convergir para o recíproco, por favor consulte http://en.wikipedia.org/wiki/Division_(digital)

Depois de ter o recíproco, o quociente Q = A * 1 / B.

O restante R = A – Q * B.

Para determinar se isso seria mais rápido que a força bruta (já que haverá muito mais multiplicações, já que estaremos usando registros de 32 bits para simular números de 64 e 128 bits, faça o perfil dele).

Se B é constante em seu código, você pode pré-calcular o recíproco e simplesmente calcular usando as duas últimas fórmulas. This, I am sure will be faster than bit-shifting.

Espero que isto ajude.

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?