我有一个 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 来除。在这种情况下，你的答案会有什么不同？
您可以使用俄罗斯农民乘法的除法版本。
要找到余数，请执行（在伪代码中）：
X = B;
while (X <= A/2)
{
X <<= 1;
}
while (A >= B)
{
if (A >= X)
A -= X;
X >>= 1;
}
模数留在 A 中。
您需要实现移位、比较和减法以对由一对 64 位数字组成的值进行运算，但这相当简单（可能您应该将左移 1 实现为X + X）。
这将最多循环 255 次（使用 128 位 A）。当然，您需要对零除数进行预检查。
也许您正在寻找一个完成的程序，但多精度算术的基本算法可以在 Knuth 的计算机编程艺术，第 2 卷中找到。您可以在此处找到在线描述的除法算法。这些算法处理任意多精度算术，因此比您需要的更通用，但您应该能够将它们简化为在 64 位或 32 位数字上完成的 128 位算术。为合理的工作量做好准备 (a) 理解算法，以及 (b) 将其转换为 C 或汇编程序。
您可能还想查看Hacker's Delight，它充满了非常聪明的汇编程序和其他低级黑客技术，包括一些多精度算术。
如果您的 B 足够小以使uint64_t +操作不换行：
给定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解码为很多微指令（与 FP 除法不同）。
对于只有 64 位但不适合 64 位div的情况，回退路径仍然只使用两条 64 位指令，因此直接会出错。BA/BA/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 位类型。不过，这可能是一个很好的起点。
解决方案取决于您要解决的问题。
例如，如果您在以 64 位整数为模的环中进行算术运算，那么使用 Montgomerys 约简是非常有效的。当然，这假设您多次使用相同的模数，并且将环的元素转换为特殊表示是值得的。
对这个蒙哥马利减少的速度给出一个非常粗略的估计：我有一个旧的基准，它在 2.4Ghz 核心 2 上以 1600 ns 的 1600 ns 执行 64 位模数和指数的模幂运算。这个幂运算大约进行 96 次模乘（和模减少），因此每次模乘需要大约 40 个周期。
我制作了两个版本的 Mod128by64 '俄罗斯农民'除法功能：经典和速度优化。速度优化可以在我的 3Ghz PC 上每秒执行超过 1000.000 次随机计算，比经典功能快三倍以上。如果我们比较计算 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;
我想分享一些想法。
恐怕它不像 MSN 建议的那么简单。
在表达式中：
(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B
乘法和加法都可能溢出。我认为人们可以考虑到这一点，并且仍然使用一般概念并进行一些修改，但有些事情告诉我它会变得非常可怕。
我很好奇 64 位模运算是如何在 MSVC 中实现的，我试图找出答案。我真的不知道汇编，我只有 Express 版本，没有 VC\crt\src\intel\llrem.asm 的来源，但我想我在玩了一会儿之后设法了解了发生了什么带有调试器和反汇编输出。我试图弄清楚在正整数和除数> = 2 ^ 32的情况下如何计算余数。当然有一些处理负数的代码，但我没有深入研究。
这是我的看法：
如果除数 >= 2^32 则被除数和除数都向右移动，以使除数适合 32 位。换句话说：如果需要 n 位以二进制形式写下除数并且 n > 32，则除数和被除数的 n-32 个最低有效位都将被丢弃。之后，使用硬件支持将 64 位整数除以 32 位整数来执行除法。结果可能不正确，但我认为可以证明，结果最多可能相差 1。除法后，除数（原始一个）乘以结果，并从被除数中减去乘积。然后在必要时通过添加或减去除数来纠正它（如果除法的结果被减一）。
利用对 64 位除以 32 位除法的硬件支持，很容易将 128 位整数除以 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 个步骤之后，变量剩余部分将保存您正在寻找的内容。（如果我弄错了字节顺序，请不要杀我。我什至不是程序员）
如果除数大于 2^32-1 我没有好消息。在我之前描述的过程中，我没有完整的证据证明移位后的结果不超过 1，我相信 MSVC 正在使用该过程。然而，我认为这与以下事实有关，即被丢弃的部分至少比除数小 2^31 倍，被除数小于 2^64，除数大于 2^32-1 ，所以结果小于 2^32。
如果被除数有 128 位，那么丢弃位的技巧将不起作用。所以在一般情况下，最好的解决方案可能是 GJ 或 caf 提出的解决方案。（好吧，即使丢弃位有效，它也可能是最好的。128 位整数的除法、乘法减法和校正可能会更慢。）
我也在考虑使用浮点硬件。x87 浮点单元使用 80 位精度格式，长度为 64 位。我认为可以得到 64 位除以 64 位的确切结果。（不是直接的余数，而是使用“MSVC 过程”中的乘法和减法的余数）。如果除数 >=2^64 且 < 2^128 以浮点格式存储它似乎类似于丢弃“MSVC 过程”中的最低有效位。也许有人可以证明这种情况下的错误是必然的，并发现它很有用。我不知道它是否有机会比 GJ 的解决方案更快，但也许值得一试。
@caf 接受的答案非常好，评价很高，但它包含一个多年未见的错误。
为了帮助测试该解决方案和其他解决方案，我发布了一个测试工具并使其成为社区 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;
}
作为一般规则，除法很慢，乘法更快，而位移更快。从到目前为止我所看到的答案来看，大多数答案都是使用位移位的蛮力方法。还有另一种方式。它是否更快还有待观察（AKA profile it）。
不是除法，而是乘以倒数。因此，要发现 A % B，首先计算 B ... 1/B 的倒数。这可以通过使用 Newton-Raphson 收敛方法的几个循环来完成。要做到这一点，将取决于表中一组好的初始值。
有关在倒数上收敛的 Newton-Raphson 方法的更多详细信息，请参阅http://en.wikipedia.org/wiki/Division_(digital)
一旦得到倒数，商 Q = A * 1/B。
余数 R = A - Q*B。
为了确定这是否会比蛮力更快（因为我们将使用 32 位寄存器来模拟 64 位和 128 位数字，所以会有更多的乘法），对其进行分析。
如果 B 在您的代码中是常数，您可以预先计算倒数并使用最后两个公式简单地计算。我相信这会比位移更快。
希望这可以帮助。
如果 128 位无符号乘 63 位无符号足够好，那么它可以在最多 63 个周期的循环中完成。
通过将其限制为 1 位来考虑这是一个建议的解决方案 MSN 的溢出问题。我们通过将问题拆分为 2、模乘并将结果相加来做到这一点。
在以下示例中，upper 对应于最高有效 64 位，lower 对应于最低有效 64 位，div 是除数。
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;
}
唯一的问题是，如果除数是 64 位，那么我们会得到 1 位的溢出（信息丢失），从而导致错误的结果。
让我感到困扰的是，我还没有找到一种巧妙的方法来处理溢出。
我不知道如何编译汇编代码，感谢任何帮助编译和测试它们。
我通过与 gmplib "mpz_mod()" 进行比较并汇总 100 万个循环结果来解决了这个问题。从减速（seedup 0.12）到加速 1.54 需要很长时间——这就是我认为这个线程中的 C 代码会很慢的原因。
此线程中包含测试工具的详细信息：
https ://www.raspberrypi.org/forums/viewtopic.php?f=33&t=311893&p=1873122#p1873122
这是“mod_256()”，比使用 gmplib“mpz_mod()”更快，使用 __builtin_clzll() 进行更长的轮班是必不可少的：
typedef __uint128_t uint256_t[2];
#define min(x, y) ((x<y) ? (x) : (y))
int clz(__uint128_t u)
{
// unsigned long long h = ((unsigned long long *)&u)[1];
unsigned long long h = u >> 64;
return (h!=0) ? __builtin_clzll(h) : 64 + __builtin_clzll(u);
}
__uint128_t mod_256(uint256_t x, __uint128_t n)
{
if (x[1] == 0) return x[0] % n;
else
{
__uint128_t r = x[1] % n;
int F = clz(n);
int R = clz(r);
for(int i=0; i<128; ++i)
{
if (R>F+1)
{
int h = min(R-(F+1), 128-i);
r <<= h; R-=h; i+=(h-1); continue;
}
r <<= 1; if (r >= n) { r -= n; R=clz(r); }
}
r += (x[0] % n); if (r >= n) r -= n;
return r;
}
}
如果您有一台最近的 x86 机器，则 SSE2+ 有 128 位寄存器。除了基本的 x86 之外，我从未尝试为任何东西编写程序集，但我怀疑那里有一些指南。
我在战斗结束 9 年后，但这里有一个有趣的 O(1) 边缘案例，值得一提的是 2 的幂。
#include <stdio.h>
// example with 32 bits and 8 bits.
int main() {
int i = 930;
unsigned char b = (unsigned char) i;
printf("%d", (int) b); // 162, same as 930 % 256
}
由于 C 中没有预定义的 128 位整数类型，因此必须在数组中表示 A 的位。尽管 B（64 位整数）可以存储在unsigned long long int变量中，但需要将 B 的位放入另一个数组中才能有效地处理 A 和 B。
之后，B 递增为 Bx2，Bx3，Bx4，...，直到它是小于 A 的最大 B。然后可以使用一些以 2 为底的减法知识计算 (AB)。
这是您正在寻找的解决方案吗？