最快的方法来计算一个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?