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

我一直在阅读 divmul 汇编操作,我决定通过在C中编写一个简单的程序来看它们的运行情况:

文件division.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 文件,它并没有't contain any div operations! Instead, it does some kind of black magic with bit shifting and magic numbers. Here'是一个计算 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?它如何产生这个神奇的数字以及为什么一切都有效?

回答(4)

2 years ago

整数除法是您可以在现代处理器上执行的最慢的算术运算之一,延迟可达数十个周期和吞吐量不佳 . (对于x86,请参阅Agner Fog's instruction tables and microarch guide) .

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

以这种方式实现C / 运算符而不是涉及 div 的多指令序列只是GCC 's default way of doing division by constants. It doesn' t需要跨操作进行优化,即使调试也不会改变任何内容 . (但是对于小代码大小使用 -Os 确实让GCC使用 div . )使用乘法逆而不是除法就像使用 lea 而不是 muladd

因此,如果在编译时未知除数,则只会在输出中看到 dividiv .

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

2 years ago

除以5与乘以1/5相同,再乘以4/5并向右移2位相同 . 有关的值是十六进制的 CCCCCCCCCCCCD ,如果放在一个十六进制点之后是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... 等于1 . 因此,我们知道 x + x/4 = 1 ,所以 5x/4 = 1x=4/5 . 然后将其表示为十六进制的 CCCCCCCCCCCCD 用于舍入(因为超出最后一个的二进制数字将是 1 ) .

2 years ago

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

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

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


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

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

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

2 years ago

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

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

在本文中,uword有N位,udword有2N位,n =分子,d =分母= 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位乘法器的公式 .

回到图6.2 . 只有当除数> 2 ^(N-1)(当ℓ== N)时,分子才能大于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 numerator = 2^(ℓ) = (upper part of 2^(N+ℓ))
; lower 8 bytes of numerator for mlow  = 0
; lower 8 bytes of numerator for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e)
numerator dq    2 dup(?)        ;16 byte numerator
divisor   dq    1 dup(?)        ; 8 byte divisor
; ...
        mov     rcx,divisor
        mov     rdx,0
        mov     rax,numerator+8    ;upper 8 bytes of numerator
        div     rcx                ;after div, rax == 1
        mov     rax,numerator      ;lower 8 bytes of numerator
        div     rcx
        mov     rdx,1              ;rdx:rax = N+1 bit value = 65 bit value

您可以使用GCC进行测试 . 你已经看到了如何处理j = i / 5 . 看看如何处理j = i / 7(应该是N 1位乘法器的情况) .