为什么在强度降低乘法到循环携带加法之后，这段代码执行得更慢？
我正在阅读 Agner Fog 的 optimization manuals，遇到了这个例子：
double data[LEN];
void compute()
{
const double A = 1.1, B = 2.2, C = 3.3;
int i;
for(i=0; i<LEN; i++) {
data[i] = A*i*i + B*i + C;
}
}
Agner 指出有一种方法可以优化此代码 - 通过意识到循环可以避免使用代价高昂的乘法，而是使用每次迭代应用的“增量”。
我用一张纸来证实这个理论，首先...
https://i.stack.imgur.com/wzzgt.png
...当然，他是对的——在每次循环迭代中，我们可以通过添加“delta”来基于旧结果计算新结果。这个增量从值“A+B”开始，然后在每一步增加“2*A”。
因此，我们将代码更新为如下所示：
void compute()
{
const double A = 1.1, B = 2.2, C = 3.3;
const double A2 = A+A;
double Z = A+B;
double Y = C;
int i;
for(i=0; i<LEN; i++) {
data[i] = Y;
Y += Z;
Z += A2;
}
}
就操作复杂度而言，这两个版本的功能差异确实是惊人的。与加法相比，乘法在我们的 CPU 中速度明显较慢而闻名。我们已经替换了 3 次乘法和 2 次加法……只有 2 次加法！
所以我继续添加一个循环来执行 compute 很多次 - 然后保持执行所需的最短时间：
unsigned long long ts2ns(const struct timespec *ts)
{
return ts->tv_sec * 1e9 + ts->tv_nsec;
}
int main(int argc, char *argv[])
{
unsigned long long mini = 1e9;
for (int i=0; i<1000; i++) {
struct timespec t1, t2;
clock_gettime(CLOCK_MONOTONIC_RAW, &t1);
compute();
clock_gettime(CLOCK_MONOTONIC_RAW, &t2);
unsigned long long diff = ts2ns(&t2) - ts2ns(&t1);
if (mini > diff) mini = diff;
}
printf("[-] Took: %lld ns.\n", mini);
}
我编译了这两个版本，运行它们......然后看到这个：
gcc -O3 -o 1 ./code1.c
gcc -O3 -o 2 ./code2.c
./1
[-] Took: 405858 ns.
./2
[-] Took: 791652 ns.
嗯，这是出乎意料的。由于我们报告了最短执行时间，我们正在丢弃由操作系统的各个部分引起的“噪音”。我们还注意在一台完全不做任何事情的机器上运行。结果或多或少是可重复的 - 重新运行两个二进制文件表明这是一个一致的结果：
for i in {1..10} ; do ./1 ; done
[-] Took: 406886 ns.
[-] Took: 413798 ns.
[-] Took: 405856 ns.
[-] Took: 405848 ns.
[-] Took: 406839 ns.
[-] Took: 405841 ns.
[-] Took: 405853 ns.
[-] Took: 405844 ns.
[-] Took: 405837 ns.
[-] Took: 406854 ns.
for i in {1..10} ; do ./2 ; done
[-] Took: 791797 ns.
[-] Took: 791643 ns.
[-] Took: 791640 ns.
[-] Took: 791636 ns.
[-] Took: 791631 ns.
[-] Took: 791642 ns.
[-] Took: 791642 ns.
[-] Took: 791640 ns.
[-] Took: 791647 ns.
[-] Took: 791639 ns.
接下来唯一要做的就是查看编译器为两个版本中的每一个创建了什么样的代码。
objdump -d -S 表明 compute 的第一个版本 - “哑”，但不知何故快代码 - 有一个看起来像这样的循环：
https://i.stack.imgur.com/JbWNJ.png
那么第二个优化版本呢——它只做了两个添加？
https://i.stack.imgur.com/HQwXx.png
现在我不了解你，但就我自己而言，我……很困惑。第二个版本的指令减少了大约 4 倍，其中两个主要的只是基于 SSE 的添加 (addsd)。第一个版本不仅有 4 倍的指令……它还包含（如预期的那样）乘法（mulpd）。
我承认我没想到会有这样的结果。不是因为我是 Agner 的粉丝（我是，但这无关紧要）。
知道我缺少什么吗？我在这里犯了什么错误，可以解释速度的差异吗？请注意，我已经对 Xeon W5580 和 Xeon E5-1620 进行了测试——在这两个版本中，第一个（哑）版本比第二个版本快得多。
为了便于重现结果，两个版本的代码有两个要点：Dumb yet somehow faster 和 optimized, yet somehow slower。
PS请不要评论浮点精度问题；这不是这个问题的重点。
4*A2 或类似的增量。可能 clang 可以使用 -ffast-math 为您做到这一点（甚至可能是 GCC，但 GCC 往往会在没有多个累加器的情况下展开。）在 Haswell 或更高版本上使用 FMA 时，Horner 的方法对于如此短的多项式非常有用，易于输出of-order exec 要隐藏，尽管它仍然需要 i 的 FP 版本
significand = sig1 * sig2; exponent = exp1+exp2），而对于浮点加法，它需要串行完成（确定结果指数，然后将两个值“移位”到匹配结果指数，然后确定结果有效位）。
addpd（和 vfma...）的 mulpd 相比，Alder Lake 将 addpd 延迟从 4 个周期降低到 3 个周期，后者是自 Skylake 以来 addpd/subpd/mulpd/vfma...pd 的延迟. AMD 在某些 CPU 上的附加值较低，但 Zen2 具有 3 周期延迟 addpd 和 mulpd 与 5c fma 相比，如 Broadwell
了解您所看到的性能差异的关键在于矢量化。是的，基于加法的解决方案在其内部循环中只有两条指令，但重要的区别不在于循环中有多少指令，而在于每条指令执行了多少工作。
在第一个版本中，输出完全取决于输入：每个 data[i] 都是 i 本身的一个函数，这意味着每个 data[i] 可以按任何顺序计算：编译器可以向前、向后执行它们，横向，无论如何，您仍然会得到相同的结果 - 除非您从另一个线程观察该内存，否则您永远不会注意到数据正在以哪种方式被处理。
在第二个版本中，输出不依赖于 i — 它依赖于上次循环中的 A 和 Z。
如果我们将这些循环的主体表示为小的数学函数，它们将具有非常不同的整体形式：
f(i) -> 迪
f(Y, Z) -> (di, Y', Z')
在后一种形式中，对 i 没有实际的依赖关系——计算函数值的唯一方法是从函数的最后一次调用中知道前面的 Y 和 Z，这意味着函数形成一个链条——在你完成前一个之前，你不能做下一个。
为什么这很重要？因为 CPU 有 vector parallel 指令，每个 可以同时执行两个、四个甚至八个算术运算！ （AVX CPU 可以并行执行更多操作。）即 4 次乘法、4 次加法、4 次减法、4 次比较 — 4 次任意操作！因此，如果您尝试计算的输出仅取决于输入，那么您可以安全地一次计算两个、四个甚至八个——它们是否向前并不重要或向后，因为结果是相同的。但是，如果输出依赖于先前的计算，那么你就会被困在串行形式中——一次一个。
这就是为什么“更长”的代码会赢得性能。尽管它有更多的设置，并且实际上做更多的工作，但大部分工作都是并行完成的：它不只是在循环的每次迭代中计算 data[i] - 它是同时计算 data[i]、data[i+1]、data[i+2] 和 data[i+3]，然后跳到下一组四个。
为了扩展我在这里的意思，编译器首先将原始代码变成了这样的东西：
int i;
for (i = 0; i < LEN; i += 4) {
data[i+0] = A*(i+0)*(i+0) + B*(i+0) + C;
data[i+1] = A*(i+1)*(i+1) + B*(i+1) + C;
data[i+2] = A*(i+2)*(i+2) + B*(i+2) + C;
data[i+3] = A*(i+3)*(i+3) + B*(i+3) + C;
}
如果你眯着眼睛看，你可以说服自己会做和原版一样的事情。之所以这样做是因为所有这些相同的垂直运算符行：所有这些 * 和 + 操作都是相同的操作，只是在不同的数据上执行 - CPU 具有特殊的内置指令，可以执行多个* 或多个 + 在同一时间对不同数据进行操作，每个操作仅在一个时钟周期内。
请注意较快解决方案中说明中的字母 p — addpd 和 mulpd — 以及较慢解决方案中说明中的字母 s — addsd。那是“添加包装双打”和“乘以包装双打”与“添加单双打”。
不仅如此，看起来编译器也部分展开了循环——循环不只是在每次迭代中执行 两个 值，实际上是 四个，并将操作交错到避免依赖和停顿，所有这些都减少了汇编代码必须测试 i < 1000 的次数。
但是，所有这些都只有在循环的迭代之间没有依赖关系时才有效：如果唯一决定每个 data[i] 发生什么的事情是 i 本身。如果存在依赖关系，如果上一次迭代的数据影响下一次迭代，那么编译器可能会受到它们的限制，以至于它根本无法更改代码——而不是编译器能够使用花哨的并行指令或聪明的优化（CSE、strength reduction、loop unrolling、重新排序等），你得到的代码正是你输入的代码——添加 Y，然后添加 Z，然后重复。
但是在这里，在代码的第一个版本中，编译器正确地识别出数据中没有依赖关系，并发现它可以并行完成工作，所以它确实做到了，这就是一切的不同之处。
这里的主要区别是循环依赖。第二种情况下的循环是依赖的——循环中的操作依赖于之前的迭代。这意味着每次迭代甚至不能在前一次迭代完成之前开始。在第一种情况下，循环体是完全独立的——循环体中的所有内容都是自包含的，仅取决于迭代计数器和常量值。这意味着循环可以并行计算——多次迭代可以同时进行。然后，这允许循环被简单地展开和矢量化，重叠许多指令。
如果您要查看性能计数器（例如，使用 perf stat ./1），您会看到第一个循环除了运行速度更快之外，每个周期还运行更多指令 (IPC)。相比之下，第二个循环有更多的依赖周期——CPU 无所事事，等待指令完成，然后才能发出更多指令的时间。
第一个可能会成为内存带宽的瓶颈，特别是如果您让编译器在 Sandy Bridge (gcc -O3 -march=native) 上使用 AVX 自动矢量化（如果它设法使用 256 位向量）。此时，IPC 将下降，尤其是对于 L3 cache 而言太大的输出数组。
注意：展开和矢量化不需要独立的循环——当（某些）循环依赖存在时，您可以执行它们。然而，它更难，回报也更少。因此，如果您想从矢量化中看到最大的加速，它有助于在可能的情况下消除循环依赖。
perf stat --all-user ./1 之类的东西来累积整个程序的计数。这很好，因为它大部分时间都在循环中。您可能希望将时间移到循环之外或将其删除以进行这种分析，也许通过将实际工作放在 __attribute__((noinline,noipa)) 函数中来对优化器隐藏重复循环，以停止过程间分析和内联。
Z0 += 8*A2（或16*A2 如果每个 Z 向量包含 4 个双精度数而不是 2)。你需要一些数学来解释一个元素跨步 8 或 16 个 i 值而不是 1，也许在某个地方有一个乘法。您几乎可以肯定比每次迭代都重做 int->FP 转换做得更好；这是一种获得独立迭代的昂贵方式。
这种 method of finite differences 强度降低优化可以比您可以为每个 i 单独重新评估多项式的最佳速度提供加速。但只有当你将它推广到更大的步幅时，循环中仍然有足够的并行度。 我的版本在我的 Skylake 上每个时钟周期存储一个向量（四个双精度数），用于适合 L1d 缓存的小数组；否则就是带宽测试。在早期的英特尔上，它还应该最大限度地提高 SIMD FP-add 吞吐量，包括带有 AVX 的 Sandy Bridge（如果您对齐输出，则每两个时钟有 1x 256 位添加/时钟和 1x 256 位存储。）
对上一次迭代中的值的依赖是致命的
这个 strength-reduction optimization（只是相加而不是从一个新的 i 开始并相乘）在循环迭代中引入了串行依赖，涉及 FP 数学而不是整数增量。
原始版本具有每个输出元素之间的数据并行性：每个元素仅依赖于常量和它自己的 i 值。编译器可以使用 SIMD 自动矢量化（SSE2，如果使用 -O3 -march=native，则为 AVX），并且 CPU 可以在循环迭代中重叠工作并无序执行。尽管有很多额外的工作，但 CPU 能够在编译器的帮助下应用足够的蛮力。
但是根据 poly(i) 计算 poly(i+1) 的版本具有非常有限的并行性；没有 SIMD 矢量化，您的 CPU 每四个周期只能运行两个标量加法，例如，其中四个周期是英特尔 Skylake 通过 Tiger Lake 的 FP 加法的延迟。 (https://uops.info/)。
huseyin tugrul buyukisik's answer 展示了如何在更现代的 CPU 上接近最大化原始版本的吞吐量，使用两个 FMA 运算来评估多项式 (Horner's scheme)，以及一个 int 到浮点转换或一个浮点点增量。 （后者创建一个 FP 依赖链，您需要展开以隐藏它。）
所以最好的情况是每个输出向量有三个浮点数学运算。 （加上一家商店）。当前的 Intel CPU 只有两个浮点执行单元，可以运行 FP 数学运算，包括 int-to-double。 （对于 512 位向量，当前 CPU 会关闭端口 1 上的向量 ALU，因此根本只有两个 SIMD ALU 端口，因此非 FP 数学运算，如 SIMD 整数增量，也会竞争 SIMD 吞吐量。除了对于只有一个 512 位 FMA 单元的 CPU，则端口 5 可用于其他工作。）
AMD 因为 Zen 2 在两个端口上有两个 FMA/mul 单元，在两个不同的端口上有两个 FP add/sub 单元，所以如果你使用 FMA 进行加法，理论上每个时钟周期最多有四个 SIMD 加法。
Haswell/Broadwell 有 2 个/时钟 FMA，但只有 1 个/时钟 FP add/sub（延迟较低）。这适用于简单代码，not great 适用于经过优化以具有大量并行性的代码。这可能就是英特尔在 Skylake 中改变它的原因。 （Alder Lake 重新引入了低延迟 FP add/sub，但 2/clock 吞吐量与乘法相同。有趣的是，非交换 latency：目标只有 2 个周期，另一个操作数只有 3 个周期，所以它适用于更长的时间依赖链。）
您的 Sandy Bridge (E5-1620) 和 Nehalem (W5580) CPU 在不同的端口上有 1/clock FP add/sub、1/clock FP mul。这就是 Haswell 的基础。以及为什么添加额外的乘法不是一个大问题：它们可以与现有的添加并行运行。 （Sandy Bridge 的宽度为 256 位，但您在未启用 AVX 的情况下编译：使用 -march=native。）
寻找并行性：以任意步幅降低强度
您的 compute2 根据前一个值计算下一个 Y 和下一个 Z。即，步幅为 1，您需要 data[i+1] 的值。所以每次迭代都依赖于前一个迭代。
如果将其推广到其他步幅，则可以提前 4、6、8 或更多个单独的 Y 和 Z 值，以便它们彼此同步跨越，彼此独立。这为编译器和/或 CPU 重新获得了足够的并行性来利用。
poly(i) = A i^2 + B i + C
poly(i+s) = A (i+s)^2 + B (i+s) + C
= A*i^2 + A*2*s*i + A*s^2 + B*i + B*s + C
= poly(i) + A*2*s*i + A*s^2 + B*s + C
所以这有点混乱，不完全清楚如何将其分解为 Y 和 Z 部分。 （这个答案的早期版本弄错了。）
可能更容易从 FP 值序列 (Method of Finite Differences) 的步幅的一阶和二阶差异向后工作。这将直接找到我们需要添加的内容； Z[] 初始化器和步骤。
这基本上就像取一阶和二阶导数，然后优化循环有效地积分以恢复原始函数。以下输出由以下基准中 main 的正确性检查部分生成。
# method of differences for stride=1, A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+1)
0 1
1 3 2 # 4-1 = 3 | 3-1 = 2
4 5 2 # 9-4 = 5 | 5-3 = 2
9 7 2 # ...
16 9 2
25 11 2
相同的多项式 (x^2)，但步幅为 3。非 2 的幂有助于显示步幅的因子/幂来自何处，与自然发生的因子 2 相比。
# for stride of 3, printing in groups. A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+3)
0 9
1 15
4 21
9 27 18 # 36- 9 = 27 | 27-9 = 18
16 33 18 # 49-16 = 33 | 33-15 = 18
25 39 18 # ...
36 45 18 # 81-36 = 45 | 45-27 = 18
49 51 18
64 57 18
81 63 18
100 69 18
121 75 18
Y[] 和 Z[] 初始化器
初始 Y[j] = poly(j) 因为它必须存储到相应位置的输出（data[i+j] = Y[j]）。
初始 Z[j] 将被添加到 Y[j] 中，并且需要将其变为 poly(j+stride)。因此，初始 Z[j] = poly(j+stride) - Y[j]，如果我们愿意，我们可以对其进行代数化简。 （对于编译时常数 A、B、C，编译器将以任何一种方式进行常数传播。） Z[j] 保存通过 poly(x) 的一阶差异，用于 poly(0..stride) 的起点-1)。这是上表的中间列。
对 Z[j] += second_difference 的必要更新是一个标量常数，正如我们从二阶差异中看到的那样。通过使用几个不同的步幅和 A 值（i^2 的系数），我们可以看到它是 A * 2 * (stride * stride)。 （使用像 3 和 5 这样的非互质值有助于解开事物。）使用更多的代数，您可以象征性地展示这一点。从微积分 PoV 来看，因子 2 是有意义的：d(A*x^2)/dx = 2Ax，二阶导数是 2A。
// Tested and correct for a few stride and coefficient values.
#include <stdalign.h>
#include <stdlib.h>
#define LEN 1024
alignas(64) double data[LEN];
//static const double A = 1, B = 0, C = 0; // for easy testing
static const double A = 5, B = 3, C = 7; // can be function args
void compute2(double * const __restrict__ data)
{
const int stride = 16; // unroll factor. 1 reduces to the original
const double diff2 = (stride * stride) * 2 * A; // 2nd-order differences
double Z[stride], Y[stride];
for (int j = 0 ; j<stride ; j++){ // this loop will fully unroll
Y[j] = j*j*A + j*B + C; // poly(j) starting values to increment
//Z[j] = (j+stride)*(j+stride)*A + (j+stride)*B + C - Y[j];
//Z[j] = 2*j*stride*A + stride*stride*A + stride*B;
Z[j] = ((2*j + stride)*A + B)*stride; // 1st-difference to next Y[j], from this to the next i
}
for(ptrdiff_t i=0; i < LEN - (stride-1); i+=stride) {
// loops that are easy(?) for a compiler to roll up into some SIMD vectors
for (int j=0 ; j<stride ; j++) data[i+j] = Y[j]; // store
for (int j=0 ; j<stride ; j++) Y[j] += Z[j]; // add
for (int j=0 ; j<stride ; j++) Z[j] += diff2; // add
}
// cleanup for the last few i values
for (int j = 0 ; j < LEN % stride ; j++) {
// let the compiler see LEN%stride to help it decide *not* to auto-vectorize this part
//size_t i = LEN - (stride-1) + j;
//data[i] = poly(i);
}
}
对于 stride=1（不展开），这些简化为原始值。但是对于较大的 stride，编译器可以将 Y[] 和 Z[] 的元素分别保存在几个 SIMD 向量中，因为每个 Y[j] 只与相应的 Z[j] 交互。
有 stride 个独立的并行依赖链供编译器 (SIMD) 和 CPU（流水线执行单元）利用，运行速度比原始速度快 stride 倍，直至达到 SIMD FP-add 吞吐量的瓶颈如果您的缓冲区不适合 L1d，而不是延迟或存储带宽。 （或者直到编译器面对并且不展开和矢量化这些循环的地步！）
这在实践中是如何编译的：非常适合 clang
(Godbolt compiler explorer) 使用 clang14 -O3 -march=skylake -ffast-math 可以很好地使用 stride=16（4 个 YMM 向量，每个 4 个 double）进行 Clang 自动矢量化。
看起来 clang 进一步展开了 2，将 Z[j] += diff2 快捷方式转换为 tmp = Z[j] + diff2; / Z[j] += 2*diff2;。这减轻了 Z 依赖链上的压力，只留下 Y[j] 来应对 Skylake 的延迟瓶颈。
所以每个 asm 循环迭代执行 2x 8 vaddpd 指令和 2x 4 存储。循环开销是 add + 宏融合 cmp/jne，所以 2 uops。 （或者对于全局数组，只有一个 add/jne uop，将负索引向上计数到零；它的索引相对于数组的末尾。）
Skylake 以每时钟周期几乎 1 个商店和 2x vaddpd 的速度运行它。这是这两件事的最大吞吐量。前端只需要跟上 3 uops / 时钟周期多一点，但自 Core2 以来它一直是 4 宽。 Sandy Bridge 系列中的 uop 缓存使这没有问题。 （除非您在 Skylake 上遇到 JCC 勘误表，所以我使用了 -mbranches-within-32B-boundaries to have clang pad instructions to avoid that。）
由于 Skylake 的 vaddpd 延迟为 4 个周期，来自 stride=16 的 4 个依赖链仅足以保持 4 个独立操作的运行。任何时候 Y[j]+= 没有运行它准备好的循环，就会产生一个气泡。由于 Clang 对 Z[] 链的额外展开，Z[j]+= 可以提前运行，因此 Z 链可以领先。显然，使用最旧的就绪优先调度，它倾向于稳定到 Yj+= uops 没有冲突的状态，因为它确实在我的 Skylake 上全速运行。如果我们可以让编译器仍然为 stride=32 生成漂亮的 asm，那将留下更多空间，但不幸的是它没有。 （以对奇数尺寸进行更多清理工作为代价。）
奇怪的是，Clang 只使用 -ffast-math 对其进行矢量化。以下完整基准测试中的模板版本不需要 -ffast-math。源代码经过精心编写，对 SIMD 友好，并按源代码顺序进行数学运算。 （不过，快速数学允许 clang 展开更多的 Z 增量。）
编写循环的另一种方法是使用一个内部循环而不是所有 Y 操作，然后是所有 Z 操作。这在下面的基准测试中很好（有时实际上更好），但即使使用 -ffast-math，它也不会矢量化。从编译器中为像这样的非平凡问题获取最佳展开的 SIMD 汇编可能是繁琐且不可靠的，并且可能需要一些时间。
我将它包含在 Godbolt 的 #if 0 / #else / #endif 块中。
// can auto-vectorize better or worse than the other way
// depending on compiler and surrounding code.
for(int i=0; i < LEN - (stride-1); i+=stride) {
for (int j = 0 ; j<stride ; j++){
data[i+j] = Y[j];
Y[j] += Z[j];
Z[j] += deriv2;
}
}
我们必须手动选择适当的展开量。太大的展开因子甚至会阻止编译器查看正在发生的事情并停止将临时数组保存在寄存器中。例如，32 或 24 是 clang 的问题，但不是 16。可能有一些调整选项可以强制编译器将循环展开到一定数量；有些 GCC 有时可以用来让它在编译时看穿某些东西。
另一种方法是使用 #include <immintrin.h> 和 __m256d Z[4] 而不是 double Z[16] 进行手动矢量化。但是这个版本可以为其他 ISA 向量化，比如 AArch64。
当问题大小不是展开的倍数时，大展开因子的其他缺点是留下更多的清理工作。 （您可以使用 compute1 策略进行清理，让编译器在执行标量之前将其向量化以进行一两次迭代。）
理论上，编译器将允许使用 -ffast-math 为您执行此操作，或者从 compute1 对原始多项式进行强度减少，或者从 compute2 查看步幅如何累积。
但在实践中，这真的很复杂，人类必须自己做一些事情。除非/直到有人开始教编译器如何寻找这样的模式并应用差异方法本身，并选择跨步！但是，即使使用 -ffast-math，也可能不希望对具有不同误差累积属性的算法进行大规模重写。 （整数不会有任何精度问题，但它仍然是一个复杂的模式匹配/替换。）
实验性能结果：
我使用 clang13.0.0 在我的桌面（i7-6700k）上进行了测试。实际上，这确实在每个时钟周期运行 1 个 SIMD 存储，其中有几种编译器选项（快速数学或非快速数学）和内部循环策略上的 #if 0 与 #if 1 的组合。我的基准/测试框架基于@huseyin tugrul buyukisik 的版本，经过改进以在 rdtsc 指令之间重复更可测量的数量，并使用测试循环来检查多项式的简单计算的正确性。
我还让它补偿了核心时钟频率和 "reference" frequency of the TSC read by rdtsc 之间的差异，在我的例子中是 3.9 GHz 与 4008 MHz。 （额定最大 turbo 为 4.2 GHz，但在 Linux 上 EPP = balance_performance，它只希望时钟高达 3.9 GHz。）
源代码 on Godbolt：使用 1 个内部循环，而不是 3 个单独的 j<16 循环，并且不使用 -ffast-math。使用 __attribute__((noinline)) 防止它内联到重复循环中。选项和来源的其他一些变体导致循环内出现一些 vpermpd 洗牌。
下面的基准数据来自以前的版本，它带有一个有缺陷的 Z[j] 初始化程序，但循环 asm 相同。 Godbolt 链接现在在定时循环之后进行了正确性测试，该测试通过了。我的桌面上的实际性能仍然相同，每个 double 仅超过 0.25 个周期，即使没有 #if 1 / -ffast-math 以允许额外展开铿锵声。
$ clang++ -std=gnu++17 -O3 -march=native -mbranches-within-32B-boundaries poly-eval.cpp -Wall
# warning about noipa, only GCC knows that attribute
$ perf stat --all-user -etask-clock,context-switches,cpu-migrations,page-faults,cycles,instructions,uops_issued.any,uops_executed.thread,fp_arith_inst_retired.256b_packed_double -r10 ./a.out
... (10 runs of the whole program, ending with)
...
0.252295 cycles per data element (corrected from ref cycles to core clocks for i7-6700k @ 3.9 GHz)
0.252109 cycles per data element (corrected from ref cycles to core clocks for i7-6700k @ 3.9 GHz)
xor=4303
min cycles per data = 0.251868
Performance counter stats for './a.out' (10 runs):
298.92 msec task-clock # 0.989 CPUs utilized ( +- 0.49% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
129 page-faults # 427.583 /sec ( +- 0.56% )
1,162,430,637 cycles # 3.853 GHz ( +- 0.49% ) # time spent in the kernel for system calls and interrupts isn't counted, that's why it's not 3.90 GHz
3,772,516,605 instructions # 3.22 insn per cycle ( +- 0.00% )
3,683,072,459 uops_issued.any # 12.208 G/sec ( +- 0.00% )
4,824,064,881 uops_executed.thread # 15.990 G/sec ( +- 0.00% )
2,304,000,000 fp_arith_inst_retired.256b_packed_double # 7.637 G/sec
0.30210 +- 0.00152 seconds time elapsed (+- 0.50%)
fp_arith_inst_retired.256b_packed_double 每条 FP add 或 mul 指令计为 1（FMA 为 2），因此每个时钟周期我们得到 1.98 个vaddpd 指令，用于整个程序，包括打印等。这非常接近理论上的最大 2/时钟，显然没有受到次优 uop 调度的影响。 （我增加了重复循环，因此程序将大部分时间都花在了那里，使整个程序上的 perf stat 有用。）
这种优化的目标是用更少的 FLOPS 完成相同的工作，但这也意味着我们基本上在不使用 FMA 的情况下最大化了 Skylake 的 8 FLOP/时钟限制。 （单核上 3.9 GHz 时为 30.58 GFLOP/s）。
非内联函数的汇编（objdump -drwC -Mintel）； clang 使用了 4 个 Y、Z 对 YMM 向量，并将循环进一步展开 3 倍，使其成为 24 KiB 大小的精确倍数，无需清理。请注意 add rax,0x30 每次迭代执行 3 * stride=0x10 双打。
0000000000001440 <void compute2<3072>(double*)>:
# just loading constants; the setup loop did fully unroll and disappear
1440: c5 fd 28 0d 18 0c 00 00 vmovapd ymm1,YMMWORD PTR [rip+0xc18] # 2060 <_IO_stdin_used+0x60>
1448: c5 fd 28 15 30 0c 00 00 vmovapd ymm2,YMMWORD PTR [rip+0xc30] # 2080
1450: c5 fd 28 1d 48 0c 00 00 vmovapd ymm3,YMMWORD PTR [rip+0xc48] # 20a0
1458: c4 e2 7d 19 25 bf 0b 00 00 vbroadcastsd ymm4,QWORD PTR [rip+0xbbf] # 2020
1461: c5 fd 28 2d 57 0c 00 00 vmovapd ymm5,YMMWORD PTR [rip+0xc57] # 20c0
1469: 48 c7 c0 d0 ff ff ff mov rax,0xffffffffffffffd0
1470: c4 e2 7d 19 05 af 0b 00 00 vbroadcastsd ymm0,QWORD PTR [rip+0xbaf] # 2028
1479: c5 fd 28 f4 vmovapd ymm6,ymm4 # buggy Z[j] initialization in this ver used the same value everywhere
147d: c5 fd 28 fc vmovapd ymm7,ymm4
1481: c5 7d 28 c4 vmovapd ymm8,ymm4
1485: 66 66 2e 0f 1f 84 00 00 00 00 00 data16 cs nop WORD PTR [rax+rax*1+0x0]
# top of outer loop. The NOP before this is to align it.
1490: c5 fd 11 ac c7 80 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x180],ymm5
1499: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
149d: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
14a1: c5 fd 11 9c c7 a0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1a0],ymm3
14aa: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
14ae: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
14b2: c5 fd 11 94 c7 c0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1c0],ymm2
14bb: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
14bf: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
14c3: c5 fd 11 8c c7 e0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1e0],ymm1
14cc: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
14d0: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
14d4: c5 fd 11 ac c7 00 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x200],ymm5
14dd: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
14e1: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
14e5: c5 fd 11 9c c7 20 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x220],ymm3
14ee: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
14f2: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
14f6: c5 fd 11 94 c7 40 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x240],ymm2
14ff: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
1503: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
1507: c5 fd 11 8c c7 60 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x260],ymm1
1510: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
1514: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
1518: c5 fd 11 ac c7 80 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x280],ymm5
1521: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
1525: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
1529: c5 fd 11 9c c7 a0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2a0],ymm3
1532: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
1536: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
153a: c5 fd 11 94 c7 c0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2c0],ymm2
1543: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
1547: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
154b: c5 fd 11 8c c7 e0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2e0],ymm1
1554: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
1558: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
155c: 48 83 c0 30 add rax,0x30
1560: 48 3d c1 0b 00 00 cmp rax,0xbc1
1566: 0f 82 24 ff ff ff jb 1490 <void compute2<3072>(double*)+0x50>
156c: c5 f8 77 vzeroupper
156f: c3 ret
有关的：
对于必须按顺序发生的操作，处理器的延迟界限和吞吐量界限 - 分析具有两个依赖链的代码，一个从另一个读取，一个在其自身之前读取。与强度降低循环相同的依赖模式，除了它的一个链是 FP 乘法。 （这也是一个多项式评估方案，但针对一个大多项式。）
从二阶导数计算出的曲线的 SIMD 优化是另一种能够跨过串行依赖的情况。
是否可以在计算中对串行依赖项使用 SIMD，例如指数移动平均滤波器？ - 如果前面有 n 步的封闭式公式，您可以使用它来回避串行依赖关系。
乱序执行，如何解决真正的依赖？ - 当一条指令依赖于一条尚未执行的指令时，CPU 必须等待。
依赖链分析链分析，来自 Agner Fog 的一个例子。
现代微处理器 90 分钟指南！ - 乱序执行和管道的一般背景。现代 CPU 风格的短向量 SIMD 以这种形式存在，通过单个 CPU 的流水线获得更多工作，而无需加宽流水线。相比之下，GPU 有许多简单的管道。
为什么 mulss 在 Haswell 上只需要 3 个周期，与 Agner 的指令表不同？ （展开具有多个累加器的 FP 循环）- 一些实验性的数字与展开以隐藏 FP 依赖链的延迟，以及寄存器重命名的一些 CPU 架构背景。
vaddpd），所以我在没有 AVX-512（我的桌面没有）的情况下每个元素获得 0.251 个周期。在测试时，我注意到您使用的是 rdtsc 数字而不是核心时钟周期，这是一个很大的假设。对于某些 Xeon，实际核心时钟在 running "heavy" 512-bit instructions 时接近 TSC 频率可能是正确的，但这是一个冒险的假设。
如果您需要此代码快速运行，或者您对此感到好奇，可以尝试以下方法：
您将 a[i] = f(i) 的计算更改为两个加法。修改它以使用两次加法计算 a[4i] = f(4i)，使用两次加法计算 a[4i+1] = f(4i+1)，依此类推。现在您有四个可以并行完成的计算。
编译器很有可能会执行相同的循环展开和矢量化，并且您具有相同的延迟，但是对于四个操作，而不是一个。
通过仅使用加法作为优化，您将丢失所有（较新 CPU 的）乘法管道的 GFLOPS，并且循环携带的依赖性通过停止自动矢量化使其变得更糟。如果它是自动向量化的，它将比乘法+加法快得多。并且每个数据的能源效率更高（仅添加优于 mul + add）。
另一个问题是数组末尾会收到更多舍入误差，因为累加的数量。但它不应该在非常大的数组之前可见（除非数据类型变为浮点数）。
当您使用 GCC 构建选项（在较新的 CPU 上）-std=c++20 -O3 -march=native -mavx2 -mprefer-vector-width=256 -ftree-vectorize -fno-math-errno 应用 Horner's scheme 时，
void f(double * const __restrict__ data){
double A=1.1,B=2.2,C=3.3;
for(int i=0; i<1024; i++) {
double id = double(i);
double result = A;
result *=id;
result +=B;
result *=id;
result += C;
data[i] = result;
}
}
编译器产生这个：
.L2:
vmovdqa ymm0, ymm2
vcvtdq2pd ymm1, xmm0
vextracti128 xmm0, ymm0, 0x1
vmovapd ymm7, ymm1
vcvtdq2pd ymm0, xmm0
vmovapd ymm6, ymm0
vfmadd132pd ymm7, ymm4, ymm5
vfmadd132pd ymm6, ymm4, ymm5
add rdi, 64
vpaddd ymm2, ymm2, ymm8
vfmadd132pd ymm1, ymm3, ymm7
vfmadd132pd ymm0, ymm3, ymm6
vmovupd YMMWORD PTR [rdi-64], ymm1
vmovupd YMMWORD PTR [rdi-32], ymm0
cmp rax, rdi
jne .L2
vzeroupper
ret
并使用 -mavx512f -mprefer-vector-width=512：
.L2:
vmovdqa32 zmm0, zmm3
vcvtdq2pd zmm4, ymm0
vextracti32x8 ymm0, zmm0, 0x1
vcvtdq2pd zmm0, ymm0
vmovapd zmm2, zmm4
vmovapd zmm1, zmm0
vfmadd132pd zmm2, zmm6, zmm7
vfmadd132pd zmm1, zmm6, zmm7
sub rdi, -128
vpaddd zmm3, zmm3, zmm8
vfmadd132pd zmm2, zmm5, zmm4
vfmadd132pd zmm0, zmm5, zmm1
vmovupd ZMMWORD PTR [rdi-128], zmm2
vmovupd ZMMWORD PTR [rdi-64], zmm0
cmp rax, rdi
jne .L2
vzeroupper
ret
由于 mul+add 加入到单个 FMA 中，所有浮点运算都是“打包”向量形式和更少的指令（它是一个两次展开的版本）。每 64 字节数据有 16 条指令（如果 AVX-512 则为 128 字节）。
Horner 方案的另一个好处是它在 FMA 指令中的计算精度更高，并且每次循环迭代只有 O(1) 次操作，因此它不会在更长的数组中累积那么多错误。
我认为 Agner Fog 的优化手册中的优化必须来自 Quake III fast inverse square root approximation 时代，这仅适用于 x87，但obsoleted by SSE。在那些时候，SIMD 的宽度不足以产生太大的差异，而且对于 sqrt() 和除法也很慢，但有 rsqrtss。手册说版权所有 2004，因此有 Celeron 个 CPU 带有 SSE 而不是 FMA。第一个 AVX 台式机 CPU 推出的时间要晚得多，而 FMA 甚至更晚。
这是另一个降低强度的版本（用于 id 值）：
void f(double * const __restrict__ data){
double B[]={2.2,2.2,2.2,2.2,2.2,2.2,2.2,2.2,
2.2,2.2,2.2,2.2,2.2,2.2,2.2,2.2};
double C[]={3.3,3.3,3.3,3.3,3.3,3.3,3.3,3.3,
3.3,3.3,3.3,3.3,3.3,3.3,3.3,3.3};
double id[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
for(long long i=0; i<1024; i+=16) {
double result[]={1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,
1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1};
// The same thing, just with explicit auto-vectorization help
for(int j=0;j<16;j++)
{
result[j] *=id[j];
result[j] +=B[j];
result[j] *=id[j];
result[j] += C[j];
data[i+j] = result[j];
}
// strength reduction
for(int j=0;j<16;j++)
{
id[j] += 16.0;
}
}
}
集会：
.L2:
vmovapd zmm3, zmm0
vmovapd zmm2, zmm1
sub rax, -128
vfmadd132pd zmm3, zmm6, zmm7
vfmadd132pd zmm2, zmm6, zmm7
vfmadd132pd zmm3, zmm5, zmm0
vfmadd132pd zmm2, zmm5, zmm1
vaddpd zmm0, zmm0, zmm4
vaddpd zmm1, zmm1, zmm4
vmovupd ZMMWORD PTR [rax-128], zmm3
vmovupd ZMMWORD PTR [rax-64], zmm2
cmp rdx, rax
jne .L2
vzeroupper
ret
当数据、A、B 和 C 数组按 alignas(64) 对齐并且数据数组大小足够小时，它会以 0.26 cycles per element 运行速度。
mulpd 在与 addpd 不同的端口上运行的 CPU 上，因此它们只竞争前端吞吐量。如果您只关心速度而不是准确性，那么在 gnasher 的回答中建议的策略（我在前面的评论中建议过）使用多个累加器来隐藏 FP 延迟，可能会更好，避免 int->FP 成本。
int->FP 成本，并且无法通过积极展开来隐藏。用 std::memcpy 而不是强制转换可能更好地代表一些魔力。当我有更多时间时，我会测试它。 （如果循环计数小于 53 位，它应该可以工作）
3 ops * 0.5c/op / 8 elements per ZMM vector = 每个元素 0.1875 个周期。但是还有两个（消除的）vmovapd 指令和两个存储，这样就填满了 Skylake-X 上的整个 4 宽管道；只有 Ice Lake 更宽的管道也可以运行循环开销。但是 Ice Lake 禁用了 mov 消除（至少对于整数，我忘记了 SIMD）所以那些 vmovapd 指令将与 FMA 竞争。
vaddpd 延迟（SKX 或 ICX，使用 Alder Lake 时只有 3 个周期）。需要更多的展开来隐藏任何当前 CPU 上的延迟；你在这里只展开 2 个 ZMM 向量。 （当然，输出阵列应该适合 L1d 缓存，因为您需要每 1.5 个时钟周期存储一次，每 3 个 FP 数学运算一个结果向量 = 每 1.5 个周期一个）4 个周期延迟，所需吞吐量为每 1.5 个时钟周期（对于 vaddpd）需要展开至少 4/1.5 = 2.666。所以也可以做4。
vaddpd dep 链。也许只是来自其他逻辑核心、嘈杂的基准测试环境的竞争。）
与加法相比，乘法在我们的 CPU 中速度明显较慢而闻名。
这在历史上可能是正确的，并且对于更简单的低功耗 CPU 可能仍然正确，但是如果 CPU 设计人员准备“解决问题”，则乘法几乎可以与加法一样快。
现代 CPU 旨在通过流水线和具有多个执行单元的组合同时处理多条指令。
但是，问题在于数据依赖性。如果一条指令依赖于另一条指令的结果，那么在它所依赖的指令完成之前，它的执行不能开始。
现代 CPU 试图通过“乱序执行”来解决这个问题。等待数据的指令可以保持排队，同时允许执行其他指令。
但即使采取了这些措施，有时 CPU 也会简单地用完新的工作来调度。
通过手动将代码并行化为以下内容，您似乎也可以吃到蛋糕：
double A4 = A+A+A+A;
double Z = 3A+B;
double Y1 = C;
double Y2 = A+B+C;
int i;
// ... setup unroll when LEN is odd...
for(i=0; i<LEN; i++) {
data[i] = Y1;
data[++i] = Y2;
Y1 += Z;
Y2 += Z;
Z += A4;
}
它可能不像写的那样完全起作用，但你明白了：展开循环，以便每个依赖数据的路径都可以并行完成。对于正在考虑的机器，四步展开应该可以实现最大性能，但是当然，您会获得在软件中对架构进行硬编码所带来的所有有趣的事情。
addpd，因此 6 标量添加在飞行中）；使用 AVX 在 Sandy Bridge 上的人数增加了一倍。
关注公众号
不定期副业成功案例分享
i++是一个循环携带的 dep，但编译器可以使用它，因为整数数学是关联的，不像没有-ffast-math的 FP）double，而 SIMD FPaddsd/addpd在 Nehalem 和 Sandy Bridge 上具有 3 个周期的延迟和 1 个周期的吞吐量。 （虽然在循环中有两个单独的添加链，这可能是每 1.5 个时钟周期添加一个标量 FP，所以是的，也许 SIMD 更重要。）-ffast-math），您希望它以一种对展开友好的方式进行，以允许矢量化。