use*_*783 53 c algorithm x86 assembly modulo
我有一个128位无符号整数A和一个64位无符号整数B.什么是最快的计算方法A % B- 即将A除以B的(64位)余数？
我希望用C或汇编语言来做这件事,但我需要针对32位x86平台.遗憾的是,我无法利用编译器对128位整数的支持,也无法利用x64架构在单条指令中执行所需操作的能力.
编辑:
谢谢你到目前为止的答案.但是,在我看来,建议的算法会非常慢 - 执行128位到64位除法的最快方法是利用处理器对64位乘32位除法的原生支持吗？有没有人知道是否有办法在一些较小的部门中执行更大的划分？
回复:B多久换一次？
主要是我对一般解决方案感兴趣 - 如果A和B每次都可能不同,你会进行什么计算？
然而,第二种可能的情况是B不会像A那样经常变化 - 每个B可能有多达200个As除以.在这种情况下,你的答案有何不同？
caf*_*caf 31
您可以使用俄罗斯农民增殖的分部版本.
要查找余数,请执行(伪代码):
X = B;
while (X <= A/2)
{
X <<= 1;
}
while (A >= B)
{
if (A >= X)
A -= X;
X >>= 1;
}
模数保留在A.
你需要实现移位,比较和减法来操作由一对64位数组成的值,但这是相当微不足道的.
这将循环最多255次(128位A).当然,您需要对零除数进行预检.
Dal*_*und 13
也许你正在寻找一个完成的程序,但是可以在Knuth的计算机编程艺术第2卷中找到多精度算术的基本算法.你可以在这里找到在线描述的除法算法.算法处理任意多精度算术,因此比您需要的更通用,但您应该能够在64位或32位数字上对128位算术进行简化.准备好合理的工作量(a)理解算法,(b)将其转换为C或汇编程序.
您可能还想查看Hacker's Delight,它充满了非常聪明的汇编程序和其他低级别的hackery,包括一些多精度算法.
MSN*_*MSN 11
鉴于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对于勇敢的人),所以我个人会这样做.
这几乎是未经测试的部分速度修改的Mod128by64'俄罗斯农民'算法功能.不幸的是我是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;
至少还有一种速度优化是可能的!在"巨大除数数移位优化"之后,我们可以测试除数高位,如果为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:
我知道问题指定了 32 位代码，但 64 位的答案可能对其他人有用或有趣。
是的，64b/32b => 32b 除法确实是 128b % 64b => 64b 的有用构建块。libgcc 的__umoddi3（源链接如下）给出了如何做那种事情的想法，但它只在 2N / N => N 除法之上实现了 2N % 2N => 2N，而不是 4N % 2N => 2N。
可以使用更广泛的多精度库，例如https://gmplib.org/manual/Integer-Division.html#Integer-Division。
64 位机器上的 GNU C确实提供了一个__int128type和 libgcc 函数来在目标架构上尽可能有效地进行乘法和除法。
x86-64 的div r/m64指令执行 128b/64b => 64b 除法（也产生余数作为第二个输出），但如果商溢出，它会出错。所以你不能直接使用它 if A/B > 2^64-1，但你可以让 gcc 为你使用它（甚至内联 libgcc 使用的相同代码）。
这将编译（Godbolt 编译器资源管理器）成一两条div指令（发生在libgcc函数调用中）。如果有更快的方法，libgcc 可能会使用它。
#include <stdint.h>
uint64_t AmodB(unsigned __int128 A, uint64_t B) {
return A % B;
}
__umodti3它调用的函数计算完整的 128b/128b 模，但该函数的实现确实检查除数的高半部分为 0 的特殊情况，如您在 libgcc 源代码中所见。（libgcc 从该代码构建函数的 si/di/ti 版本，适用于目标架构。udiv_qrnnd是一个内联 asm 宏，它为目标架构 执行无符号 2N/N => N 除法。
对于 x86-64（以及其他具有硬件除法指令的体系结构），快速路径（当high_half(A) < B; 保证div不会出错时）只是两个未采用的分支，一些用于无序 CPU 咀嚼的绒毛，以及div r64根据Agner Fog 的 insn 表，单条指令在现代 x86 CPU 上大约需要 50-100 个周期1。其他一些工作可以与 并行发生div，但整数除法单元不是非常流水线化并且div解码为很多 uops（与 FP 除法不同）。
对于只有 64 位div的情况，回退路径仍然只使用两个 64 位指令B，但A/B不适合 64 位，因此A/B直接会出错。
请注意，libgcc__umodti3只是内联__udivmoddi4到仅返回余数的包装器中。
脚注 1：32 位div在 Intel CPU 上快 2 倍以上。在 AMD CPU 上，性能仅取决于实际输入值的大小，即使它们是 64 位寄存器中的小值。如果小值很常见，那么在进行 64 位或 128 位除法之前，可能值得将分支基准测试为简单的 32 位除法版本。
B这可能是值得考虑的计算定点乘法逆的B，如果存在。例如，对于编译时常量，gcc 会针对小于 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
x86 的mul r64指令做 64b*64b => 128b (rdx:rax) 乘法，并且可以用作构建块来构造 128b * 128b => 256b 乘法来实现相同的算法。由于我们只需要完整 256b 结果的高一半，因此可以节省一些乘法。
现代英特尔 CPU 具有非常高的性能mul：3c 延迟，每个时钟吞吐量一个。但是，所需的移位和加法的确切组合因常数而异，因此在运行时计算乘法逆的一般情况在每次用作 JIT 编译或静态编译版本（甚至在预计算开销之上）。
IDK 盈亏平衡点所在。对于 JIT 编译，它会高于约 200 次重用，除非您为常用B值缓存生成的代码。对于“正常”方式，它可能在 200 次重用的范围内，但 IDK 为 128 位 / 64 位除法找到模块化乘法逆是多么昂贵。
libdivide可以为您执行此操作，但仅适用于 32 位和 64 位类型。不过，这可能是一个很好的起点。
| 归档时间： |
|
| 查看次数： |
15876 次 |
| 最近记录： |