为啥GCC在实现整数除法时使用乘以一个奇怪的数字？

Posted 2023-02-14

【中文标题】为啥GCC在实现整数除法时使用乘以一个奇怪的数字？

【英文标题】：Why does GCC use multiplication by a strange number in implementing integer division?为什么GCC在实现整数除法时使用乘以一个奇怪的数字？

【发布时间】：2017-05-02 05:21:16

【问题描述】：

我一直在阅读有关 div 和 mul 汇编操作的文章，我决定通过用 C 编写一个简单的程序来看看它们的实际效果：

文件分割.c

#include <stdlib.h>
#include <stdio.h>

int main()

    size_t i = 9;
    size_t j = i / 5;
    printf("%zu\n",j);
    return 0;

然后生成汇编语言代码：

gcc -S division.c -O0 -masm=intel

但是查看生成的division.s 文件，它不包含任何 div 操作！相反，它使用位移和幻数执行某种黑魔法。这是一个计算i/5的代码sn-p：

mov     rax, QWORD PTR [rbp-16]   ; Move i (=9) to RAX
movabs  rdx, -3689348814741910323 ; Move some magic number to RDX (?)
mul     rdx                       ; Multiply 9 by magic number
mov     rax, rdx                  ; Take only the upper 64 bits of the result
shr     rax, 2                    ; Shift these bits 2 places to the right (?)
mov     QWORD PTR [rbp-8], rax    ; Magically, RAX contains 9/5=1 now, 
                                  ; so we can assign it to j

这里发生了什么？为什么 GCC 根本不使用 div？它是如何生成这个幻数的，为什么一切都有效？

【问题讨论】：

gcc 优化常量除法，尝试除法 2,3,4,5,6,7,8，您很可能会看到每种情况下非常不同的代码。

注意：幻数 -3689348814741910323 转换为 CCCCCCCCCCCCCCCD 作为 uint64_t 或大约 (2^64)*4/5。

@qiubit ：编译器不会因为优化被禁用而反常地生成低效的代码。例如，将执行不涉及代码重新排序或变量消除的琐碎“优化”。本质上，单个源语句将单独转换为该操作的最有效代码。编译器优化会考虑周围的代码，而不仅仅是单个语句。

阅读这篇精彩的文章：Labor of Division

一些编译器实际上会因为优化被禁用而反常地生成低效的代码。特别是，他们会这样做以简化调试，例如在各个代码行上设置断点的能力。事实上，GCC 非常不寻常，因为它没有真正的“无优化”模式，因为它的许多优化都是本构打开的。这是一个使用 GCC 可以看到的示例。另一方面，Clang 和 MSVC将在 -O0 处发出 div 指令。 （抄送@克利福德）

【参考方案1】：

除以 5 与乘以 1/5 相同，这又与乘以 4/5 并右移 2 位相同。相关的值是十六进制的CCCCCCCCCCCCCCCD，如果放在十六进制点之后，它是 4/5 的二进制表示（即五分之四的二进制是 0.110011001100 重复出现 - 原因见下文）。我想你可以从这里拿走它！您可能想查看fixed point arithmetic（但请注意它在末尾四舍五入为整数）。

至于为什么，乘法比除法快，而当除数固定时，这是一条更快的路线。

请参阅Reciprocal Multiplication, a tutorial 了解有关其工作原理的详细文章，并根据定点进行解释。它展示了求倒数的算法是如何工作的，以及如何处理有符号的除法和模数。

让我们考虑一下为什么0.CCCCCCCC...（十六进制）或0.110011001100...二进制是4/5。将二进制表示除以 4（右移 2 位），我们将得到 0.001100110011...，通过简单的检查可以将原始添加得到 0.111111111111...，这显然等于 1，与 0.9999999... 中的方法相同十进制等于一。因此，我们知道x + x/4 = 1，所以5x/4 = 1，x=4/5。然后将其表示为十六进制的 CCCCCCCCCCCCD 以进行舍入（因为存在的最后一位之外的二进制数字将是 1）。

【讨论】：

@user2357112 随时发布您自己的答案，但我不同意。您可以将乘法视为 64.0 位乘以 0.64 位乘以给出 128 位定点答案，其中最低 64 位被丢弃，然后除以 4（正如我在第一段中指出的那样）。您可能会想出一个替代的模算术答案，它同样可以很好地解释位移动，但我很确定这可以作为解释。

这个值实际上是“CCCCCCCCCCCCCCCD” 最后一个D很重要，它确保当结果被截断时，精确的除法得到正确的答案。

没关系。我没有看到他们取的是 128 位乘法结果的高 64 位；这不是你可以用大多数语言做的事情，所以我最初并没有意识到它正在发生。通过明确提及获取 128 位结果的高 64 位如何等效于乘以定点数并向下舍入，这个答案将得到很大改善。 （另外，最好解释一下为什么它必须是 4/5 而不是 1/5，以及为什么我们必须向上舍入而不是向下舍入 4/5。）

因为您必须计算出在舍入边界上向上除以 5 需要多大的错误，然后将其与计算中的最坏情况错误进行比较。大概 gcc 开发人员已经这样做了，并得出结论，它总是会给出正确的结果。

