c-用手在浮点除法上减少强度

Question

在本学期计算机科学的最后一项作业中，我们必须对一些代码片段进行强度折减。其中大多数都是直截了当的，尤其是在查看编译器输出时。但其中一个问题我无法解决，即使有编译器的帮助。

我们的教授给了我们以下提示：

提示：查询IEEE 754单精度浮点数在存储器中表示。

以下是代码片段：(a属于double*类型)

for (int i = 0; i < N; ++i) {
a[i] += i / 5.3;   
}

起初，我试图查看这个被剪断的godbolt的编译器输出。我试图在没有任何优化的情况下编译它：(注意：我只复制了for循环中的相关部分)

mov     eax, DWORD PTR [rbp-4]
cdqe
lea     rdx, [0+rax*8]
mov     rax, QWORD PTR [rbp-16]
add     rax, rdx
movsd   xmm1, QWORD PTR [rax]
cvtsi2sd        xmm0, DWORD PTR [rbp-4]    //division relevant
movsd   xmm2, QWORD PTR .LC0[rip]          //division relevant
divsd   xmm0, xmm2                         //division relevant
mov     eax, DWORD PTR [rbp-4]
cdqe
lea     rdx, [0+rax*8]
mov     rax, QWORD PTR [rbp-16]
add     rax, rdx
addsd   xmm0, xmm1
movsd   QWORD PTR [rax], xmm0

和-O3:

.L2:
pshufd  xmm0, xmm2, 238             //division relevant
cvtdq2pd        xmm1, xmm2          //division relevant
movupd  xmm6, XMMWORD PTR [rax]
add     rax, 32
cvtdq2pd        xmm0, xmm0          //division relevant
divpd   xmm1, xmm3                  //division relevant
movupd  xmm5, XMMWORD PTR [rax-16]
paddd   xmm2, xmm4
divpd   xmm0, xmm3                  //division relevant
addpd   xmm1, xmm6
movups  XMMWORD PTR [rax-32], xmm1
addpd   xmm0, xmm5
movups  XMMWORD PTR [rax-16], xmm0
cmp     rax, rbp
jne     .L2

我对程序集代码的除法部分进行了注释。但是这个输出并不能帮助我理解如何对代码片段应用强度降低。(可能有太多的优化正在进行，无法完全理解输出)

其次，我试图理解浮点部分5.3的位表示。即：

0   |10000001|01010011001100110011010
Sign|Exponent|Mantissa

但这对我也没有帮助。

Answer 1

如果我们采用维基百科的定义，

强度降低是一种编译器优化，其中昂贵的操作被等价但较低成本的操作取代

然后我们可以在这里通过将昂贵的浮点除法转换为浮点乘加两个浮点乘加(FMA)来应用强度降低。假设double映射到IEEE-754binary64，浮点计算的默认舍入模式是四舍五入到最近或偶数，并且int是32位类型，我们可以通过简单的穷举测试来证明转换的正确性：

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
int main (void)
{
const double rcp_5p3 = 1.0 / 5.3; // 0x1.826a439f656f2p-3
int i = INT_MAX;
do {
double ref = i / 5.3;
double res = fma (fma (-5.3, i * rcp_5p3, i), rcp_5p3, i * rcp_5p3);
if (res != ref) {
printf ("error: i=%2d  res=%23.13a  ref=%23.13an", i, res, ref);
return EXIT_FAILURE;
}
i--;
} while (i >= 0);
return EXIT_SUCCESS;
}

大多数常见处理器体系结构(如x86-64和ARM64)的现代实例都支持FMA的硬件，因此fma()可以直接映射到适当的硬件指令。这应该通过查看生成的二进制文件的反汇编来确认。在缺乏对FMA的硬件支持的情况下，显然不应该应用转换，因为fma()的软件实现很慢，有时功能不正确。

这里的基本思想是，在数学上，除法等价于倒数的乘法。然而，对于有限精度浮点运算，这不一定是真的。上面的代码试图通过在FMA的帮助下确定朴素方法中的误差并在必要时应用校正来提高比特精确计算的可能性。有关背景(包括参考文献)，请参阅前面的问题。

据我所知，目前还没有一个通用的数学证明算法来确定上面的变换对于哪些除数与哪些除数配对是安全的(也就是说，提供了比特精确的结果)，这就是为什么穷举测试对于证明变换是有效的是绝对必要的。

