最快的方法来计算一个128位整数模64位整数
我有一个128位无符号整数A和一个64位无符号整数B.计算A % B的最快方法是什么？这是将A除以B得到的（64位）余数？
我正在寻找以C或汇编语言来做到这一点，但我需要针对32位的x86平台。 这不幸意味着我不能利用128位整数的编译器支持，也不能利用x64体系结构在单个指令中执行所需操作的能力。
编辑：
感谢您迄今的答复。 但是，在我看来，推荐的algorithm会非常慢 – 不是执行128位乘64位除法的最快方法是利用处理器对64位乘32位的本机支持吗？ 有没有人知道是否有办法在几个较小的部门执行更大的部门？
Re：B多久换一次？
主要是我对一个通用的解决scheme感兴趣 – 如果A和B每次都可能会有所不同，你会进行什么样的计算？
但是，第二种可能的情况是B不会像A那么频繁地变化 – 可能有多达200个A被B分开。在这种情况下你的答案会有什么不同？
你可以使用俄文农民增殖版本。
要find余数，执行（在伪代码中）：
X = B; while (X <= A/2) { X <<= 1; } while (A >= B) { if (A >= X) A -= X; X >>= 1; }
模量留在A.
您需要执行移位，比较和减法操作，由一对64位数组成的值进行操作，但这相当简单。
这将最多循环255次（128位A）。 当然，你需要做一个零除数的预检。
也许你正在寻找一个完成的程序，但是多精度algorithm的基本algorithm可以在Knuth的“计算机编程艺术”第2卷中find。你可以在这里find在线描述的除法algorithm。 这些algorithm处理任意的多精度算术，所以比你需要的更一般化，但是你应该能够简化它们在64位或32位数字上完成的128位算术。 准备合理的工作量（a）理解algorithm，（b）将其转换为C或汇编程序。
您可能还想看看Hacker's Delight ，这个Hacker's Delight有很多非常聪明的汇编和其他低级hackery，包括一些多精度algorithm。
给定A = AH*2^64 + AL ：
A % B == (((AH % B) * (2^64 % B)) + (AL % B)) % B == (((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
如果你的编译器支持64位整数，那么这可能是最简单的方法。 MSVC在32位x86上的64位模数的实现是一些多毛的循环填充程序集（勇敢的VC\crt\src\intel\llrem.asm ），所以我个人认为是这样的。
这几乎是未经testing的部分速度修改Mod128by64'俄罗斯农民'algorithm的function。 不幸的是我是一个Delphi用户，所以这个函数在Delphi下工作。 :)但汇编程序几乎是一样的，所以…
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;
至less还有一个速度优化是可能的！ 在'巨大除数转换优化'后，我们可以testing除数高位，如果它是0，我们不需要使用额外的bh寄存器作为第65位存储在它。 所以循环的展开部分可以看起来像：
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:
我想分享一些想法。
这恐怕不像MSN提出的那样简单。
在expression式中：
(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
乘法和加法都可能溢出。 我认为可以考虑到这一点，仍然使用一般的概念，但有一些事情告诉我它会变得非常可怕。
我很好奇64位模数操作是如何在MSVC中实现的，我试图找出一些东西。 我真的不知道程序集，所有我可用的是Express版本，没有VC \ crt \ src \ intel \ llrem.asm的源代码，但是我想我在devise一些游戏后得到了一些想法与debugging器和反汇编输出。 我试图找出余数是如何计算的正整数和除数> = 2 ^ 32。 有一些代码负责处理当然，但我没有深入了解。
以下是我的看法：
如果除数> = 2 ^ 32，那么除法和除数都被正确地移位以使除数符合32位。 换句话说，如果需要n位数字以二进制forms写入除数，且n> 32，则除数和除数的n-32个最低有效位数将被丢弃。 之后，使用硬件支持将64位整数除以32位进行除法。 结果可能是不正确的，但我认为可以certificate，结果最多可以是1。除法之后，除数（原来的）乘以结果，并从股息中扣除乘积。 然后根据需要加上或减去除数来校正（如果除法结果为1）。
很容易将128位整数除以32位，利用64位乘32位除法的硬件支持。 在除数<2 ^ 32的情况下，可以如下计算余数只执行4个分割：
我们假设股息存储在：
DWORD dividend[4] = ...
其余的将进入：
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.
经过这4个步骤后，variables余数将保留你正在寻找的东西。 （如果我弄错了，请不要杀我，我甚至不是程序员）
如果除数大于2 ^ 32-1，我没有好消息。 我没有一个完全的证据certificate，在我之前描述的我认为MSVC正在使用的程序中，转换后的结果不超过1。 但我认为这与被抛弃的部分至less比除数less2 ^ 31倍有关，被除数小于2 ^ 64，除数大于2 ^ 32-1 ，结果小于2 ^ 32。
如果红利有128位，丢弃位的技巧将不起作用。 所以在一般情况下，最好的解决办法可能是GJ或caf提出的scheme。 （即使丢弃位可能是最好的，128位整数的乘法减法和校正可能会变慢）。
我也在考虑使用浮点硬件。 x87浮点单元使用长度为64位的80位精度格式。 我认为可以得到64位除以64位的确切结果。 （不是直接的余数，而是“MSVC过程”中使用乘法和减法的余数）。 如果以“浮点格式”存储的除数> = 2 ^ 64和<2 ^ 128似乎与丢弃“MSVC过程”中的最低有效位相似。 也许有人可以certificate在这种情况下的错误是有约束力的，并find它的用处。 我不知道它是否有机会比GJ的解决scheme更快，但也许值得尝试。
解决scheme取决于你正在试图解决什么。
例如，如果你正在做一个64位整数环模算术，那么使用Montgomer减less是非常有效的。 当然，这假定你多次使用相同的模数，并且将环的元素转换成特殊的表示。
对Montgomer减less的速度给出一个非常粗略的估计：我有一个旧的基准，它在2.4Ghz核心2上执行64位模数和1600 ns的指数的模幂运算。这个指数约有96个模乘（和模块化减less），因此每个模块化乘法需要大约40个周期。
我做了两个版本的Mod128by64'俄罗斯农民'分工function：经典和速度优化。 速度优化可以在我的3Ghz PC上进行超过每秒1000.000次的随机计算，比传统function快三倍以上。 如果我们将计算128的执行时间乘以64，并用64位的模数来计算64，那么比这个函数慢了大约50％。
经典的俄罗斯农民：
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;
速度优化的俄罗斯农民：
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;
@caf接受的答案真的很好，评价很高，但是它包含了一个多年没见的bug。
为了帮助testing这个和其他解决scheme，我发布了一个testing工具并使其成为社区wiki。
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; }
我知道指定32位代码的问题，但64位的答案可能对别人有用或有趣。
是的，64b / 32b => 32b划分确实为128b％64b => 64b构build了一个有用的构build块。 libgcc的__umoddi3 （源链接如下）给出了如何做这种事情的想法，但是它只在2N / N => N分区的顶部实现了2N％2N => 2N，而不是4N％2N => 2N。
可以使用更宽的多精度库，例如https://gmplib.org/manual/Integer-Division.html#Integer-Division 。
64位机器上的GNU C确实提供了__int128types和libgcc函数，以便在目标架构上尽可能高效地进行乘法和除法运算。
x86-64的div r/m64指令执行128b / 64b => 64b除法（也产生余数作为第二个输出），但是如果商溢出则会发生故障。 所以如果A/B > 2^64-1 ，你不能直接使用它，但是你可以得到gcc来为你使用它（或者甚至内联与libgcc使用的相同的代码）。
这编译（ Godbolt编译器资源pipe理器 ）一个或两个div指令（发生在libgcc函数调用内）。 如果有更快的方法，libgcc可能会使用它。
#include <stdint.h> uint64_t AmodB(unsigned __int128 A, uint64_t B) { return A % B; }
它所调用的__umodti3函数计算完整的128b / 128b模数，但是该函数的实现确实检查了除数的高一半是0的特殊情况，就像在libgcc源代码中看到的一样 。 （libgcc从该代码构build函数的si / di / ti版本，适合于目标体系结构udiv_qrnnd是一个内联asmmacros，对目标体系结构执行无符号2N / N => N分割。
对于x86-64 （以及其他具有硬件除法指令的体系结构）， 快速path （当high_half(A) < B ;保证div不会错误） 只是两个未被采用的分支 ，根据Agner Fog's insn表格 ，在现代x86 CPU上需要大约50-100个周期的单个div r64指令 。 一些其他的工作可以与div并行发生，但整数除法单元不是很stream水线， div解码到很多的uops（不像FP分区）。
在B只有64位的情况下，回退path仍然只使用两个64位的div指令，但A/B不适合64位，所以A/B直接发生故障。
请注意，libgcc的__umodti3只是将__umodti3联到一个只返回余数的包装器中。
由同一个B重复模
如果存在的话，可以考虑计算B的模乘法逆 。 例如，使用编译时常量，gcc对比128b更窄的types进行优化。
uint64_t modulo_by_constant64(uint64_t A) { return A % 0x12345678ABULL; } movabs rdx, -2233785418547900415 mov rax, rdi mul rdx mov rax, 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指令执行64b * 64b => 128b（rdx：rax）乘法，可以用作构build128b * 128b => 256b乘法的构build块来实现相同的algorithm。 由于我们只需要完整256b结果的高一半，所以可以节省几个倍数。
现代英特尔CPU具有非常高的性能，每个时钟的吞吐量为1个。 然而，所需的移位和相加的确切组合要随着常量而变化，所以在运行时计算乘法倒数的一般情况在每次用作JIT编译或静态编译版本（偶数在预计算开销之上）。
IDK的盈亏平衡点。 对于JIT编译，除非您为常用的B值caching生成的代码，否则它将高于〜200次重用。 对于“正常”的方式，它有可能在200次重用的范围内，但是IDK在128位/ 64位除法中find一个模乘的逆是多么昂贵。
libdivide可以为你做这个，但只适用于32位和64位的types。 不过，这可能是一个很好的起点。
一般来说，划分速度慢，乘法速度快，移位速度快。 从迄今为止我所见到的答案来看，大部分的答案都是使用位移的暴力方法。 还有另一种方式。 是否更快还有待观察（AKA简介）。
而不是分裂，乘以倒数。 因此，要发现A％B，首先计算B … 1 / B的倒数。 这可以通过使用牛顿 – 拉夫逊收敛方法的几个循环来完成。 要做到这一点，将取决于表中的一组好的初始值。
关于牛顿 – 拉夫逊方法在互逆上的更多细节，请参考http://en.wikipedia.org/wiki/Division_(digital）;
一旦你有倒数，商Q = A * 1 / B。
余数R = A-Q * B。
为了确定这是否会比蛮力更快（因为将会有更多的乘法，因为我们将使用32位寄存器来模拟64位和128位数字，请对其进行分析。
如果在代码中B是常量，则可以预先计算倒数，并使用最后两个公式进行简单计算。 这一点，我相信会比移位更快。
希望这可以帮助。
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?
