检查 64 位常量乘法参数
Checking parameters of multiplication by constant in 64 bit
对于我的 BigInteger 代码，对于非常大的 BigInteger，输出结果很慢。所以现在我使用递归的分而治之算法，仍然需要2'30"将当前已知最大的素数转换为超过2200万位的十进制字符串（但只需135毫秒即可将其转换为十六进制字符串） .
我仍然想减少时间，所以我需要一个可以非常快速地将 NativeUInt（即 32 位平台上的 UInt32，64 位平台上的 UInt64）除以 100 的例程。所以我用常量乘法。这在 32 位代码中工作正常，但我不能 100% 确定 64 位代码。
所以我的问题是：有没有办法检查无符号 64 位值与常量相乘结果的可靠性？我通过简单地尝试使用 UInt32 (0..$FFFFFFFF) 的所有值来检查 32 位值。这花了大约。 3分钟。检查所有 UInt64 所花的时间比我的一生要长得多。有没有办法检查使用的参数（常数，post-shift）是否可靠？
我注意到，如果选择的参数错误（但接近），DivMod100() 总是会因 00004B 这样的值而失败。是否有特殊值或范围来检查 64 位，所以我不必检查 all 值？
我当前的代码：
const
{$IF DEFINED(WIN32)}
// Checked
Div100Const = UInt32(UInt64(FFFFFFFFF) div 100 + 1);
Div100PostShift = 5;
{$ELSEIF DEFINED(WIN64)}
// Unchecked!!
Div100Const = $A3D70A3D70A3D71;
// UInt64(UInt128( FFFF FFFF FFFF FFFF) div 100 + 1);
// UInt128 is fictive type.
Div100PostShift = 2;
{$IFEND}
// Calculates X div 100 using multiplication by a constant, taking the
// high part of the 64 bit (or 128 bit) result and shifting
// right. The remainder is calculated as X - quotient * 100;
// This was tested to work safely and quickly for all values of UInt32.
function DivMod100(var X: NativeUInt): NativeUInt;
{$IFDEF WIN32}
asm
// EAX = address of X, X is UInt32 here.
PUSH EBX
MOV EDX,Div100Const
MOV ECX,EAX
MOV EAX,[ECX]
MOV EBX,EAX
MUL EDX
SHR EDX,Div100PostShift
MOV [ECX],EDX // Quotient
// Slightly faster than MUL
LEA EDX,[EDX + 4*EDX] // EDX := EDX * 5;
LEA EDX,[EDX + 4*EDX] // EDX := EDX * 5;
SHL EDX,2 // EDX := EDX * 4; 5*5*4 = 100.
MOV EAX,EBX
SUB EAX,EDX // Remainder
POP EBX
end;
{$ELSE WIN64}
asm
.NOFRAME
// RCX is address of X, X is UInt64 here.
MOV RAX,[RCX]
MOV R8,RAX
XOR RDX,RDX
MOV R9,Div100Const
MUL R9
SHR RDX,Div100PostShift
MOV [RCX],RDX // Quotient
// Faster than LEA and SHL
MOV RAX,RDX
MOV R9D,100
MUL R9
SUB R8,RAX
MOV RAX,R8 // Remainder
end;
{$ENDIF WIN32}
我找到了 libdivide.h 的解决方案。这是 Win64 稍微复杂的部分：
{$ELSE WIN64}
asm
.NOFRAME
MOV RAX,[RCX]
MOV R8,RAX
XOR RDX,RDX
MOV R9,Div100Const // New: AE147AE147AE15
MUL R9 // Preliminary result Q in RDX
// Additional part: add/shift
ADD RDX,R8 // Q := Q + X shr 1;
RCR RDX,1
SHR RDX,Div100PostShift // Q := Q shr 6;
MOV [RCX],RDX // X := Q;
// Faster than LEA and SHL
MOV RAX,RDX
MOV R9D,100
MUL R9
SUB R8,RAX
MOV RAX,R8 // Remainder
end;
{$ENDIF WIN32}
@Rudy 的回答中的代码来自以下步骤：
- 将1/100写成二进制形式：
0.000000(10100011110101110000);
- 计数小数点后的前导零：
S = 6；
- 72 个第一个有效位是：
1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101
- 四舍五入为65位；这种舍入的执行方式有某种魔力；通过对 Rudy 的答案中的常量进行逆向工程，正确的四舍五入是：
1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 1
- 删除前导
1 位：
0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0101
- 以十六进制形式写入（取回报复常数）：
A = 47 AE 14 7A E1 47 AE 15
X div 100 = (((uint128(X) * uint128(A)) shr 64) + X) shr 7 (7 = 1 + S)
像往常一样编写优化代码时，使用编译器输出作为提示/起点。可以安全地假设它所做的任何优化在一般情况下都是安全的。错误代码编译器错误很少见。
gcc 使用常数 0x28f5c28f5c28f5c3 实现无符号 64 位 divmod。我没有详细研究为除法生成常量，但是有生成它们的算法可以给出已知的良好结果（因此不需要详尽的测试）。
该代码实际上有一些重要的区别：它使用的常量与 OP 的常量不同。
请参阅评论以分析这实际上在做什么：首先除以 4，因此它可以使用一个常数，该常数仅在被除数足够小时才可用于除以 25。这也避免了以后需要添加。
#include <stdint.h>
// rem, quot ordering takes one extra instruction
struct divmod { uint64_t quotient, remainder; }
div_by_100(uint64_t x) {
struct divmod retval = { x%100, x/100 };
return retval;
}
compiles to (gcc 5.3 -O3 -mtune=haswell):
movabs rdx, 2951479051793528259
mov rax, rdi ; Function arg starts in RDI (SysV ABI)
shr rax, 2
mul rdx
shr rdx, 2
lea rax, [rdx+rdx*4] ; multiply by 5
lea rax, [rax+rax*4] ; multiply by another 5
sal rax, 2 ; imul rax, rdx, 100 is better here (Intel SnB).
sub rdi, rax
mov rax, rdi
ret
; return values in rdx:rax
使用“二进制”选项查看十六进制常量，因为反汇编程序输出是那样做的，这与 gcc 的 asm 源输出不同。
乘以 100 的部分。
gcc 使用上述 lea/lea/shl 序列，与您的问题相同。您的答案是使用 mov imm/mul 序列。
您的评论都说他们选择的版本更快。如果是这样，那是因为一些微妙的指令对齐或其他次要影响：在英特尔 SnB 系列上，它是 the same number of uops (3)，并且相同的关键路径延迟（mov imm 不在关键路径上，并且 mul 是 3 个周期）。
clang uses 我认为最好的选择 (imul rax, rdx, 100)。之前看到clang选的就想到了，不重要。那是 1 个融合域 uop（只能在 p0 上执行），仍然有 3c 延迟。因此，如果您使用此例程进行多精度处理时受延迟限制，它可能无济于事，但它是最佳选择。 （如果您受延迟限制，将您的代码内联到一个循环中而不是通过内存传递其中一个参数可以节省很多周期。）
imul 有效，因为 。 mul 没有 2 或 3 操作数形式，因为无论输入的有符号或无符号解释如何，结果的低半部分都是相同的。
顺便说一句，带有 -march=native 的 clang 使用 mulx 作为 64x64->128，而不是 mul，但它没有获得任何好处。根据 Agner Fog 的表格，它比 mul.
延迟一个周期。
AMD 的 imul r,r,i（尤其是 64b 版本）的延迟比 3c 还差，这也许就是 gcc 避免它的原因。 IDK 多少 gcc 维护者投入了多少工作来调整成本，所以像 -mtune=haswell 这样的设置工作得很好，但是 lot 的代码没有用任何 -mtune 设置编译（即使-march) 暗示了一个），所以当 gcc 做出最适合旧 CPU 或 AMD 的选择时，我并不感到惊讶。
clang 仍然使用 imul r64, r64, imm 和 -mtune=bdver1 (Bulldozer)，这节省了 m-ops，但比使用 lea/lea/shl 的延迟多了 1c。 （比例>1 的 lea 在 Bulldozer 上是 2c 延迟）。
对于我的 BigInteger 代码，对于非常大的 BigInteger，输出结果很慢。所以现在我使用递归的分而治之算法，仍然需要2'30"将当前已知最大的素数转换为超过2200万位的十进制字符串（但只需135毫秒即可将其转换为十六进制字符串） .
我仍然想减少时间，所以我需要一个可以非常快速地将 NativeUInt（即 32 位平台上的 UInt32，64 位平台上的 UInt64）除以 100 的例程。所以我用常量乘法。这在 32 位代码中工作正常，但我不能 100% 确定 64 位代码。
所以我的问题是：有没有办法检查无符号 64 位值与常量相乘结果的可靠性？我通过简单地尝试使用 UInt32 (0..$FFFFFFFF) 的所有值来检查 32 位值。这花了大约。 3分钟。检查所有 UInt64 所花的时间比我的一生要长得多。有没有办法检查使用的参数（常数，post-shift）是否可靠？
我注意到，如果选择的参数错误（但接近），DivMod100() 总是会因 00004B 这样的值而失败。是否有特殊值或范围来检查 64 位，所以我不必检查 all 值？
我当前的代码：
const
{$IF DEFINED(WIN32)}
// Checked
Div100Const = UInt32(UInt64(FFFFFFFFF) div 100 + 1);
Div100PostShift = 5;
{$ELSEIF DEFINED(WIN64)}
// Unchecked!!
Div100Const = $A3D70A3D70A3D71;
// UInt64(UInt128( FFFF FFFF FFFF FFFF) div 100 + 1);
// UInt128 is fictive type.
Div100PostShift = 2;
{$IFEND}
// Calculates X div 100 using multiplication by a constant, taking the
// high part of the 64 bit (or 128 bit) result and shifting
// right. The remainder is calculated as X - quotient * 100;
// This was tested to work safely and quickly for all values of UInt32.
function DivMod100(var X: NativeUInt): NativeUInt;
{$IFDEF WIN32}
asm
// EAX = address of X, X is UInt32 here.
PUSH EBX
MOV EDX,Div100Const
MOV ECX,EAX
MOV EAX,[ECX]
MOV EBX,EAX
MUL EDX
SHR EDX,Div100PostShift
MOV [ECX],EDX // Quotient
// Slightly faster than MUL
LEA EDX,[EDX + 4*EDX] // EDX := EDX * 5;
LEA EDX,[EDX + 4*EDX] // EDX := EDX * 5;
SHL EDX,2 // EDX := EDX * 4; 5*5*4 = 100.
MOV EAX,EBX
SUB EAX,EDX // Remainder
POP EBX
end;
{$ELSE WIN64}
asm
.NOFRAME
// RCX is address of X, X is UInt64 here.
MOV RAX,[RCX]
MOV R8,RAX
XOR RDX,RDX
MOV R9,Div100Const
MUL R9
SHR RDX,Div100PostShift
MOV [RCX],RDX // Quotient
// Faster than LEA and SHL
MOV RAX,RDX
MOV R9D,100
MUL R9
SUB R8,RAX
MOV RAX,R8 // Remainder
end;
{$ENDIF WIN32}
我找到了 libdivide.h 的解决方案。这是 Win64 稍微复杂的部分：
{$ELSE WIN64}
asm
.NOFRAME
MOV RAX,[RCX]
MOV R8,RAX
XOR RDX,RDX
MOV R9,Div100Const // New: AE147AE147AE15
MUL R9 // Preliminary result Q in RDX
// Additional part: add/shift
ADD RDX,R8 // Q := Q + X shr 1;
RCR RDX,1
SHR RDX,Div100PostShift // Q := Q shr 6;
MOV [RCX],RDX // X := Q;
// Faster than LEA and SHL
MOV RAX,RDX
MOV R9D,100
MUL R9
SUB R8,RAX
MOV RAX,R8 // Remainder
end;
{$ENDIF WIN32}
@Rudy 的回答中的代码来自以下步骤：
- 将1/100写成二进制形式：
0.000000(10100011110101110000); - 计数小数点后的前导零：
S = 6； - 72 个第一个有效位是：
1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101
- 四舍五入为65位；这种舍入的执行方式有某种魔力；通过对 Rudy 的答案中的常量进行逆向工程，正确的四舍五入是：
1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 1
- 删除前导
1位：
0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0101
- 以十六进制形式写入（取回报复常数）：
A = 47 AE 14 7A E1 47 AE 15
X div 100 = (((uint128(X) * uint128(A)) shr 64) + X) shr 7(7 = 1 + S)
像往常一样编写优化代码时，使用编译器输出作为提示/起点。可以安全地假设它所做的任何优化在一般情况下都是安全的。错误代码编译器错误很少见。
gcc 使用常数 0x28f5c28f5c28f5c3 实现无符号 64 位 divmod。我没有详细研究为除法生成常量，但是有生成它们的算法可以给出已知的良好结果（因此不需要详尽的测试）。
该代码实际上有一些重要的区别：它使用的常量与 OP 的常量不同。
请参阅评论以分析这实际上在做什么：首先除以 4，因此它可以使用一个常数，该常数仅在被除数足够小时才可用于除以 25。这也避免了以后需要添加。
#include <stdint.h>
// rem, quot ordering takes one extra instruction
struct divmod { uint64_t quotient, remainder; }
div_by_100(uint64_t x) {
struct divmod retval = { x%100, x/100 };
return retval;
}
compiles to (gcc 5.3 -O3 -mtune=haswell):
movabs rdx, 2951479051793528259
mov rax, rdi ; Function arg starts in RDI (SysV ABI)
shr rax, 2
mul rdx
shr rdx, 2
lea rax, [rdx+rdx*4] ; multiply by 5
lea rax, [rax+rax*4] ; multiply by another 5
sal rax, 2 ; imul rax, rdx, 100 is better here (Intel SnB).
sub rdi, rax
mov rax, rdi
ret
; return values in rdx:rax
使用“二进制”选项查看十六进制常量，因为反汇编程序输出是那样做的，这与 gcc 的 asm 源输出不同。
乘以 100 的部分。
gcc 使用上述 lea/lea/shl 序列，与您的问题相同。您的答案是使用 mov imm/mul 序列。
您的评论都说他们选择的版本更快。如果是这样，那是因为一些微妙的指令对齐或其他次要影响：在英特尔 SnB 系列上，它是 the same number of uops (3)，并且相同的关键路径延迟（mov imm 不在关键路径上，并且 mul 是 3 个周期）。
clang uses 我认为最好的选择 (imul rax, rdx, 100)。之前看到clang选的就想到了，不重要。那是 1 个融合域 uop（只能在 p0 上执行），仍然有 3c 延迟。因此，如果您使用此例程进行多精度处理时受延迟限制，它可能无济于事，但它是最佳选择。 （如果您受延迟限制，将您的代码内联到一个循环中而不是通过内存传递其中一个参数可以节省很多周期。）
imul 有效，因为 mul 没有 2 或 3 操作数形式，因为无论输入的有符号或无符号解释如何，结果的低半部分都是相同的。
顺便说一句，带有 -march=native 的 clang 使用 mulx 作为 64x64->128，而不是 mul，但它没有获得任何好处。根据 Agner Fog 的表格，它比 mul.
AMD 的 imul r,r,i（尤其是 64b 版本）的延迟比 3c 还差，这也许就是 gcc 避免它的原因。 IDK 多少 gcc 维护者投入了多少工作来调整成本，所以像 -mtune=haswell 这样的设置工作得很好，但是 lot 的代码没有用任何 -mtune 设置编译（即使-march) 暗示了一个），所以当 gcc 做出最适合旧 CPU 或 AMD 的选择时，我并不感到惊讶。
clang 仍然使用 imul r64, r64, imm 和 -mtune=bdver1 (Bulldozer)，这节省了 m-ops，但比使用 lea/lea/shl 的延迟多了 1c。 （比例>1 的 lea 在 Bulldozer 上是 2c 延迟）。