Pairwise

成对求和：按位对齐 polars 无过滤的 .sum()

polars，float_sum，Pairwise summation。

polars 中无过滤的整列求和走的是 float_sum，一个固定的成对求和实现。它没有 Kahan 补偿，也不做朴素顺序累加，而是把输入切成固定大小的块，块内按 lane 分组累加，块之间递归对半合并。每一步都是普通 IEEE 754 加法，所以它不承诺更小的误差，只承诺一个确定、SIMD 友好的求和序。浮点加法不满足结合律，同一批数用不同顺序相加，末位 ulp 会不一样；要让 C++ 和 polars 的 .sum() 逐位一致，最省事的办法不是换一个更准的算法，而是把这条求和序原样复刻下来。

整体求和序

入口 sum 先算 \(\text{remainder} = \text{len} \bmod 128\)，把最前面的 \(\text{remainder}\) 个元素拿出来做顺序累加，余下部分的长度一定是 \(128\) 的倍数，交给 pairwiseSum。最后返回 \(\text{main} + \text{rest}\)。这最后一步也有顺序，是 \(\text{main}\) 加 \(\text{rest}\)，反过来写末位可能又不对。长度不足 \(128\) 时整段都落在 \(\text{rest}\) 里，\(\text{main}\) 记 \(0\)。

inline double sum(const double* f, std::size_t len) {
    const std::size_t remainder = len % kRecursionLimit;
    double rest = 0.0;
    for (std::size_t i = 0; i < remainder; ++i) {
        rest += f[i];
    }
    double main = 0.0;
    if (len > remainder) {
        main = pairwiseSum(f + remainder, len - remainder);
    }
    return main + rest;
}

块间递归

pairwiseSum 的递归不是按元素对半分，而是按块数对半分。长度正好是 \(128\) 的倍数且非 \(0\) 时，先看是不是只剩一个块：是就直接进块内求和；否则算 \(\text{blocks} = \frac{\text{len}}{128}\)，左半取 \(\lfloor \frac{\text{blocks}}{2} \rfloor \times 128\) 个元素，右半取剩下的，两边各自递归再相加。

inline double pairwiseSum(const double* f, std::size_t len) {
    if (len == kRecursionLimit) return sumBlock128(f);
    const std::size_t blocks = len / kRecursionLimit;
    const std::size_t leftLen = (blocks / 2) * kRecursionLimit;
    return pairwiseSum(f, leftLen) + pairwiseSum(f + leftLen, len - leftLen);
}

因为 \(\text{blocks} \ge 2\)，左半和右半的长度仍然都是 \(128\) 的倍数，递归到底时叶子永远是一整个 \(128\) 元素块。这个分法比按元素二分少一层偏移量的麻烦，深度大约是块数的对数而不是元素数的对数。

块内归约

块内求和 sumBlock128 是最需要照抄的部分。它把 \(128\) 个元素看成 \(8\) 个 stripe，每个 stripe 宽 \(16\)。先维护一个 \(16\) 宽的累加数组，按 stripe 依次把每个 lane 加到对应的累加槽上。这个顺序等价于上游 chunks_exact(16) 之后逐个向量相加，所以 SIMD 版本和标量版本会得到同一位模式。

inline double sumBlock128(const double* f) {
    double vsum[16] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    for (int s = 0; s < 8; ++s) {
        const double* stripe = f + s * static_cast<int>(kStripe);
        for (int j = 0; j < 16; ++j) {
            vsum[j] += stripe[j];
        }
    }
    return horizontalSum(vsum);
}

随后的水平归约也固定顺序：先做 \(16 \to 8\)，\(t_{j} = v_{j} + v_{8+j}\)；再做 \(8 \to 4\)，\(u_{j} = t_{j} + t_{4+j}\)；最后返回 \((u_{0} + u_{2}) + (u_{1} + u_{3})\)。

