c++ - 是否有一种无分支方法可以快速找到两个 double 浮点值的最小值/最大值？

Question

我有两个 double ，a和b，它们都在[0,1]中。由于性能原因，我希望a和b的最小值/最大值而不进行分支。

假设a和b均为正且小于1，是否有一种有效的方法来获取两者的最小值/最大值？理想情况下，我不希望分支。

Answer 1

是的，有一种方法可以计算两个double的最大值或最小值，而无需任何分支。这样做的C++代码如下所示:

#include <algorithm>

double FindMinimum(double a, double b)
{
    return std::min(a, b);
}

double FindMaximum(double a, double b)
{
    return std::max(a, b);
}

我敢打赌，您以前见过。唯恐您不相信这是无分支的check out the disassembly:

FindMinimum(double, double):
    minsd   xmm1, xmm0
    movapd  xmm0, xmm1
    ret

FindMaximum(double, double):
    maxsd   xmm1, xmm0
    movapd  xmm0, xmm1
    ret

这就是从所有针对x86的流行编译器中获得的。使用SSE2指令集，特别是minsd/maxsd指令，该指令无分支地评估两个 double 浮点值的最小值/最大值。
所有64位x86处理器都支持SSE2； AMD64扩展需要它。即使是大多数不带64位的x86处理器也支持SSE2。它于2000年发布。您必须走很长一段路才能找到不支持SSE2的处理器。但是，如果您呢？好吧，即使在那里you get branchless code on most popular compilers:

FindMinimum(double, double):
    fld      QWORD PTR [esp + 12]
    fld      QWORD PTR [esp + 4]
    fucomi   st(1)
    fcmovnbe st(0), st(1)
    fstp     st(1)
    ret

FindMaximum(double, double):
    fld      QWORD PTR [esp + 4]
    fld      QWORD PTR [esp + 12]
    fucomi   st(1)
    fxch     st(1)
    fcmovnbe st(0), st(1)
    fstp     st(1)
    ret

fucomi指令执行比较，设置标志，然后fcmovnbe指令根据这些标志的值执行条件移动。这一切都是完全无分支的，并依赖于1995年Pentium Pro引入x86 ISA的指令，该指令自Pentium II以来在所有x86芯片上均受支持。
此处唯一不会生成无分支代码的编译器是MSVC，因为it doesn't take advantage of the FCMOVxx instruction。相反，您得到:

double FindMinimum(double, double) PROC
    fld     QWORD PTR [a]
    fld     QWORD PTR [b]
    fcom    st(1)            ; compare "b" to "a"
    fnstsw  ax               ; transfer FPU status word to AX register
    test    ah, 5            ; check C0 and C2 flags
    jp      Alt
    fstp    st(1)            ; return "b"
    ret
Alt:
    fstp    st(0)            ; return "a"
    ret
double FindMinimum(double, double) ENDP

double FindMaximum(double, double) PROC
    fld     QWORD PTR [b]
    fld     QWORD PTR [a]
    fcom    st(1)            ; compare "b" to "a"
    fnstsw  ax               ; transfer FPU status word to AX register
    test    ah, 5            ; check C0 and C2 flags
    jp      Alt
    fstp    st(0)            ; return "b"
    ret
Alt:
    fstp    st(1)            ; return "a"
    ret
double FindMaximum(double, double) ENDP

注意分支JP指令(如果设置了奇偶校验位则跳转)。 FCOM指令用于进行比较，这是基本x87 FPU指令集的一部分。不幸的是，这会在FPU状态字中设置标志，因此为了分支这些标志，需要将其提取。这就是FNSTSW指令的目的，该指令将x87 FPU状态字存储到通用的AX寄存器中(它也可以存储到内存中，但是……为什么？)。然后，该代码TEST为适当的位，并进行相应分支以确保返回正确的值。除了分支之外，检索FPU状态字也将相对较慢。这就是Pentium Pro引入FCOM指令的原因。
但是，通过位旋转操作确定最小/最大，您不太可能提高任何代码的速度。有两个基本原因:

唯一生成低效率代码的编译器是MSVC，没有什么好的方法来强制它生成所需的指令。尽管MSVC支持内联汇编用于32位x86目标it is a fool's errand when seeking performance improvements。我还将引用自己:

Inline assembly disrupts the optimizer in rather significant ways, so unless you're writing significant swaths of code in inline assembly, there is unlikely to be a substantial net performance gain. Furthermore, Microsoft's inline assembly syntax is extremely limited. It trades flexibility for simplicity in a big way. In particular, there is no way to specify input values, so you're stuck loading the input from memory into a register, and the caller is forced to spill the input from a register to memory in preparation. This creates a phenomenon I like to call "a whole lotta shufflin' goin' on", or for short, "slow code". You don't drop to inline assembly in cases where slow code is acceptable. Thus, it is always preferable (at least on MSVC) to figure out how to write C/C++ source code that persuades the compiler to emit the object code you want. Even if you can only get close to the ideal output, that's still considerably better than the penalty you pay for using inline assembly.

为了访问浮点值的原始位，您必须进行域转换，从浮点到整数，然后再回到浮点。这很慢，尤其是在没有SSE2的情况下，因为从x87 FPU到ALU中的通用整数寄存器获取值的唯一方法是间接通过内存。

如果您仍然想采用这种策略(例如，对其进行基准测试)，则可以利用以下事实:浮点值按照其IEEE 754表示法按字典顺序排序，除了符号位。因此，由于您假设两个值都是正值:

FindMinimumOfTwoPositiveDoubles(double a, double b):
    mov   rax, QWORD PTR [a]
    mov   rdx, QWORD PTR [b]
    sub   rax, rdx              ; subtract bitwise representation of the two values
    shr   rax, 63               ; isolate the sign bit to see if the result was negative
    ret

FindMaximumOfTwoPositiveDoubles(double a, double b):
    mov   rax, QWORD PTR [b]    ; \ reverse order of parameters
    mov   rdx, QWORD PTR [a]    ; /  for the SUB operation
    sub   rax, rdx
    shr   rax, 63
    ret

或者，为避免内联汇编:

bool FindMinimumOfTwoPositiveDoubles(double a, double b)
{
    static_assert(sizeof(a) == sizeof(uint64_t),
                  "A double must be the same size as a uint64_t for this bit manipulation to work.");
    const uint64_t aBits = *(reinterpret_cast<uint64_t*>(&a));
    const uint64_t bBits = *(reinterpret_cast<uint64_t*>(&b));
    return ((aBits - bBits) >> ((sizeof(uint64_t) * CHAR_BIT) - 1));
}

bool FindMaximumOfTwoPositiveDoubles(double a, double b)
{
    static_assert(sizeof(a) == sizeof(uint64_t),
                  "A double must be the same size as a uint64_t for this bit manipulation to work.");
    const uint64_t aBits = *(reinterpret_cast<uint64_t*>(&a));
    const uint64_t bBits = *(reinterpret_cast<uint64_t*>(&b));
    return ((bBits - aBits) >> ((sizeof(uint64_t) * CHAR_BIT) - 1));
}

请注意，此实现存在一些严重警告。特别是，如果两个浮点值具有不同的符号，或者两个值都为负，则它将中断。如果两个值均为负，则可以修改代码以翻转其符号，进行比较，然后返回相反的值。要处理两个值具有不同符号的情况，可以添加代码以检查符号位。

    // ...

    // Enforce two's-complement lexicographic ordering.
    if (aBits < 0)
    {
        aBits = ((1 << ((sizeof(uint64_t) * CHAR_BIT) - 1)) - aBits);
    }
    if (bBits < 0)
    {
        bBits = ((1 << ((sizeof(uint64_t) * CHAR_BIT) - 1)) - bBits);
    }

    // ...

处理负零也将是一个问题。 IEEE 754表示+0.0等于-0.0，因此您的比较函数必须决定是否要将这些值视为不同，或者向比较例程添加特殊代码以确保将负零和正零视为等效。
添加所有这些特殊情况的代码肯定会降低性能，以至于我们无法通过简单的浮点比较来达到收支平衡，并且很可能最终会变得更慢。