为什么GCC在实现整数除法时使用乘以奇数的乘法？

Question

我一直在阅读div和mul组装操作,我决定通过在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:

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

这里发生了什么？为什么海湾合作委员会根本不使用div？它如何产生这个神奇的数字以及为什么一切都有效？

Answer 1

整数除法是您可以在现代处理器上执行的最慢的算术运算之一,延迟可达数十个周期和吞吐量不佳.(对于x86,请参阅Agner Fog的说明表和微指南).

如果您提前知道除数,则可以通过将其替换为具有相同效果的一组其他运算(乘法,加法和移位)来避免除法.即使需要进行多次操作,它通常仍然比整数除法本身快得多.

以/这种方式实现C 运算符而不是涉及多指令序列div只是GCC默认的常量除法.它不需要跨操作进行优化,即使是调试也不会改变任何内容.(虽然使用-Os小代码大小确实可以使用GCC div.)使用乘法逆而不是除法就像使用lea而不是mul和add

因此,如果在编译时不知道除数,则只倾向于看到div或idiv在输出中.

有关编译器如何生成这些序列的信息,以及允许您自己生成它们的代码(几乎肯定是不必要的,除非您使用脑死亡编译器),请参阅libdivide.

Answer 2

除以5与乘以1/5相同,再乘以4/5并向右移2位相同.有关的值是CCCCCCCCCCCCCCCD十六进制,如果放在一个十六进制点后面是4/5的二进制表示(即四分之五的二进制0.110011001100重复出现 - 见下面的原因).我想你可以从这里拿走它!您可能想要检查定点算术(尽管注意它在最后舍入为整数.

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

请参阅Reciprocal Multiplication,这是一个关于它如何工作的详细文章的教程,用定点来解释.它显示了查找倒数的算法如何工作,以及如何处理有符号的除法和模数.

让我们考虑一下为什么0.CCCCCCCC...(十六进制)或0.110011001100...二进制是4/5.将二进制表示除以4(右移2位),我们将得到0.001100110011...通过平凡检查可以添加原始得到的0.111111111111...,显然等于1,0.9999999...十进制中的相同方式等于1.因此,我们知道x + x/4 = 1,所以5x/4 = 1,x=4/5.然后将其表示为CCCCCCCCCCCCD十六进制以进行舍入(因为超出最后一个的二进制数字将是a 1).

Answer 3

通常,乘法比除法快得多.因此,如果我们可以通过乘以倒数来逃避,我们可以通过常数显着加快除法

皱纹是我们不能准确地表示倒数(除非除法是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 ^64,因此乘法后的误差小于1/5.除以4后,误差小于(1/5)*2 ^-2.

(1/5)*2 ^-2 <1/5所以答案总是等于做一个精确的除法并向零舍入.

不幸的是,这对所有除数都不起作用.

如果我们试图将4/7表示为0.64定点数,并且从零开始舍入,则最终得到(6/7)*2 ^-64的误差.乘以一个略低于2 ⁶⁴的i值后,我们最终得到一个不到6/7的误差,除以4后,我们最终得到的误差略低于1.5/7,大于1/7.

因此,为了正确地实现7,我们需要乘以0.65的固定点数.我们可以通过乘以固定点数的低64位来实现,然后加上原始数字(这可能会溢出到进位)然后通过进位进行旋转.

Answer 4

这里是一个算法文档的链接,它生成我在Visual Studio中看到的值和代码(在大多数情况下),并且我假设仍然在GCC中用于将变量整数除以常数整数.

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

在文章中,uword有N位,udword有2N位,n = numerator = dividend,d = denominator = divisor,l最初设置为ceil(log2(d)),shpre是pre-shift(在乘法之前使用) )= 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位乘法器的情况).