inline double horizontalSum(const double* v) {
    double t[8];
    for (int j = 0; j < 8; ++j) {
        t[j] = v[j] + v[8 + j];
    }
    double u[4];
    for (int j = 0; j < 4; ++j) {
        u[j] = t[j] + t[4 + j];
    }
    return (u[0] + u[2]) + (u[1] + u[3]);
}

和 Kahan 的边界

“差不多”的求和算法都会差 ulp，原因还是结合律：\((a + b) + c\) 和 \(a + (b + c)\) 的舍入点不同。拿 \(4801\) 段真实明细验证过，按上面的序列逐位模拟能全部对上；同样的数据换成普通顺序累加会差 \(593\) 段，Kahan 差 \(646\) 段，普通 pairwise 树差 \(593\) 段。所以这个实现的目标不是精度，而是进入 polars 那条求和序。

工程侧用起来很薄：把需要对齐的样本先收集进一个数组，收尾时调一次 sum，数组为空时入口自然返回 \(0\)，不用额外判空。和 Kahan 那条路径的关系也要分清：Kahan 对应的是带 filter 的 group_by 路径，无 filter 的 .sum() 才是这里复刻的成对求和。两者对同一批数据可能差 \(1\) 到 \(2\) 个 ulp，平时看不出问题；一旦这个和进入两个大数相减的表达式，这点差异可能被放大数千倍，对拍就过不去。

抛开对齐看复杂度

如果不考虑 py 对齐，单把它当成一个 C++ 求和实现，加法数量其实比顺序累加更多，但依赖链更短。这两个量要分开算。

记 \(B = \lfloor \frac{n}{128} \rfloor\)，\(r = n \bmod 128\)，主块部分长度是 \(128B\)。先数加法。主块的每个元素在 lane 累加里恰好被加一次，共 \(128B\) 次；每个块的水平归约是 \(16 \to 8 \to 4 \to 2 \to 1\)，合计 \(8 + 4 + 3 = 15\) 次，共 \(15B\) 次；块间合并是一棵有 \(B\) 个叶子的二叉树，内部节点 \(B - 1\) 个，对应 \(B - 1\) 次加法；余数串行累加 \(r\) 次；最后 \(\text{main} + \text{rest}\) 一次。\(B \ge 1\) 时总加法数是

\[128B + 15B + (B - 1) + r + 1 = 144B + r \]

当 \(n\) 是 \(128\) 的倍数时 \(r = 0\)，总加法数正好是 \(\frac{144}{128}n = \frac{9}{8}n\)，渐近比顺序累加的 \(n - 1\) 多 \(12.5\%\)。\(n = 10^{7} = 128 \times 78125\)，总加法数是 \(144 \times 78125 = 1.125 \times 10^{7}\)，顺序累加是 \(9999999\) 次。\(n < 128\) 时没有主块，只有 \(n\) 次余数累加加收尾一次，是另一种边界。

再看依赖深度，也就是一条加法链最多串多少次。块内每个 lane 先做 \(8\) 次串行加法，水平归约 \(4\) 层，所以一个块的深度是 \(8 + 4 = 12\) 次加法，不是 \(11\)：最后 \((u_{0} + u_{2}) + (u_{1} + u_{3})\) 里 \((u_{0} + u_{2})\) 和 \((u_{1} + u_{3})\) 可以并行，但它们到最终结果还差一层。这里的深度只数浮点加法：return 不会额外产生一次加法，函数是否内联也不会改变括号结构，所以 horizontalSum 返回值的数据依赖深度就是 \(12\)；未内联时的调用和返回在真实墙钟上有额外延迟，但那是另一笔开销。块与块之间，合并树有 \(B\) 个叶子，深度是 \(\lceil \log_2 B \rceil\)。余数那 \(r\) 个元素是另一条串行链。最后 \(\text{main} + \text{rest}\) 再加一层，所以关键路径是

\[\max(12 + \lceil \log_2 B \rceil, r) + 1 \]

