当我第一次得到一个 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 可能会再增加两倍。
但后来我找到了关于 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 寄存器中打包一个 double-double 值，而是将 4 个 double-double 值打包到一个 AVX 高寄存器和一个 AVX 低寄存器中。
对于 Mandelbrot 集，我真正需要的是双倍乘法和加法。在那篇论文中，这些是df64_add和df64_mult函数。下图显示了我df64_mult的软件 FMA（左）和硬件 FMA（右）功能的程序集。这清楚地表明，硬件 FMA 是对双倍乘法的一大改进。
那么硬件FMA在双双Mandelbrot集计算中表现如何呢？答案是这仅比使用软件 FMA 快 15%。这比我希望的要少得多。双双 Mandelbrot 计算需要 4 次双双加法和四次双双乘法（x*x、y*y、x*y和2*(x*y)）。但是，2*(x*y)对于 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