为什么这段代码在强度降低乘法到循环携带加法之后执行得更慢？
Why does this code execute more slowly after strength-reducing multiplications to loop-carried additions?
我正在阅读 Agner Fog 的优化手册，我遇到了这个例子：
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 表示有一种方法可以优化此代码 - 通过意识到循环可以避免使用代价高昂的乘法，而是使用每次迭代应用的“增量”。
我用一张纸来证实理论，首先...
...当然，他是对的 - 在每次循环迭代中，我们可以通过添加“增量”来根据旧结果计算新结果。此增量从值“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 中以慢得多而著称。我们用 2 次加法代替了 3 次乘法和 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.
嗯，这出乎意料。由于我们报告的是最短执行时间，因此我们正在丢弃由 OS 的各个部分引起的“噪音”。我们还注意 运行 在一台什么都不做的机器上。结果或多或少是可重复的——重新运行宁这两个二进制文件表明这是一个一致的结果：
# 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 的第一个版本——“愚蠢”但不知何故速度很快的代码——有一个看起来像这样的循环：
第二个优化版本怎么样 - 只做了两次加法？
现在我不了解你，但就我自己而言，我很……困惑。第二个版本的指令少了大约 4 倍，其中两个主要指令只是基于 SSE 的添加 (addsd)。第一个版本，不仅有 4 倍多的指令......它也充满了（正如预期的那样）乘法 (mulpd)。
我承认我没想到这个结果。不是因为我是 Agner 的粉丝（我是，但这无关紧要）。
知道我错过了什么吗？我在这里犯了什么错误，可以解释速度上的差异吗？请注意，我已经在 Xeon W5580 和 Xeon E5 1620 上进行了测试 - 在这两个版本中，第一个（哑）版本比第二个快得多。
编辑：为了便于重现结果，我在两个版本的代码中添加了两个要点：Dumb yet somehow faster and optimized, yet somehow slower.
P.S。请不要评论浮点精度问题；这不是这个问题的重点。
理解您所看到的性能差异的关键在于矢量化。是的，addition-based 解决方案在其内部循环中只有两条指令，但重要的区别不在于 循环中有多少条 指令，而是在 每条指令执行了多少工作。
在第一个版本中，输出完全依赖于输入：每个 data[i] 都是 i 本身的函数，这意味着每个 data[i] 都可以在任何顺序：编译器可以向前、向后、向侧面执行它们，无论如何，您仍然会得到相同的结果 — 除非您从另一个线程观察该内存，否则您永远不会注意到数据被处理的方式。
在第二个版本中，输出不依赖于 i — 它依赖于上次循环的 A 和 Z。
如果我们将这些循环的主体表示为小的数学函数，它们将具有非常不同的整体形式：
- f(i) -> di
- f(Y, Z) -> (di, Y', Z')
在后一种形式中，没有对 i 的实际依赖 — 计算函数值的唯一方法是通过了解先前的 Y 和 Z函数的最后一次调用，这意味着函数形成一个链 - 在完成上一个函数之前不能执行下一个。
为什么这很重要？因为CPU有vector parallel指令，each可以同时进行2个，4个，甚至8个算术运算！ （AVX CPUs 可以并行执行更多操作。）这是四次乘法、四次加法、四次减法、四次比较——四次随便！因此，如果您尝试计算的输出 仅 取决于输入，那么您可以安全地一次计算两个、四个甚至八个——它们是否无关紧要'向前或向后，因为结果是一样的。但是，如果输出依赖于 先前的计算 ，那么您将不得不以串行形式进行操作 — 一次一个。
这就是为什么“更长”的代码在性能方面胜出的原因。尽管它有更多的设置，而且它实际上 做 更多的工作，但大部分工作是并行完成的：它不是在每次迭代中只计算 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 有特殊的built-in 指令可以同时对不同数据执行多个 * 或多个 + 操作，每个操作只需要一个时钟周期。
注意较快解决方案中说明中的字母 p — addpd 和 mulpd — 以及较慢解决方案中说明中的字母 s — addsd。那是“添加打包双打”和“多打包双打”与“添加单双打”。
不仅如此，看起来编译器也部分展开了循环 — 循环不仅在每次迭代中执行 两个 个值，而且实际上 四个，交错操作以避免依赖和停顿，所有这些都减少了汇编代码必须测试的次数 i < 1000。
不过，所有这些仅在循环迭代之间没有依赖关系的情况下才有效：如果唯一确定每个 data[i] 会发生什么的事情是i 本身。如果存在依赖关系，如果上一次迭代的数据影响下一次迭代，那么编译器可能会受到它们的限制，以至于它根本无法更改代码——而不是编译器能够使用花哨的并行指令或巧妙的优化（CSE、强度降低、循环展开、重新排序等），您得到的代码正是您输入的代码 — 添加 Y，然后添加 Z，然后重复。
但是在这里，在第一个版本的代码中，编译器正确地识别出数据中没有依赖关系，并发现它可以并行完成工作，所以它做到了，这就是完全不同。
这里的主要区别在于循环依赖性。第二种情况下的循环是 dependent——循环中的操作依赖于前一次迭代。这意味着每次迭代甚至不能在上一次迭代完成之前开始。在第一种情况下，循环体是完全 独立的 -- 循环体中的所有内容都是自包含的，仅取决于迭代计数器和常量值。这意味着循环可以并行计算——多个迭代可以同时进行。然后，这允许循环被平凡地展开和矢量化，重叠许多指令。
如果您查看性能计数器（例如使用 perf stat ./1），您会发现第一个循环除了 运行 更快之外还 运行 每个指令更多周期（IPC）。相比之下，第二个循环有更多的依赖周期——CPU 无所事事，等待指令完成，然后才能发出更多指令的时间。
第一个可能会成为内存带宽瓶颈，特别是如果您让编译器 auto-vectorize 在您的 Sandybridge (gcc -O3 -march=native) 上使用 AVX，如果它设法使用 256 位向量。那时 IPC 将下降，特别是对于对于 L3 缓存来说太大的输出数组。
请注意，展开和矢量化不需要 独立循环——当存在（某些）循环依赖项时，您可以执行它们。然而，它更难 并且 回报更少。因此，如果您想看到矢量化的最大加速，它有助于尽可能消除循环依赖性。
如果您需要此代码 运行 快速，或者如果您很好奇，可以尝试以下操作：
您将 a[i] = f(i) 的计算更改为两次加法。修改它以使用两次加法计算 a[4i] = f(4i)，使用两次加法计算 a[4i+1] = f(4i+1)，等等。现在您有四个可以并行完成的计算。
编译器很有可能会执行相同的循环展开和矢量化，并且您有相同的延迟，但对于四个操作，而不是一个。
通过仅使用加法作为优化，您将错过（较新的 CPU）乘法管道的所有 gflops，并且循环携带的依赖性通过停止 auto-vectorization 使情况变得更糟。如果它被自动矢量化，它会比乘法+加法快得多。每个数据的能源效率更高（仅添加比 mul+add 更好）。
另一个问题是，由于加法次数的累积，数组的末尾会收到更多的舍入误差。但在非常大的数组之前它不应该是可见的（除非数据类型变为 float）。
当您应用带有 GCC 构建选项的 Horner Scheme 时（在较新的 CPU 上）-std=c++20 -O3 -march=native -mavx2 -mprefer-vector-width=256 -ftree-vectorize -fno-math-errno，
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
所有 FP 操作都是“压缩”向量形式和更少的指令（它是 twice-unrolled 版本），因为 mul+add 加入单个 FMA。每 64 字节数据 16 条指令（如果是 AVX512，则为 128 字节）。
Horner Scheme 的另一个好处是它在 FMA 指令中的计算精度更高，并且每次循环迭代只需 O(1) 次操作，因此它不会在较长的数组中累积那么多错误。
我认为Agner Fog的优化手册中的优化必须来自Quake-3快速逆平方根近似时代。在那个时候，SIMD 不够广泛，无法产生太大的差异，而且缺乏对 sqrt 函数的支持。手册上写着版权 2004，所以赛扬有 SSE 而不是 FMA。第一个 AVX 桌面 cpu 推出得晚得多，FMA 甚至更晚。
这是另一个带有 strength-reduction 的版本（对于 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};
// 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 周期运行 每个元素 速度。
Multiplications have a reputation for being significantly slower in our CPUs, compared to additions.
这在历史上可能是正确的，对于更简单的低功耗 CPUs 可能仍然是正确的，但如果 CPU 设计人员准备好“解决问题”，则可以进行乘法运算几乎和加法一样快。
现代 CPU 旨在通过流水线和多个执行单元的组合同时处理多条指令。
这个问题虽然是数据依赖性。如果一条指令依赖于另一条指令的结果，那么在它所依赖的指令完成之前，它的执行不能开始。
现代 CPUs 试图通过“无序执行”来解决这个问题。等待数据的指令可以保持排队，同时允许执行其他指令。
但即使采取这些措施，有时 CPU 也可以简单地 运行 安排新工作。
这个 method of finite differences strength-reduction 优化 可以 提供比你能做的最好的加速 re-evaluating 分别为每个 i.但只有当你将它推广到更大的步幅时，循环中仍然有足够的并行性。 我的版本在我的 Skylake 上每个时钟周期存储 1 个矢量（4 个双精度），用于适合 L1d 缓存的小型阵列，否则就是带宽测试。在较早的英特尔上，它还应该最大化 SIMD FP-add 吞吐量，包括带有 AVX 的 Sandybridge（1x 256 位 add/clock 和每 2 个时钟 1x 256 位存储）。
对上一次迭代的值的依赖是致命的
这个strength-reduction optimization（只是添加而不是从一个新的i开始并乘法）引入了跨循环迭代的串行依赖性，涉及FP数学而不是比整数增量。
原件在每个输出元素上都有数据并行性：每个元素只依赖于常量和它自己的i值。编译器可以 auto-vectorize 与 SIMD（SSE2，或 AVX，如果你使用 -O3 -march=native），并且 CPUs 可以将跨循环迭代的工作与 out-of-order 执行重叠。尽管有大量的额外工作，CPU 还是能够在编译器的帮助下应用足够的蛮力。
但是根据poly(i)计算poly(i+1)的版本并行度非常有限；没有 SIMD 矢量化，并且您的 CPU 每 4 个周期只能 运行 两次标量加法，例如，其中 4 个周期是通过 Tiger Lake 在 Intel Skylake 上进行 FP 加法的延迟。 (https://uops.info/).
@huseyin tugrul buyukisik 的回答显示了如何在更现代的 CPU 上接近最大化原始版本的吞吐量，其中有两个 FMA 运算来评估多项式（Horner 的方案），加上一个int->FP 转换或 FP 增量。 （后者创建了一个 FP dep 链，您需要展开才能隐藏它。）
所以最好的情况是每个 SIMD 输出向量有 3 个 FP 数学运算。 （加上商店）。目前的 Intel CPUs 只有两个 FP 执行单元可以 运行 FP 数学运算，包括 int->double。 （对于 512 位向量，当前 CPUs 关闭了端口 1 上的向量 ALU，因此只有 2 个 SIMD ALU 端口，因此 non-FP-math ops like SIMD-integer increment 也会竞争用于 SIMD 吞吐量。除了 CPUs 只有一个 512 位 FMA 单元外，端口 5 可用于其他工作。）
AMD因为Zen2在两个端口上有两个FMA/mul单元，在两个不同的端口上有两个FP add/sub单元，所以如果你用FMA做加法，理论上最多有四个每个时钟周期的 SIMD 添加。
Haswell/Broadwell 有 2/clock FMA，但只有 1/clock FP add/sub（延迟较低）。这对原始代码很有用，not great 对于经过优化以具有大量并行性的代码。这可能就是英特尔在 Skylake 中更改它的原因。
您的 Sandybridge (E5-1620) 和 Nehalem (W5580) CPUs 在不同的端口上有 1/clock FP add/sub、1/clock FP mul。这就是 Haswell 的基础。为什么添加额外的乘法不是大问题：它们可以 运行 与现有的加法并行。 （Sandybridge 的宽度为 256 位，但您在未启用 AVX 的情况下进行编译：使用 -march=native。）
寻找平行度：strength-reduction 任意步幅
您的 compute2 根据前一个值计算下一个 Y 和下一个 Z。即，步幅为 1，即 data[i+1] 所需的值。所以每次迭代都依赖于前一次迭代。
如果将其推广到其他步幅，则可以推进 4、6、8 或更多独立的 Y 和 Z 值，以便它们彼此同步跳跃，彼此独立。 这为编译器 and/or 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 的差异。non-power-of-2 有助于显示步幅的 factors/powers 与 naturally-occurring 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]，然后我们可以根据需要对其进行代数化简。 （对于 compile-time 常量 A、B、C，编译器将 constant-propagate 任一种方式。）
Z[j] 在跨越 poly(x) 时保持 first-order 的差异，起点为 poly(0..stride-1)。这是上面table.
中的中间一栏
Z[j] += second_difference 的必要更新是一个标量常数，正如我们从 second-order 差异相同可以看出的那样。
通过研究几个不同的 stride 和 A 值（i^2 的系数），我们可以看到它是 A * 2 * (stride * stride). （使用像 3 和 5 这样的 non-coprime 值有助于理清事物。）有了更多的代数知识，您就可以象征性地展示这一点。 2 的因子从微积分 PoV 来看是有意义的：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].[= 交互113=]
有 stride 个独立的并行 dep 链供编译器 (SIMD) 和 CPU（流水线执行单元）利用，运行宁 stride比原始速度快 1 倍，直到您在 SIMD FP-add 吞吐量而不是延迟上出现瓶颈，或者如果您的缓冲区不适合 L1d，则存储带宽。 （或者直到编译器面对并且没有很好地/根本没有展开和矢量化这些循环的地步！）
这在实践中是如何编译的：很好地使用 clang
(Godbolt compiler explorer) Clang auto-vectorize 与 stride=16 （4x YMM 向量，每个向量 4 doubles）与clang14 -O3 -march=skylake -ffast-math.
看起来 clang 进一步展开了 2，将 Z[j] += diff2 缩短为 tmp = Z[j] + diff2; / Z[j] += 2*diff2;。这减轻了 Z dep 链上的压力，只留下 Y[j] 来应对 Skylake 上的延迟瓶颈。
因此每个 asm 循环迭代执行 2x 8 vaddpd 指令和 2x 4 存储。循环开销是 add + macro-fused cmp/jne，所以 2 微指令。 （或者对于全局数组，只有一个 add/jne uop，将负索引向上计数为零；它相对于数组的末尾进行索引。）
Skylake 运行 在几乎 1 个存储和每个时钟周期 2x vaddpd 中实现。这是这两件事的最大吞吐量。 front-end 只需要跟上 3 微码/时钟周期多一点，但自 Core2 以来它一直是 4 宽。
Sandybridge-family 中的 uop 缓存使这没问题。 （除非你 运行 进入 Skylake 上的 JCC 勘误表，所以我使用 -mbranches-within-32B-boundaries to have clang pad instructions to avoid that。）
Skylake 的 vaddpd 延迟为 4 个周期，来自 stride=16 的 4 个 dep 链仅勉强足以保持 4 个独立操作在运行中。任何时候 Y[j]+= 没有 运行 循环就绪，就会产生气泡。由于 clang 对 Z[] 链的额外展开，Z[j]+= 可以提前 运行，因此 Z 链可以领先。使用 oldest-ready-first 调度，它往往会稳定到 Yj+= uops 没有冲突的状态，显然，因为它确实 运行 在我的 Skylake 上全速运行。如果我们能让编译器仍然为 stride=32 生成漂亮的 asm，那将留下更多空间，但不幸的是它没有。 （以对奇数尺寸进行更多清理工作为代价。）
奇怪的是，Clang 仅使用 -ffast-math 对其进行矢量化。下面的完整基准测试中的模板版本不需要 --fast-math。源代码被仔细编写为 SIMD-friendly，并按源代码顺序进行数学运算。 （不过，Fast-math 允许 clang 进一步展开 Z 增量。）
编写循环的另一种方法是使用一个内部循环而不是所有 Y 操作，然后是所有 Z 操作。这在下面的基准测试中很好（有时实际上更好），但这里即使使用 -ffast-math 也不会矢量化。对于像这样的 non-trivial 问题，从编译器中获得最佳展开的 SIMD asm 可能是繁琐且不可靠的，并且可能需要一些时间。
我将它包含在 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] 进行手动矢量化。但是这个版本可以为 AArch64 等其他 ISA 向量化。
当 problem-size 不是展开的倍数时，大展开因子的其他缺点是留下更多的清理工作。 （您可以使用 compute1 策略进行清理，让编译器在执行标量之前将其矢量化一两次迭代。）
理论上编译器是 allowe 使用 -ffast-math 为您执行此操作，或者从 compute1 对原始多项式执行 strength-reduction，或者从 compute2 查看步幅如何累积。
但在实践中，这真的很复杂，需要人类自己来做。除非/直到有人开始教编译器如何寻找这样的模式并应用差异方法本身，并选择大步前进！但是，即使使用 -ffast-math，对具有不同 error-accumulation 属性的算法进行大规模重写也可能是不可取的。 （整数不会有任何精度问题，但它仍然是一个复杂的 pattern-match / 替换。）
实验性能结果：
我用 clang13.0.0 在我的台式机 (i7-6700k) 上进行了测试。这实际上是 运行 每个时钟周期 1 个 SIMD 存储，具有多种编译器选项组合（fast-math 或不）和 #if 0 与 #if 1 在内循环策略上。我的基准/测试框架基于@huseyin tugrul buyukisik 的版本，改进后在 rdtsc 指令之间重复更可测量的数量，并使用测试循环来检查多项式简单计算的正确性。
我还让它补偿了核心时钟频率与 "reference" frequency of the TSC read by rdtsc 之间的差异，在我的例子中是 3.9GHz 与 4008 MHz。 （额定最大睿频频率为 4.2GHz，但在 Linux 上 EPP = balance_performance，它只想达到 3.9 GHz。）
源代码 on Godbolt：使用一个内部循环，而不是 3 个单独的 j<16 循环，并且 不 使用 -ffast-math。使用 __attribute__((noinline)) 来防止它内联到重复循环中。选项和源的一些其他变体导致循环内的一些 vpermpd 洗牌。
下面的基准数据来自以前的版本，Z[j] 初始值设定项有问题，但循环 asm 相同。 Godbolt link 现在有一个正确性测试在定时循环之后，它通过了。实际性能在我的桌面上仍然相同，每 double 仅超过 0.25 个周期，即使没有 #if 1 / -ffast-math 允许 clang 额外展开。
$ 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.9GHz)
0.252109 cycles per data element (corrected from ref cycles to core clocks for i7-6700k @ 3.9GHz)
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/时钟，显然没有受到 sub-optimal uop 调度的影响。 （我增加了重复循环，所以程序的大部分时间都花在那里，使整个程序的 perf stat 很有用。）
此优化的目标是以更少的 FLOPS 完成相同的工作，但这也意味着我们实质上是在不使用 FMA 的情况下最大化 Skylake 的 8 FLOP/clock 限制。 （单核 3.9GHz 时 30.58 GFLOP/s）。
Asm 的non-inline 函数(objdump -drwC -Mintel); clang 使用 4 对 Y，Z 对 YMM 向量，并将循环进一步展开 3 倍，使其成为 24KiB 大小的精确倍数，无需清理。注意 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
相关：
- 分析具有两个 dep 链的代码，一个从另一个读取，一个更早。与 strength-reduced 循环相同的依赖模式，除了它的一个链是一个 FP 乘法。 （这也是一种多项式评估方案，但是对于一个大多项式。） 又一个跨越串行依赖的案例 - 如果前面有 n 步的 closed-form 公式，您可以使用它来回避串行依赖性。 - Out of Order Execution, How to Solve True Dependency? - CPU当指令依赖于尚未执行的指令时必须等待。
non-loop-carried 依赖链分析，来自 Agner Fog 的一个例子。 - 现代微处理器
90 分钟指南！ - out-of-order 执行和管道的一般背景。现代 CPU-style short-vector SIMD 以这种形式存在，以便在不加宽管道的情况下通过单个 CPU 的管道获得更多工作。相比之下，GPU 有很多简单的流水线。
- 一些实验数字展开以隐藏 FP 依赖链的延迟，以及一些 CPU-architecture 寄存器重命名背景。
看起来你可以吃蛋糕也可以吃，通过手动将代码并行化为这样的东西：
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;
}
可能不完全像写的那样起作用，但你明白了：展开循环，这样数据相关的路径就可以并行完成。对于正在考虑的机器，4 步展开应该可以实现最佳性能，但是当然，您会得到 hard-coding 软件架构带来的所有有趣的东西。
我正在阅读 Agner Fog 的优化手册，我遇到了这个例子：
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 表示有一种方法可以优化此代码 - 通过意识到循环可以避免使用代价高昂的乘法，而是使用每次迭代应用的“增量”。
我用一张纸来证实理论，首先...
...当然，他是对的 - 在每次循环迭代中，我们可以通过添加“增量”来根据旧结果计算新结果。此增量从值“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 中以慢得多而著称。我们用 2 次加法代替了 3 次乘法和 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.
嗯，这出乎意料。由于我们报告的是最短执行时间，因此我们正在丢弃由 OS 的各个部分引起的“噪音”。我们还注意 运行 在一台什么都不做的机器上。结果或多或少是可重复的——重新运行宁这两个二进制文件表明这是一个一致的结果：
# 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 的第一个版本——“愚蠢”但不知何故速度很快的代码——有一个看起来像这样的循环：
第二个优化版本怎么样 - 只做了两次加法？
现在我不了解你，但就我自己而言，我很……困惑。第二个版本的指令少了大约 4 倍，其中两个主要指令只是基于 SSE 的添加 (addsd)。第一个版本，不仅有 4 倍多的指令......它也充满了（正如预期的那样）乘法 (mulpd)。
我承认我没想到这个结果。不是因为我是 Agner 的粉丝（我是，但这无关紧要）。
知道我错过了什么吗？我在这里犯了什么错误，可以解释速度上的差异吗？请注意，我已经在 Xeon W5580 和 Xeon E5 1620 上进行了测试 - 在这两个版本中，第一个（哑）版本比第二个快得多。
编辑：为了便于重现结果，我在两个版本的代码中添加了两个要点：Dumb yet somehow faster and optimized, yet somehow slower.
P.S。请不要评论浮点精度问题；这不是这个问题的重点。
理解您所看到的性能差异的关键在于矢量化。是的，addition-based 解决方案在其内部循环中只有两条指令，但重要的区别不在于 循环中有多少条 指令，而是在 每条指令执行了多少工作。
在第一个版本中，输出完全依赖于输入：每个 data[i] 都是 i 本身的函数，这意味着每个 data[i] 都可以在任何顺序：编译器可以向前、向后、向侧面执行它们，无论如何，您仍然会得到相同的结果 — 除非您从另一个线程观察该内存，否则您永远不会注意到数据被处理的方式。
在第二个版本中，输出不依赖于 i — 它依赖于上次循环的 A 和 Z。
如果我们将这些循环的主体表示为小的数学函数，它们将具有非常不同的整体形式：
- f(i) -> di
- f(Y, Z) -> (di, Y', Z')
在后一种形式中，没有对 i 的实际依赖 — 计算函数值的唯一方法是通过了解先前的 Y 和 Z函数的最后一次调用，这意味着函数形成一个链 - 在完成上一个函数之前不能执行下一个。
为什么这很重要？因为CPU有vector parallel指令，each可以同时进行2个，4个，甚至8个算术运算！ （AVX CPUs 可以并行执行更多操作。）这是四次乘法、四次加法、四次减法、四次比较——四次随便！因此，如果您尝试计算的输出 仅 取决于输入，那么您可以安全地一次计算两个、四个甚至八个——它们是否无关紧要'向前或向后，因为结果是一样的。但是，如果输出依赖于 先前的计算 ，那么您将不得不以串行形式进行操作 — 一次一个。
这就是为什么“更长”的代码在性能方面胜出的原因。尽管它有更多的设置，而且它实际上 做 更多的工作，但大部分工作是并行完成的：它不是在每次迭代中只计算 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 有特殊的built-in 指令可以同时对不同数据执行多个 * 或多个 + 操作，每个操作只需要一个时钟周期。
注意较快解决方案中说明中的字母 p — addpd 和 mulpd — 以及较慢解决方案中说明中的字母 s — addsd。那是“添加打包双打”和“多打包双打”与“添加单双打”。
不仅如此，看起来编译器也部分展开了循环 — 循环不仅在每次迭代中执行 两个 个值，而且实际上 四个，交错操作以避免依赖和停顿，所有这些都减少了汇编代码必须测试的次数 i < 1000。
不过，所有这些仅在循环迭代之间没有依赖关系的情况下才有效：如果唯一确定每个 data[i] 会发生什么的事情是i 本身。如果存在依赖关系，如果上一次迭代的数据影响下一次迭代，那么编译器可能会受到它们的限制，以至于它根本无法更改代码——而不是编译器能够使用花哨的并行指令或巧妙的优化（CSE、强度降低、循环展开、重新排序等），您得到的代码正是您输入的代码 — 添加 Y，然后添加 Z，然后重复。
但是在这里，在第一个版本的代码中，编译器正确地识别出数据中没有依赖关系，并发现它可以并行完成工作，所以它做到了，这就是完全不同。
这里的主要区别在于循环依赖性。第二种情况下的循环是 dependent——循环中的操作依赖于前一次迭代。这意味着每次迭代甚至不能在上一次迭代完成之前开始。在第一种情况下，循环体是完全 独立的 -- 循环体中的所有内容都是自包含的，仅取决于迭代计数器和常量值。这意味着循环可以并行计算——多个迭代可以同时进行。然后，这允许循环被平凡地展开和矢量化，重叠许多指令。
如果您查看性能计数器（例如使用 perf stat ./1），您会发现第一个循环除了 运行 更快之外还 运行 每个指令更多周期（IPC）。相比之下，第二个循环有更多的依赖周期——CPU 无所事事，等待指令完成，然后才能发出更多指令的时间。
第一个可能会成为内存带宽瓶颈，特别是如果您让编译器 auto-vectorize 在您的 Sandybridge (gcc -O3 -march=native) 上使用 AVX，如果它设法使用 256 位向量。那时 IPC 将下降，特别是对于对于 L3 缓存来说太大的输出数组。
请注意，展开和矢量化不需要 独立循环——当存在（某些）循环依赖项时，您可以执行它们。然而，它更难 并且 回报更少。因此，如果您想看到矢量化的最大加速，它有助于尽可能消除循环依赖性。
如果您需要此代码 运行 快速，或者如果您很好奇，可以尝试以下操作：
您将 a[i] = f(i) 的计算更改为两次加法。修改它以使用两次加法计算 a[4i] = f(4i)，使用两次加法计算 a[4i+1] = f(4i+1)，等等。现在您有四个可以并行完成的计算。
编译器很有可能会执行相同的循环展开和矢量化，并且您有相同的延迟，但对于四个操作，而不是一个。
通过仅使用加法作为优化，您将错过（较新的 CPU）乘法管道的所有 gflops，并且循环携带的依赖性通过停止 auto-vectorization 使情况变得更糟。如果它被自动矢量化，它会比乘法+加法快得多。每个数据的能源效率更高（仅添加比 mul+add 更好）。
另一个问题是，由于加法次数的累积，数组的末尾会收到更多的舍入误差。但在非常大的数组之前它不应该是可见的（除非数据类型变为 float）。
当您应用带有 GCC 构建选项的 Horner Scheme 时（在较新的 CPU 上）-std=c++20 -O3 -march=native -mavx2 -mprefer-vector-width=256 -ftree-vectorize -fno-math-errno，
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
所有 FP 操作都是“压缩”向量形式和更少的指令（它是 twice-unrolled 版本），因为 mul+add 加入单个 FMA。每 64 字节数据 16 条指令（如果是 AVX512，则为 128 字节）。
Horner Scheme 的另一个好处是它在 FMA 指令中的计算精度更高，并且每次循环迭代只需 O(1) 次操作，因此它不会在较长的数组中累积那么多错误。
我认为Agner Fog的优化手册中的优化必须来自Quake-3快速逆平方根近似时代。在那个时候，SIMD 不够广泛，无法产生太大的差异，而且缺乏对 sqrt 函数的支持。手册上写着版权 2004，所以赛扬有 SSE 而不是 FMA。第一个 AVX 桌面 cpu 推出得晚得多，FMA 甚至更晚。
这是另一个带有 strength-reduction 的版本（对于 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};
// 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 周期运行 每个元素 速度。
Multiplications have a reputation for being significantly slower in our CPUs, compared to additions.
这在历史上可能是正确的，对于更简单的低功耗 CPUs 可能仍然是正确的，但如果 CPU 设计人员准备好“解决问题”，则可以进行乘法运算几乎和加法一样快。
现代 CPU 旨在通过流水线和多个执行单元的组合同时处理多条指令。
这个问题虽然是数据依赖性。如果一条指令依赖于另一条指令的结果，那么在它所依赖的指令完成之前，它的执行不能开始。
现代 CPUs 试图通过“无序执行”来解决这个问题。等待数据的指令可以保持排队，同时允许执行其他指令。
但即使采取这些措施，有时 CPU 也可以简单地 运行 安排新工作。
这个 method of finite differences strength-reduction 优化 可以 提供比你能做的最好的加速 re-evaluating 分别为每个 i.但只有当你将它推广到更大的步幅时，循环中仍然有足够的并行性。 我的版本在我的 Skylake 上每个时钟周期存储 1 个矢量（4 个双精度），用于适合 L1d 缓存的小型阵列，否则就是带宽测试。在较早的英特尔上，它还应该最大化 SIMD FP-add 吞吐量，包括带有 AVX 的 Sandybridge（1x 256 位 add/clock 和每 2 个时钟 1x 256 位存储）。
对上一次迭代的值的依赖是致命的
这个strength-reduction optimization（只是添加而不是从一个新的i开始并乘法）引入了跨循环迭代的串行依赖性，涉及FP数学而不是比整数增量。
原件在每个输出元素上都有数据并行性：每个元素只依赖于常量和它自己的i值。编译器可以 auto-vectorize 与 SIMD（SSE2，或 AVX，如果你使用 -O3 -march=native），并且 CPUs 可以将跨循环迭代的工作与 out-of-order 执行重叠。尽管有大量的额外工作，CPU 还是能够在编译器的帮助下应用足够的蛮力。
但是根据poly(i)计算poly(i+1)的版本并行度非常有限；没有 SIMD 矢量化，并且您的 CPU 每 4 个周期只能 运行 两次标量加法，例如，其中 4 个周期是通过 Tiger Lake 在 Intel Skylake 上进行 FP 加法的延迟。 (https://uops.info/).
@huseyin tugrul buyukisik 的回答显示了如何在更现代的 CPU 上接近最大化原始版本的吞吐量，其中有两个 FMA 运算来评估多项式（Horner 的方案），加上一个int->FP 转换或 FP 增量。 （后者创建了一个 FP dep 链，您需要展开才能隐藏它。）
所以最好的情况是每个 SIMD 输出向量有 3 个 FP 数学运算。 （加上商店）。目前的 Intel CPUs 只有两个 FP 执行单元可以 运行 FP 数学运算，包括 int->double。 （对于 512 位向量，当前 CPUs 关闭了端口 1 上的向量 ALU，因此只有 2 个 SIMD ALU 端口，因此 non-FP-math ops like SIMD-integer increment 也会竞争用于 SIMD 吞吐量。除了 CPUs 只有一个 512 位 FMA 单元外，端口 5 可用于其他工作。）
AMD因为Zen2在两个端口上有两个FMA/mul单元，在两个不同的端口上有两个FP add/sub单元，所以如果你用FMA做加法，理论上最多有四个每个时钟周期的 SIMD 添加。
Haswell/Broadwell 有 2/clock FMA，但只有 1/clock FP add/sub（延迟较低）。这对原始代码很有用，not great 对于经过优化以具有大量并行性的代码。这可能就是英特尔在 Skylake 中更改它的原因。
您的 Sandybridge (E5-1620) 和 Nehalem (W5580) CPUs 在不同的端口上有 1/clock FP add/sub、1/clock FP mul。这就是 Haswell 的基础。为什么添加额外的乘法不是大问题：它们可以 运行 与现有的加法并行。 （Sandybridge 的宽度为 256 位，但您在未启用 AVX 的情况下进行编译：使用 -march=native。）
寻找平行度：strength-reduction 任意步幅
您的 compute2 根据前一个值计算下一个 Y 和下一个 Z。即，步幅为 1，即 data[i+1] 所需的值。所以每次迭代都依赖于前一次迭代。
如果将其推广到其他步幅，则可以推进 4、6、8 或更多独立的 Y 和 Z 值，以便它们彼此同步跳跃，彼此独立。 这为编译器 and/or 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 的差异。non-power-of-2 有助于显示步幅的 factors/powers 与 naturally-occurring 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]，然后我们可以根据需要对其进行代数化简。 （对于 compile-time 常量 A、B、C，编译器将 constant-propagate 任一种方式。）
中的中间一栏Z[j]在跨越poly(x)时保持 first-order 的差异，起点为poly(0..stride-1)。这是上面table.Z[j] += second_difference的必要更新是一个标量常数，正如我们从 second-order 差异相同可以看出的那样。通过研究几个不同的
stride和A值（i^2 的系数），我们可以看到它是A * 2 * (stride * stride). （使用像 3 和 5 这样的 non-coprime 值有助于理清事物。）有了更多的代数知识，您就可以象征性地展示这一点。 2 的因子从微积分 PoV 来看是有意义的：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].[= 交互113=]
有 stride 个独立的并行 dep 链供编译器 (SIMD) 和 CPU（流水线执行单元）利用，运行宁 stride比原始速度快 1 倍，直到您在 SIMD FP-add 吞吐量而不是延迟上出现瓶颈，或者如果您的缓冲区不适合 L1d，则存储带宽。 （或者直到编译器面对并且没有很好地/根本没有展开和矢量化这些循环的地步！）
这在实践中是如何编译的：很好地使用 clang
(Godbolt compiler explorer) Clang auto-vectorize 与 stride=16 （4x YMM 向量，每个向量 4 doubles）与clang14 -O3 -march=skylake -ffast-math.
看起来 clang 进一步展开了 2，将 Z[j] += diff2 缩短为 tmp = Z[j] + diff2; / Z[j] += 2*diff2;。这减轻了 Z dep 链上的压力，只留下 Y[j] 来应对 Skylake 上的延迟瓶颈。
因此每个 asm 循环迭代执行 2x 8 vaddpd 指令和 2x 4 存储。循环开销是 add + macro-fused cmp/jne，所以 2 微指令。 （或者对于全局数组，只有一个 add/jne uop，将负索引向上计数为零；它相对于数组的末尾进行索引。）
Skylake 运行 在几乎 1 个存储和每个时钟周期 2x vaddpd 中实现。这是这两件事的最大吞吐量。 front-end 只需要跟上 3 微码/时钟周期多一点，但自 Core2 以来它一直是 4 宽。
Sandybridge-family 中的 uop 缓存使这没问题。 （除非你 运行 进入 Skylake 上的 JCC 勘误表，所以我使用 -mbranches-within-32B-boundaries to have clang pad instructions to avoid that。）
Skylake 的 vaddpd 延迟为 4 个周期，来自 stride=16 的 4 个 dep 链仅勉强足以保持 4 个独立操作在运行中。任何时候 Y[j]+= 没有 运行 循环就绪，就会产生气泡。由于 clang 对 Z[] 链的额外展开，Z[j]+= 可以提前 运行，因此 Z 链可以领先。使用 oldest-ready-first 调度，它往往会稳定到 Yj+= uops 没有冲突的状态，显然，因为它确实 运行 在我的 Skylake 上全速运行。如果我们能让编译器仍然为 stride=32 生成漂亮的 asm，那将留下更多空间，但不幸的是它没有。 （以对奇数尺寸进行更多清理工作为代价。）
奇怪的是，Clang 仅使用 -ffast-math 对其进行矢量化。下面的完整基准测试中的模板版本不需要 --fast-math。源代码被仔细编写为 SIMD-friendly，并按源代码顺序进行数学运算。 （不过，Fast-math 允许 clang 进一步展开 Z 增量。）
编写循环的另一种方法是使用一个内部循环而不是所有 Y 操作，然后是所有 Z 操作。这在下面的基准测试中很好（有时实际上更好），但这里即使使用 -ffast-math 也不会矢量化。对于像这样的 non-trivial 问题，从编译器中获得最佳展开的 SIMD asm 可能是繁琐且不可靠的，并且可能需要一些时间。
我将它包含在 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] 进行手动矢量化。但是这个版本可以为 AArch64 等其他 ISA 向量化。
当 problem-size 不是展开的倍数时，大展开因子的其他缺点是留下更多的清理工作。 （您可以使用 compute1 策略进行清理，让编译器在执行标量之前将其矢量化一两次迭代。）
理论上编译器是 allowe 使用 -ffast-math 为您执行此操作，或者从 compute1 对原始多项式执行 strength-reduction，或者从 compute2 查看步幅如何累积。
但在实践中，这真的很复杂，需要人类自己来做。除非/直到有人开始教编译器如何寻找这样的模式并应用差异方法本身，并选择大步前进！但是，即使使用 -ffast-math，对具有不同 error-accumulation 属性的算法进行大规模重写也可能是不可取的。 （整数不会有任何精度问题，但它仍然是一个复杂的 pattern-match / 替换。）
实验性能结果：
我用 clang13.0.0 在我的台式机 (i7-6700k) 上进行了测试。这实际上是 运行 每个时钟周期 1 个 SIMD 存储，具有多种编译器选项组合（fast-math 或不）和 #if 0 与 #if 1 在内循环策略上。我的基准/测试框架基于@huseyin tugrul buyukisik 的版本，改进后在 rdtsc 指令之间重复更可测量的数量，并使用测试循环来检查多项式简单计算的正确性。
我还让它补偿了核心时钟频率与 "reference" frequency of the TSC read by rdtsc 之间的差异，在我的例子中是 3.9GHz 与 4008 MHz。 （额定最大睿频频率为 4.2GHz，但在 Linux 上 EPP = balance_performance，它只想达到 3.9 GHz。）
源代码 on Godbolt：使用一个内部循环，而不是 3 个单独的 j<16 循环，并且 不 使用 -ffast-math。使用 __attribute__((noinline)) 来防止它内联到重复循环中。选项和源的一些其他变体导致循环内的一些 vpermpd 洗牌。
下面的基准数据来自以前的版本，Z[j] 初始值设定项有问题，但循环 asm 相同。 Godbolt link 现在有一个正确性测试在定时循环之后，它通过了。实际性能在我的桌面上仍然相同，每 double 仅超过 0.25 个周期，即使没有 #if 1 / -ffast-math 允许 clang 额外展开。
$ 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.9GHz)
0.252109 cycles per data element (corrected from ref cycles to core clocks for i7-6700k @ 3.9GHz)
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/时钟，显然没有受到 sub-optimal uop 调度的影响。 （我增加了重复循环，所以程序的大部分时间都花在那里，使整个程序的 perf stat 很有用。）
此优化的目标是以更少的 FLOPS 完成相同的工作，但这也意味着我们实质上是在不使用 FMA 的情况下最大化 Skylake 的 8 FLOP/clock 限制。 （单核 3.9GHz 时 30.58 GFLOP/s）。
Asm 的non-inline 函数(objdump -drwC -Mintel); clang 使用 4 对 Y，Z 对 YMM 向量，并将循环进一步展开 3 倍，使其成为 24KiB 大小的精确倍数，无需清理。注意 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
相关：
- Out of Order Execution, How to Solve True Dependency? - CPU当指令依赖于尚未执行的指令时必须等待。
- 现代微处理器 90 分钟指南！ - out-of-order 执行和管道的一般背景。现代 CPU-style short-vector SIMD 以这种形式存在，以便在不加宽管道的情况下通过单个 CPU 的管道获得更多工作。相比之下，GPU 有很多简单的流水线。
看起来你可以吃蛋糕也可以吃，通过手动将代码并行化为这样的东西：
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;
}
可能不完全像写的那样起作用，但你明白了：展开循环，这样数据相关的路径就可以并行完成。对于正在考虑的机器，4 步展开应该可以实现最佳性能，但是当然，您会得到 hard-coding 软件架构带来的所有有趣的东西。