6

对于我的 BigInteger 代码，对于非常大的 BigInteger，输出结果很慢。所以现在我使用递归分治算法，它仍然需要 2'30" 才能将当前最大的已知素数转换为超过 2200 万位的十进制字符串（但只需 135 毫秒即可将其转换为十六进制字符串） .

我仍然想减少时间，所以我需要一个可以非常快地将 NativeUInt（即 32 位平台上的 UInt32，64 位平台上的 UInt64）除以 100 的例程。所以我用常数乘法。这在 32 位代码中可以正常工作，但我不能 100% 确定 64 位。

所以我的问题是：有没有办法检查无符号 64 位值乘以常数的结果的可靠性？我通过简单地尝试 UInt32 (0..$FFFFFFFF) 的所有值来检查 32 位值。这花了大约。3分钟。检查所有 UInt64 所花费的时间比我有生之年要长得多。有没有办法检查所使用的参数（常量、后移）是否可靠？

我注意到，如果选择的参数错误（但接近） ，DivMod100()总是会失败。$4000004B是否有特殊值或范围来检查 64 位，所以我不必检查所有值？

我当前的代码：

const
{$IF DEFINED(WIN32)}
  // Checked
  Div100Const = UInt32(UInt64($1FFFFFFFFF) div 100 + 1);
  Div100PostShift = 5;
{$ELSEIF DEFINED(WIN64)}
  // Unchecked!!
  Div100Const = $A3D70A3D70A3D71; 
  // UInt64(UInt128($3 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}
4

3 回答 3

2

像往常一样编写优化代码时，使用编译器输出作为提示/起点。可以安全地假设它所做的任何优化在一般情况下都是安全的。错误代码的编译器错误很少见。

gcc 使用常量0x28f5c28f5c28f5c3. 我没有详细研究生成除法常量，但是有一些算法可以生成它们，它们会给出已知的良好结果（因此不需要详尽的测试）。

代码实际上有一些重要的区别：它使用的常量与 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;
}

编译为（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 系列上，它是相同数量的 uops (3)和相同的关键路径延迟（mov imm不在关键路径上，并且mul是 3 个周期） .

clang 使用我认为最好的选择 ( imul rax, rdx, 100)。我在看到clang选择它之前就想到了它，这并不重要。那是 1 个融合域 uop（只能在 p0 上执行），仍然有 3c 延迟。因此，如果您在使用此例程进行多精度时受到延迟限制，它可能无济于事，但它是最佳选择。（如果您受延迟限制，将代码内联到循环中而不是通过内存传递参数之一可以节省大量循环。）

imul有效，因为您只使用 result 的低 64b。没有 2 或 3 操作数形式，mul因为无论输入的有符号或无符号解释如何，结果的低半部分都是相同的。

顺便说一句，-march=native使用mulx64x64->128，而不是mul，但不会从中获得任何好处。根据 Agner Fog 的表格，延迟比mul.

AMD 的延迟低于 3c imul r,r,i（尤其是 64b 版本），这可能是 gcc 避免它的原因。IDK gcc 维护人员在调整成本方面投入了多少工作，以使设置运行-mtune=haswell良好，但是很多代码没有使用任何-mtune设置编译（即使是 隐含的设置-march），所以当 gcc 做出最适合旧版本的选择时，我并不感到惊讶CPU，或用于 AMD。

clang 仍然使用imul r64, r64, immwith -mtune=bdver1(Bulldozer)，它可以节省 m-ops，但比使用 lea/lea/shl 的延迟要多 1c。（在 Bulldozer 上 scale>1 的 lea 是 2c 延迟）。

于 2016-01-31T23:25:43.367 回答
1

我找到了libdivide.h的解决方案。这是 Win64 稍微复杂的部分：

{$ELSE WIN64}
asm
        .NOFRAME

        MOV     RAX,[RCX]
        MOV     R8,RAX
        XOR     RDX,RDX
        MOV     R9,Div100Const       // New: $47AE147AE147AE15
        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}
于 2016-01-30T15:42:34.877 回答
1

@Rudy 的答案中的代码来自以下步骤：

  1. 以二进制形式写入 1/100：0.000000(10100011110101110000);
  2. 计算小数点后的前导零：S = 6;
  3. 72 个第一有效位是：

1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101

  1. 舍入到 65 位；这种舍入的方式有某种魔力；通过从鲁迪的答案中对常数进行逆向工程，正确的舍入是：

1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 0011 1101 0111 0000 1010 1

  1. 删除前导1位：

0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0100 0111 1010 1110 0001 0101

  1. 以十六进制形式写（取回被报复的常数）：

A = 47 AE 14 7A E1 47 AE 15

  1. X div 100 = (((uint128(X) * uint128(A)) shr 64) + X) shr 7 (7 = 1 + S)
于 2016-01-30T18:15:01.510 回答