为何循环展开的霍纳多项式求值在多项式长度N<90时呈现超出预期的准二次型更快性能？

我最近花了不少时间优化用霍纳法求值的多项式代码，针对的是长度N<32的中等规模多项式。我搞了个循环展开的实现，在-O2及以上优化级别下表现特别好，但在N<90时呈现O(N²)的时间复杂度，N超过90后就变成线性了——而且实际性能比理论预测的还快。下面是我遇到的各种奇怪现象和至今没搞懂的问题：

1. 特定长度多项式的“填充优化”反超现象 #

我一开始注意到这个问题，是因为发现某些特定长度N的多项式（只有运行时才知道），用普通霍纳循环的话，居然不如把数组零填充到更长的长度再用更长的求值函数快：

• 在Intel ICX编译器上，当5<N<9时，把系数数组补到10个元素再求值，性能反而更高；
• 在MSVC编译器上，当14<N<18时，补到18或19个元素，速度会更快。

至今我都没搞懂为啥会出现这种“越补越长反而越快”的情况。

2. 时间耗时的台阶式跳变 #

另外，不同编译器的霍纳循环耗时会出现明显的台阶式上升：

• MSVC在N约120-140时，耗时会有个~1秒的台阶式增长；
• Intel编译器则在N约100时出现类似的台阶。
  有意思的是，台阶前后的时间斜率（也就是单位N增加带来的耗时增量）几乎保持不变。

3. 理论预期与实际测试的矛盾：线性vs二次型 #

从理论上来说，霍纳法的Horner<N+1>比Horner<N>只多一条fma指令，而且每条指令都有严格的数据依赖（后一条指令必须等待前一条的计算结果），所以时间复杂度应该和N严格线性相关才对。

但我的基准测试结果却完全打破了这个预期：在N<90的中等规模下，居然呈现出近乎纯二次的时间特性——虽然二次项的系数很小，但在Intel和MS编译器上都能稳定复现（对应测试图里的橙色曲线和绿色点）。

4. 循环展开实现的现状与临界点 #

我的循环展开版实现（命名为Fast_Horner）用了Duff装置做跳转表，只需要一次分支就能处理0<N<67的多项式。虽然搞不懂N<90时的二次型特性，但好在二次项系数极小，而且在N<90时这种O(N²)的实现依然比普通循环快，所以我也就暂时接受了这个结果。

我原本估计这种O(N²)实现和线性实现的临界点在N=140，但实际测试发现，N超过80后就完全变成线性时间复杂度了。

相关测试代码 #

#include <time.h>
#include <immintrin.h>
#include <string>

#define NMAX (1024)
int loop_npoly=1;

// 指数函数的测试系数
const double p64_64[] = {
    1.0, 1.0, 0.5, 0.16666666666666666, 0.041666666666666664,
    0.008333333333333333, 0.001388888888888889, 0.0001984126984126984,
    2.48015873015873e-5, 2.7557319223985922e-6, 2.7557319223981325e-7,
    2.5052108386015238e-8, 2.0876756928192412e-9, 1.6059049047531318e-10,
    1.1470359463231895e-11, 7.671678445264243e-13, 3.4354082706487605e-14,
    6.688002060416919e-14, -2.6687364314970597e-13, 9.777568536777767e-13,
    -3.157534966066221e-12, 9.021749954970395e-12, -2.286148350388897e-11,
    5.1466153063012067e-11, -1.0302563351906206e-10, 1.834070966309193e-10,
    -2.9010550208538685e-10, 4.0688348289884633e-10, -5.041786204830028e-10,
    5.487382269177827e-10, -5.197645055648691e-10, 4.2209407690378606e-10,
    -2.863215889230545e-10, 1.539917550473329e-10, -5.717048308017362e-11,
    5.897550884821984e-12, 9.77179038552452e-12, -7.753082213173471e-12,
    2.266846857748925e-12, 6.363027193335637e-13, -8.936078596015253e-13,
    2.733705604052118e-13, 1.1867318243797467e-13, -1.5074126738677328e-13,
    5.824953355299476e-14, 1.5783407359001581e-15, -1.2251540014686882e-14,
    5.6051913099606426e-15, -8.471028446853139e-16, -4.3384880187935466e-17,
    -1.6191683721154725e-16, 1.6236128793390834e-16, -1.50096740297477e-17,
    -6.664289711051805e-17, 6.403015433527563e-17, -3.504767151862092e-17,
    1.3667212310208614e-17, -4.056415610397064e-18, 9.374879361646055e-19,
    -1.6906743792437812e-19, 2.3477536707437213e-20, -2.4353365725377663e-21,
    1.7828467564653198e-22, -8.23668196718495e-24, 2.0e-26, 1e-28
};

double plog[NMAX], d[NMAX/2];

// 初始化测试用多项式系数
void init(double* plog, int n) {
    double term = 1.0;
    plog[0] = 0;
    for (int i = 1; i <= n; i++) {
        plog[i] = term / i;
        term = -term;
    }
}

// 空函数用于基准测试对比
double null(double x, const double* a) {
    return x;
}

// 标准霍纳循环实现
__inline double hornern(double x, const double* a, int n) {
    double sum = a[n];
    for (int i = n-1; i >= 0; i--) {
        sum = sum * x + a[i];
    }
    return sum;
}

内容来源于stack exchange