【发布时间】：2015-08-14 21:49:27
【问题描述】：
当我第一次获得 Haswell 处理器时，我尝试实施 FMA 来确定 Mandelbrot 集。主要算法是这样的：
intn = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn x2 = square(x), y2 = square(y); //square(x) = x*x
floatn r2 = x2 + y2;
booln mask = r2<cut; //booln is in the float domain non integer domain
if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask)
n -= mask
floatn t = x*y; mul2(t); //mul2(t): t*=2
x = x2 - y2 + cx;
y = t + cy;
}
这将确定 n 像素是否在 Mandelbrot 集中。所以对于双浮点它运行超过 4 个像素（floatn = __m256d，intn = __m256i）。这需要 4 次 SIMD 浮点乘法和 4 次 SIMD 浮点加法。
然后我修改它以像这样与 FMA 一起使用
intn n = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn r2 = mul_add(x,x,y*y);
booln mask = r2<cut;
if(!horizontal_or(mask)) break;
add_mask(n,mask);
floatn t = x*y;
x = mul_sub(x,x, mul_sub(y,y,cx));
y = mul_add(2.0f,t,cy);
}
其中 mul_add 调用 _mm256_fmad_pd 和 mul_sub 调用 _mm256_fmsub_pd。此方法使用 4 个 FMA SIMD 操作和两个 SIMD 乘法，这比没有 FMA 的算术操作少了两次。此外，FMA 和乘法可以使用两个端口，而加法只能使用一个。
为了减少我的测试偏差，我放大了一个完全在 Mandelbrot 集中的区域，因此所有值都是 maxiter。在这种情况下使用 FMA 的方法大约快 27%。 这当然是一种改进，但是从 SSE 到 AVX 使我的性能翻了一番，所以我希望使用 FMA 可能再提高两倍。
但后来我找到了this 关于 FMA 的答案，上面写着
fused-multiply-add 指令的重要方面是（实际上）中间结果的无限精度。这有助于提高性能，但不是因为两个操作被编码在一条指令中 - 它有助于提高性能，因为中间结果的几乎无限精度有时很重要，并且在这种级别的普通乘法和加法中恢复非常昂贵精度确实是程序员所追求的。
后面给出了double*double到double-double乘法的例子
high = a * b; /* double-precision approximation of the real product */
low = fma(a, b, -high); /* remainder of the real product */
据此，我得出结论，我实施 FMA 不是最优的，因此我决定实施 SIMD 双双。我根据论文Extended-Precision Floating-Point Numbers for GPU Computation实现了double-double。该纸是用于双浮动的，所以我将其修改为双双。此外，我没有将一个双精度值封装到 SIMD 寄存器中，而是将 4 个双精度值封装到一个 AVX 高寄存器和一个 AVX 低寄存器中。
对于 Mandelbrot 集，我真正需要的是双双乘法和加法。在那篇论文中，这些是df64_add 和df64_mult 函数。
下图显示了software FMA（左）和硬件 FMA（右）的 df64_mult 函数的程序集。这清楚地表明，硬件 FMA 对双倍乘法是一个很大的改进。
那么硬件 FMA 在双双 Mandelbrot 集计算中的表现如何呢？ 答案是，这仅比使用软件 FMA 快 15%。这比我希望的要少得多。 双双 Mandelbrot 计算需要 4 次双双加法和四次双双乘法（x*x、y*y、x*y 和 2*(x*y)）。但是，2*(x*y) multiplication is trivial for double-double 所以这个乘法可以忽略成本。因此，我认为使用硬件 FMA 的改进如此之小的原因是计算以缓慢的双双加法为主（见下面的汇编）。
过去，乘法比加法慢（程序员使用了几种技巧来避免乘法），但在 Haswell 中，情况似乎正好相反。不仅因为 FMA，还因为乘法可以使用两个端口，而加法只能使用一个。
所以我的问题（最后）是：
- 当加法比乘法慢时如何优化？
- 有没有一种代数方法可以改变我的算法以使用更多的乘法
和更少的添加？我知道有方法可以做相反的事情，例如
(x+y)*(x+y) - (x*x+y*y) = 2*x*y使用两次加法来减少一次乘法。 - 有没有办法简化 df64_add 函数（例如使用 FMA）？
如果有人想知道 double-double 方法比 double 慢十倍左右。我认为这还不错，就好像有一个硬件四精度类型一样，它的速度可能至少是双精度类型的两倍，所以我的软件方法比我对硬件的预期慢五倍（如果它存在的话）。
df64_add大会
vmovapd 8(%rsp), %ymm0
movq %rdi, %rax
vmovapd 72(%rsp), %ymm1
vmovapd 40(%rsp), %ymm3
vaddpd %ymm1, %ymm0, %ymm4
vmovapd 104(%rsp), %ymm5
vsubpd %ymm0, %ymm4, %ymm2
vsubpd %ymm2, %ymm1, %ymm1
vsubpd %ymm2, %ymm4, %ymm2
vsubpd %ymm2, %ymm0, %ymm0
vaddpd %ymm1, %ymm0, %ymm2
vaddpd %ymm5, %ymm3, %ymm1
vsubpd %ymm3, %ymm1, %ymm6
vsubpd %ymm6, %ymm5, %ymm5
vsubpd %ymm6, %ymm1, %ymm6
vaddpd %ymm1, %ymm2, %ymm1
vsubpd %ymm6, %ymm3, %ymm3
vaddpd %ymm1, %ymm4, %ymm2
vaddpd %ymm5, %ymm3, %ymm3
vsubpd %ymm4, %ymm2, %ymm4
vsubpd %ymm4, %ymm1, %ymm1
vaddpd %ymm3, %ymm1, %ymm0
vaddpd %ymm0, %ymm2, %ymm1
vsubpd %ymm2, %ymm1, %ymm2
vmovapd %ymm1, (%rdi)
vsubpd %ymm2, %ymm0, %ymm0
vmovapd %ymm0, 32(%rdi)
vzeroupper
ret
【问题讨论】：
-
正如 Agner Fog 在他的一些优化资源中指出的（agner.org/optimize，但我不记得确切的文件或页面），有时 FMA 会使事情变慢。我想我记得那个例子是
x*x + y*y，对于某些英特尔处理器型号，如果实现为 mul-fma，则延迟会使时序比包含更多并行性的 mul-mul-add 更差。 -
@PascalCuoq，到目前为止，我一直在使用IACA，然后查看时间。对于双倍来说，它主要是一个猜谜游戏，因为对于 Mandelbrot 集，FMA 有很多排列。使用 double-double 很清楚在哪里使用它。但主要问题是双双的加法。与乘法相比，它是如此缓慢。可以做什么？
-
一个正确编程的双双加法包含 20 个基本的浮点运算。您可能会遇到使用较少指令的版本，但是当您需要和想要它时，这些指令会精确地失去准确性，这是对数量几乎相同的操作数的有效减法。双倍加法的加速需要硬件改进。最大幅度和最小幅度运算有一点帮助，单次舍入的三输入加法有很大帮助。一些处理器架构提供前者，我不知道有哪个处理器架构提供后者，
-
如果您的准确度比
double高一点，您可以通过对x2hi和y2hi（6 次失败）进行 Knuth 2-sum 来计算x2hi + x2lo + y2hi + y2lo + cx，另一个Knuth 2-sum 与cx（另外 6 次失败），然后将四个余数相加。您可以使用基于 FMA 的 double*double -> double double 乘积计算2*x*y+cy，使用 cy 进行 Knuth 2-sum，然后将余项相加。 （我没有对此进行测试或认真分析，所以这个建议可能是垃圾。） -
@Zboson：按幅度的最小值和最大值允许您使用“Fast2Sum”；如果
x的幅度大于y，你可以做hi = x+y; lo = (x-hi)+y；lo的计算中的差和和是精确的。 “浮点算术手册”包含所有这些技巧以及更多内容。我发现 Ogita、Rump 和 Oishi 的“精确和和点积”很有教育意义，但它绝不是 Fast2Sum 的原始参考。
标签： assembly x86 floating-point fma double-double-arithmetic