实际上您可能只需要检查 5 个可能的最高输入值，如果这些值正确，其他所有值也应该如此。

【参考方案2】：

我会从稍微不同的角度回答：因为是允许的。

C 和 C++ 是针对抽象机器定义的。编译器按照as-if 规则将该程序从抽象机器转换为具体机器。

只要不改变抽象机指定的可观察行为，编译器就可以进行任何更改。没有合理的期望编译器会以最直接的方式转换你的代码（即使很多 C 程序员都这么认为）。通常，它这样做是因为编译器希望与直接方法相比优化性能（如其他答案中详细讨论的那样）。 如果在任何情况下编译器将正确的程序“优化”为具有不同可观察行为的东西，这就是编译器错误。 我们的代码中的任何未定义行为（有符号整数溢出是一个经典示例）并且此合同无效。

【参考方案3】：

这里是一个算法文档的链接，该算法生成我在 Visual Studio 中看到的值和代码（在大多数情况下），我假设它仍然在 GCC 中用于将可变整数除以常量整数。

http://gmplib.org/~tege/divcnst-pldi94.pdf

文章中，一个uword有N位，一个udword有2N位，n=分子=被除数，d=分母=除数，ℓ初始设置为ceil(log2(d))，shpre为预移位(在乘法之前使用）= e = d 中尾随零位的数量，shpost 是移位后（在乘法之后使用），prec 是精度 = N - e = N - shpre。目标是使用预移位、乘法和后移位来优化 n/d 的计算。

向下滚动到图 6.2，它定义了如何生成 udword 乘数（最大大小为 N+1 位），但没有清楚地解释该过程。我将在下面解释这一点。

图 4.2 和图 6.2 显示了对于大多数除数，乘法器如何减少到 N 位或更少的乘法器。公式 4.5 解释了用于处理图 4.1 和 4.2 中的 N+1 位乘法器的公式是如何推导出来的。

在现代 X86 和其他处理器的情况下，乘法时间是固定的，因此预移位对这些处理器没有帮助，但它仍然有助于将乘法器从 N+1 位减少到 N 位。我不知道 GCC 或 Visual Studio 是否已经消除了 X86 目标的预移位。

回到图 6.2。只有当分母（除数）> 2^(N-1)（当 ℓ == N => mlow = 2^(2N)）时，mlow 和 mhigh 的分子（除数）才能大于 udword，在这种情况下n/d 的优化替换是比较（如果 n>=d，q = 1，否则 q = 0），因此不会生成乘数。 mlow 和 mhigh 的初始值将是 N+1 位，并且可以使用两个 udword/uword 除法来产生每个 N+1 位值（mlow 或 mhigh）。以64位模式下的X86为例：

; upper 8 bytes of dividend = 2^(ℓ) = (upper part of 2^(N+ℓ))
; lower 8 bytes of dividend for mlow  = 0
; lower 8 bytes of dividend for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e)
dividend  dq    2 dup(?)        ;16 byte dividend
divisor   dq    1 dup(?)        ; 8 byte divisor

; ...
        mov     rcx,divisor
        mov     rdx,0
        mov     rax,dividend+8     ;upper 8 bytes of dividend
        div     rcx                ;after div, rax == 1
        mov     rax,dividend       ;lower 8 bytes of dividend
        div     rcx
        mov     rdx,1              ;rdx:rax = N+1 bit value = 65 bit value

您可以使用 GCC 进行测试。您已经了解了 j = i/5 的处理方式。看看 j = i/7 是如何处理的（应该是 N+1 位乘法器的情况）。

在大多数当前处理器上，乘法具有固定的时序，因此不需要预移位。对于 X86，最终结果是大多数除数的两个指令序列，以及像 7 这样的除数的五个指令序列（为了模拟 N+1 位乘法器，如等式 4.5 和 pdf 文件的图 4.2 所示）。示例 X86-64 代码：

;       rax = dividend, rbx = 64 bit (or less) multiplier, rcx = post shift count
;       two instruction sequence for most divisors:

        mul     rbx                     ;rdx = upper 64 bits of product
        shr     rdx,cl                  ;rdx = quotient
;
;       five instruction sequence for divisors like 7
;       to emulate 65 bit multiplier (rbx = lower 64 bits of multiplier)

        mul     rbx                     ;rdx = upper 64 bits of product
        sub     rbx,rdx                 ;rbx -= rdx
        shr     rbx,1                   ;rbx >>= 1
        add     rdx,rbx                 ;rdx = upper 64 bits of corrected product
        shr     rdx,cl                  ;rdx = quotient
;       ...

【讨论】：

那篇论文描述了在 gcc 中实现它，所以我认为仍然使用相同的算法是一个安全的假设。

那篇 1994 年的论文描述了在 gcc 中实现它，因此 gcc 有时间更新其算法。以防其他人没有时间查看该 URL 中的 94 的含义。

【参考方案4】：