\(n = 10^{7}\) 时 \(r = 0\)，\(\lceil \log_2 78125 \rceil = 17\)，关键路径是 \(12 + 17 + 1 = 30\) 次加法；顺序累加的依赖链是 \(n - 1 \approx 10^{7}\) 次。这里的 \(30\) 只是深度，不是总加法数：每层递归里的块可以并行处理，但单核总时间仍由总加法数决定，所以总量是 \(O(n)\)，\(\log\) 只出现在深度里。把单块成本和 \(\log\) 相乘会把总量和深度混在一起。

深度短的好处是指令级并行和乱序执行能利用 \(16\) 个 lane 的独立性。总的加法多出八分之一，换来的关键路径从 \(O(n)\) 降到 \(O(\log n)\)；单核墙钟时间不一定按这个比例缩短，因为加法和读内存的吞吐量仍是上限，但至少不会卡在一条长度 \(n\) 的依赖链上。误差方面，它没有补偿项，不跟 Kahan 比最坏误差；树形归约通常把误差界从 \(O(\varepsilon n)\) 降到 \(O(\varepsilon \log n)\) 量级，这里先固定 \(128\) 块、块内 \(16\) lane，增长阶仍是 \(\log n\) 级。缓存上，块内按 stripe 访问，但一个 stripe 是连续的 \(16\) 个 double，\(128\) 个元素整体在一小段连续内存里，实际多不了几次 cache miss。

和按元素二分的完全 pairwise 树比，它的深度同阶，都是 \(O(\log n)\)，但总加法多出约 \(\frac{1}{8}n\)，换来的是固定的 \(16\) lane 布局。和 Kahan 比，Kahan 每个元素大约 \(4\) 次浮点运算且串行，总量约 \(4n\)、深度也约 \(4n\)；这里总量是 \(\frac{9}{8}n\)、深度是 \(O(\log n)\)，代价是没有补偿项，误差界从 Kahan 接近 \(O(\varepsilon)\) 退到 \(O(\varepsilon \log n)\)。

递归那部分在 C++ 里也可以改成迭代，但迭代版本必须保持完全相同的加法括号结构，否则位级对齐又会漂；前提是编译时不开 -ffast-math，否则整个括号结构都可能被重结合改掉。

完整实现

完整的成对求和标量实现：

点击查看代码

#include <cstddef>

constexpr std::size_t kStripe = 16;
constexpr std::size_t kRecursionLimit = 128;

// 16 宽水平归约，顺序与上游 vector_horizontal_sum 一致。
inline double horizontalSum(const double* v) {
    double t[8];
    for (int j = 0; j < 8; ++j) {
        t[j] = v[j] + v[8 + j];
    }
    double u[4];
    for (int j = 0; j < 4; ++j) {
        u[j] = t[j] + t[4 + j];
    }
    return (u[0] + u[2]) + (u[1] + u[3]);
}

// 128 元素块求和：8 个 16 宽 stripe 逐 lane 累加。
inline double sumBlock128(const double* f) {
    double vsum[16] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    for (int s = 0; s < 8; ++s) {
        const double* stripe = f + s * static_cast<int>(kStripe);
        for (int j = 0; j < 16; ++j) {
            vsum[j] += stripe[j];
        }
    }
    return horizontalSum(vsum);
}

// 块间递归二分，要求 len 为 128 的倍数且非 0。
inline double pairwiseSum(const double* f, std::size_t len) {
    if (len == kRecursionLimit) {
        return sumBlock128(f);
    }
    const std::size_t blocks = len / kRecursionLimit;
    const std::size_t leftLen = (blocks / 2) * kRecursionLimit;
    return pairwiseSum(f, leftLen) + pairwiseSum(f + leftLen, len - leftLen);
}

// 入口：前 remainder 个元素顺序累加，其余进入成对求和。
inline double sum(const double* f, std::size_t len) {
    const std::size_t remainder = len % kRecursionLimit;
    double rest = 0.0;
    for (std::size_t i = 0; i < remainder; ++i) {
        rest += f[i];
    }
    double main = 0.0;
    if (len > remainder) {
        main = pairwiseSum(f + remainder, len - remainder);
    }
    return main + rest;
}