Pascal Cuoq在评论中指出，有一种替代算法可以通过将除数的倒数预先计算到超过本机精度的精度，特别是双倍精度，来利用编译时常数除数来潜在地减少浮点除法的强度。背景参见N.Brisebarre和J.-M.M.Muller；用阿贝里精度常数"正确取整乘法">IEEE计算机汇刊，57(2)：162-1742008年2月，它还提供了如何确定这种转换对任何特定常数是否安全的指导。由于目前的情况很简单，我再次使用详尽的测试来证明它是安全的。在适用的情况下，这将把除法减少到一个FMA加一个乘法：

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <mathimf.h>
int main (void)
{
const double rcp_5p3_hi =  1.8867924528301888e-1; //  0x1.826a439f656f2p-3
const double rcp_5p3_lo = -7.2921377017921457e-18;// -0x1.0d084b1883f6e0p-57
int i = INT_MAX;
do {
double ref = i / 5.3;
double res = fma (i, rcp_5p3_hi, i * rcp_5p3_lo);
if (res != ref) {
printf ("i=%2d  res=%23.13a  ref=%23.13an", i, res, ref);
return EXIT_FAILURE;
}
i--;
} while (i >= 0);
return EXIT_SUCCESS;
}

Answer 2

为了涵盖另一个方面：由于类型为int的所有值都可以精确地表示为double(而不是float)，因此可以通过引入从0.0到N:的浮点变量来消除循环中计算i / 5.3时发生的int到double转换

double fp_i = 0;
for (int i = 0; i < N; fp_i += 1, i++)
a[i] += fp_i / 5.3;

然而，这会扼杀自动矢量化，并引入一系列相关的浮点加法。浮点加法通常是3或4个周期，因此最后一次迭代将在至少(N-1)*3个周期后退出，即使CPU可以更快地调度循环中的指令。值得庆幸的是，浮点除法并不是完全流水线化的，x86 CPU调度浮点除法的速率大致匹配或超过加法指令的延迟。

这就留下了终止矢量化的问题。可以通过手动展开循环并引入两个独立的链来将其恢复，但对于AVX，您需要四个链来进行完全矢量化：

double fp_i0 = 0, fp_i1 = 1;
int i = 0;
for (; i+1 < N; fp_i0 += 2, fp_i1 += 2, i += 2) {
double t0 = a[i], t1 = a[i+1];
a[i]   = t0 + fp_i0 / 5.3;
a[i+1] = t1 + fp_i1 / 5.3;
}
if (i < N)
a[i] += i / 5.3;

Answer 3

CAVEAT：几天后，我意识到这个答案是不正确的，因为它忽略了o / 5.3计算中下溢(低于正常值或为零)的后果。在这种情况下，将结果乘以2的幂是"精确的"，但不会产生较大整数除以5.3的结果

CCD_ 15只需要针对CCD_。对于i的偶数值，您可以简单地将(i/2)/5.3的值乘以2.0，该值在循环的前面已经计算过了。

剩下的困难是重新排序迭代，使得0和N-1之间的每个索引只处理一次，并且程序不需要记录任意数量的除法结果。

实现这一点的一种方法是对小于N的所有奇数o进行迭代，并且在计算o / 5.3以便处理索引o之后，还处理形式为o * 2**p的所有索引。

if (N > 0) {
a[0] += 0.0; // this is needed for strict IEEE 754 compliance lol
for (int o = 1; o < N; o += 2) {
double d = o / 5.3;
int i = o;
do { 
a[i] += d;
i += i;
d += d;
} while (i < N);
}
}

注意：这不使用提供的提示"查询IEEE 754单精度浮点数在内存中的表示方式"。我想我很清楚单精度浮点数是如何在内存中表示的，但我不知道这有什么关系，尤其是因为代码中没有单精度值或计算需要优化。我认为这个问题的表达方式有一个错误，但从技术上讲，以上仍然是对所用问题的部分回答。

在上面的代码片段中，我还忽略了接近INT_MAX的N值的溢出问题，因为代码已经足够复杂了。

另外需要注意的是，上面的转换只替换了两个除法中的一个。它通过使代码不可分解(并且对缓存不太友好)来实现这一点。在你的问题中，gcc -O3已经表明，自动矢量化可以应用于你的教授建议的起点，这可能比抑制一半的除法更有益。这个答案中的转换唯一的好处是，它是一种强度降低，这是你的教授要求的。