整数除法是您可以在现代处理器上执行的最慢的算术运算之一，延迟高达几十个周期并且吞吐量很差。 （对于 x86，请参阅 Agner Fog's instruction tables and microarch guide）。

如果您提前知道除数，则可以通过将其替换为具有等效效果的一组其他运算（乘法、加法和移位）来避免除法。即使需要多次运算，它通常仍然比整数除法本身快很多。

以这种方式实现 C / 运算符，而不是使用涉及 div 的多指令序列只是 GCC 进行常量除法的默认方式。它不需要跨操作进行优化，即使是调试也不会改变任何东西。 （不过，对小代码使用 -Os 确实让 GCC 使用 div。）使用乘法逆而不是除法就像使用 lea 而不是 mul 和 add

因此，如果除数在编译时未知，您只会在输出中看到 div 或 idiv。

有关编译器如何生成这些序列的信息，以及让您自己生成它们的代码（几乎肯定没有必要，除非您使用的是 Braindead 编译器），请参阅libdivide。

【讨论】：

我不确定在速度比较中将 FP 和整数运算混为一谈是否公平，@fuz。也许 Sneftel 应该说 division 是您可以在现代处理器上执行的最慢的 integer 操作？此外，在 cmets 中提供了一些链接，以进一步解释这种“魔法”。您认为它们适合收集在您的答案中以提高知名度吗？ 1, 2, 3

因为操作顺序在功能上是相同的......这始终是一个要求，即使在-O3 也是如此。编译器必须编写为所有可能的输入值提供正确结果的代码。这仅适用于 -ffast-math 的浮点数，并且 AFAIK 没有“危险”整数优化。 （启用优化后，编译器可能能够证明可能的值范围，例如，它可以使用仅适用于非负有符号整数的东西。）

真正的答案是 gcc -O0 still transforms code through internal representations as part of turning C into machine code。碰巧即使在-O0（但不是-Os）默认情况下也会启用模乘逆。其他编译器（如 clang）将 DIV 用于-O0 处的非 2 次幂常量。相关：我想我在my Collatz-conjecture hand-written asm answer

中包含了一段关于此的内容

@PeterCordes 是的，我认为 GCC（以及许多其他编译器）已经忘记为“禁用优化时应用什么样的优化”提供一个很好的理由。花了一天的大部分时间来追踪一个不起眼的代码生成错误，我现在对此有点恼火。

@Sneftel：这可能只是因为主动向编译器开发者抱怨他们的代码运行速度比预期快的应用程序开发者的数量相对较少。

【参考方案5】：

一般来说，乘法比除法快得多。因此，如果我们可以避免乘以倒数，而是可以显着加快除以常数

一个问题是我们不能精确地表示倒数（除非除法是 2 的幂，但在这种情况下，我们通常可以将除法转换为位移）。因此，为了确保正确答案，我们必须小心，倒数中的错误不会导致最终结果中的错误。

-3689348814741910323 是 0xCCCCCCCCCCCCCCCD 是 4/5 以上的值，用 0.64 定点表示。

当我们将 64 位整数乘以 0.64 定点数时，我们得到 64.64 的结果。我们将值截断为 64 位整数（实际上将其舍入到零），然后执行进一步的移位，除以 4 并再次截断通过查看位级别，很明显我们可以将两个截断视为单个截断。

这显然给了我们至少一个除以 5 的近似值，但它给我们一个准确的答案，正确地接近零吗？

要获得准确的答案，误差必须足够小，以免将答案推到舍入边界。

除以 5 的准确答案总是小数部分为 0、1/5、2/5、3/5 或 4/5。因此，相乘和移位结果中小于 1/5 的正误差永远不会使结果超出舍入边界。

我们常量的误差是 (1/5) * 2^-64。 i 的值小于 2⁶⁴，所以相乘后的误差小于 1/5。除以 4 后，误差小于 (1/5) * 2^-2。

(1/5) * 2⁻²

---

不幸的是，这不适用于所有除数。

如果我们尝试将 4/7 表示为 0.64 定点数并从零四舍五入，我们最终会得到 (6/7) * 2^-64 的错误。乘以一个略低于 2⁶⁴ 的 i 值后，我们最终得到一个略低于 6/7 的误差，除以四之后，我们最终得到一个略低于 1.5/7 的误差，大于1/7。

所以要正确实现除以 7，我们需要乘以 0.65 的定点数。我们可以通过乘以定点数的低 64 位来实现，然后加上原始数（这可能会溢出进位），然后通过进位进行循环。

【讨论】：

这个答案将模乘逆从“看起来比我想花时间更复杂的数学”变成有意义的东西。 +1 易于理解的版本。除了使用编译器生成的常量之外，我从不需要做任何事情，所以我只浏览了其他解释数学的文章。

我在代码中根本看不出与模运算有任何关系。不知道其他评论者是从哪里得到的。

它是模 2^n，就像寄存器中的所有整数数学一样。 en.wikipedia.org/wiki/…

@PeterCordes 模乘逆用于精确除法，afaik 它们对一般除法没有用

@PeterCordes 乘以定点倒数？我不知道大家怎么称呼它，但我可能会这样称呼它，它相当具有描述性