在我们的cpu中，与加法相比，乘法的速度要慢得多。

这可能在历史上是正确的，对于更简单的低功耗CPU可能仍然是正确的，但如果CPU设计者准备“在问题上扔门”，乘法几乎可以和加法一样快。

现代cpu被设计成同时处理多条指令，通过流水线和多个执行单元的组合。

但问题是数据依赖关系。如果一条指令依赖于另一条指令的结果，那么在它所依赖的那条指令完成之前，它的执行不能开始。

现代cpu试图通过“乱序执行”来解决这个问题。等待数据的指令可以保持队列，而允许执行其他指令。

但是，即使采用这些措施，有时CPU也会耗尽要调度的新工作。

在我们的cpu中，与加法相比，乘法的速度要慢得多。

这可能在历史上是正确的，对于更简单的低功耗CPU可能仍然是正确的，但如果CPU设计者准备“在问题上扔门”，乘法几乎可以和加法一样快。

现代cpu被设计成同时处理多条指令，通过流水线和多个执行单元的组合。

但问题是数据依赖关系。如果一条指令依赖于另一条指令的结果，那么在它所依赖的那条指令完成之前，它的执行不能开始。

现代cpu试图通过“乱序执行”来解决这个问题。等待数据的指令可以保持队列，而允许执行其他指令。

但是，即使采用这些措施，有时CPU也会耗尽要调度的新工作。

理解您所看到的性能差异的关键在于向量化。是的，基于加法的解决方案在其内部循环中只有两条指令，但重要的区别不在于循环中有多少条指令，而在于每条指令执行了多少工作。

在第一个版本中，输出完全依赖于输入:每个数据[i]只是i本身的函数，这意味着每个数据[i]可以以任何顺序计算:编译器可以向前、向后、横向，无论如何，你仍然会得到相同的结果——除非你从另一个线程观察内存，否则你永远不会注意到数据是如何被处理的。

在第二个版本中，输出不依赖于i——它依赖于最后一次循环中的A和Z。

如果我们将这些循环的主体表示为小的数学函数，它们将有非常不同的整体形式:

f(i) - f(Y, Z) -> (di, Y, Z)

在后一种形式中，没有实际依赖于i -你可以计算函数值的唯一方法是通过从函数的最后一次调用中知道前一个Y和Z，这意味着函数形成了一个链-你不能做下一个，直到你做了前一个。

Why does that matter? Because the CPU has vector parallel instructions that each can perform two, four, or even eight arithmetic operations at the same time! (AVX CPUs can do even more in parallel.) That's four multiplies, four adds, four subtracts, four comparisons — four whatevers! So if the output you're trying to compute is only dependent on the input, then you can safely do two, four, or even eight at a time — it doesn't matter if they're forward or backward, since the result is the same. But if the output is dependent on previous computation, then you're stuck doing it in serial form — one at a time.

这就是为什么“较长的”代码在性能方面更胜一筹。尽管它有更多的设置，而且实际上做了更多的工作，但大部分工作是并行完成的:它不仅仅是在循环的每次迭代中计算数据[i]，而是同时计算数据[i]，数据[i+1]，数据[i+2]和数据[i+3]，然后跳转到下一个4个数据的集合。

为了扩展一下我在这里的意思，编译器首先将原始代码转换成这样:

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有特殊的内置指令，可以在同一时间对不同的数据执行多个*或多个+操作，每个操作只需一个时钟周期。

注意快速解决方案addpd和mulpd指令中的字母p，以及慢速解决方案addsd指令中的字母s。这是“添加打包双打”和“增加打包双打”，而不是“添加单个双打”。

不仅如此，它看起来像编译器也部分展开了循环-循环不是每次迭代只执行两个值，而是实际上执行了四个值，并且交错操作以避免依赖关系和停顿，所有这些都减少了汇编代码测试i < 1000的次数。

All of this only works, though, if there are no dependencies between iterations of the loop: If the only thing that determines what happens for each data[i] is i itself. If there are dependencies, if data from the last iteration influences the next one, then the compiler may be so constrained by them that it can't alter the code at all — instead of the compiler being able to use fancy parallel instructions or clever optimizations (CSE, strength reduction, loop unrolling, reordering, et al.), you get out code that's exactly what you put in — add Y, then add Z, then repeat.

但是这里，在代码的第一个版本中，编译器正确地识别出数据中没有依赖关系，并发现它可以并行地完成工作，所以它做到了，这就是所有不同的地方。

The main difference here is loop dependencies. The loop in the second case is dependent—operations in the loop depend on the previous iteration. This means that each iteration can't even start until the previous iteration finishes. In the first case, the loop body is fully independent—everything in the loop body is self contained, depending solely on the iteration counter and constant values. This means that the loop can be computed in parallel—multiple iterations can all work at the same time. This then allows the loop to be trivially unrolled and vectorized, overlapping many instructions.

如果您要查看性能计数器(例如，perf stat ./1)，您将看到第一个循环除了运行得更快之外，每个周期(IPC)还运行了更多的指令。相比之下，第二个循环有更多的依赖周期——当CPU无所事事，等待指令完成，然后才能发出更多指令时。

第一个可能会导致内存带宽瓶颈，特别是如果你让编译器在Sandy Bridge (gcc -O3 -march=native)上使用AVX自动向量化，如果它设法使用256位向量的话。此时，IPC将会中断，特别是对于一个对于L3缓存来说太大的输出数组。

注意:展开和向量化不需要独立的循环——当存在(某些)循环依赖关系时可以执行它们。然而，这样做难度更大，回报也更少。因此，如果你想看到从向量化得到的最大加速，它有助于在可能的情况下删除循环依赖。

通过手动将代码并行化，你似乎可以鱼与熊掌兼得:

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;
}

它可能不完全像所写的那样具有功能性，但您可以理解它的思想:展开循环，以便每个依赖于数据的路径都可以并行执行。对于所考虑的机器，四步展开应该能够实现最大的性能，当然，在软件中对体系结构进行硬编码时，您可以获得所有有趣的东西。

如果只使用加法作为优化，就会错过(新cpu的)乘法管道的所有GFLOPS，而循环携带的依赖关系会因为停止自动向量化而使情况变得更糟。如果它是自向量化的，它将比乘法+加法快得多。并且每个数据的能源效率更高(只添加比mul + add更好)。

另一个问题是数组的末尾会收到更多的舍入错误，因为累积的加法数量。但在非常大的数组(除非数据类型变成float)之前，它不应该是可见的。

当你应用Horner的方案与GCC构建选项(在较新的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

由于mul+add连接到单个FMA，所有浮点运算都是“打包”的向量形式和更少的指令(它是一个两次展开的版本)。每64字节数据16条指令(AVX-512则为128字节)。

Horner方案的另一个优点是，它在FMA指令中计算精度更高，每次循环迭代只有O(1)个操作，因此对于较长的数组，它不会累积那么多错误。

我认为Agner Fog的优化手册中的优化一定来自《雷神之锤III》的快速平方根近似，这在x87中是有意义的，但在SSE中被淘汰了。在那个时候，SIMD不够宽，不能产生太大的差异，而且对于sqrt()和除法来说也很慢，但是有rsqrtss。手册上写着2004年的版权，所以赛扬的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个周期运行。