太长不看

template <auto First, auto Last>
constexpr auto constexpr_for = [](auto &&f) {
    [&]<auto I>(this auto &&self) {
        if constexpr (I < Last) {
            f.template operator()<I>();
            self.template operator()<I + 1>();
        }
    }.template operator()<First>();
};

简单结论：

• 使用 SIMD intrinsic 编程时，编译器 -O3 并不提供充分的 ILP 优化，因此需要显式声明。
• 使用 constexpr_for 即可做到显式 ILP。
• 在编译器 -O3 优化质量较低时，显式 ILP 可以掩盖性能问题。
• 注意 std::ranges 的抽象税。

---

下面会用 SIMD 求和的示例程序来说明。限于篇幅，本文不会解释 ILP 的基本概念。

SIMD sum

#include <x86intrin.h>
#include <algorithm>
#include <ranges>

// 使用 C++23 的 range_adaptor_closure 进行快速定制适配器
// 从而方便构造 SIMD 的批量视图（非完整块会丢弃）
template <size_t Lane>
struct simdify_t: std::ranges::range_adaptor_closure<simdify_t<Lane>> {
    constexpr auto operator()(auto &&view) const noexcept {
        auto size = std::ranges::size(view);
        return std::forward<decltype(view)>(view)
            | std::views::stride(Lane) // 每 Lane 个元素取一个
            | std::views::take(size / Lane); // 总共只取 size / Lane 个元素
    }
};
template <auto i>
constexpr simdify_t<i> simdify;

int sum_avx2(std::ranges::range auto &&rng) {
    constexpr auto lane = sizeof(__m256i) / sizeof(int);
    __m256i partial_sum {};

    auto simd_view = rng
                   | simdify<lane>;
    for(auto &&simd_v : simd_view) {
        auto addr = &simd_v;
        auto loadu = _mm256_loadu_si256((__m256i*)addr);
        partial_sum = _mm256_add_epi32(partial_sum, loadu);
    }

    auto scalar_view = rng
                     | std::views::drop(lane * std::ranges::size(simd_view));
    int sum = std::ranges::fold_left(scalar_view, 0, std::plus());

    int temp[lane];
    // Clang 编译器会优化为寄存器内归约操作
    _mm256_storeu_si256((__m256i*)std::ranges::data(temp), partial_sum);
    sum = std::ranges::fold_left(temp, sum, std::plus());

    return sum;
}

我们看一个 AVX2 版本的 SIMD 求和实现。非常直接，示例程序将输入的整个 range 拆分为 SIMD 视图和 scalar 视图并分别计算，计算完成后合并求和结果。

SIMD+ILP sum

#include <x86intrin.h>
#include <algorithm>
#include <ranges>

// 显式 ILP 的工具
template <auto First, auto Last>
constexpr auto constexpr_for = [](auto &&f) {
    [&]<auto I>(this auto &&self) {
        if constexpr (I < Last) {
            f.template operator()<I>();
            self.template operator()<I + 1>();
        }
    }.template operator()<First>();
};

template <size_t Lane>
struct simdify_t: std::ranges::range_adaptor_closure<simdify_t<Lane>> {
    constexpr auto operator()(auto &&view) const noexcept {
        auto size = std::ranges::size(view);
        return std::forward<decltype(view)>(view)
            | std::views::stride(Lane)
            | std::views::take(size / Lane);
    }
};
template <auto i>
constexpr simdify_t<i> simdify;

template <size_t ILP = 4>
int sum_avx2_ilp(std::ranges::range auto &&rng) {
    constexpr auto lane = sizeof(__m256i) / sizeof(int);
    constexpr auto bulk = lane * ILP;
    __m256i partial_sum[ILP] {};
    auto process_simd = [&](auto konstexpr_for, auto simd_view) {
        for(auto &&simd_v : simd_view) {
            // ILP 展开，最高一次并行处理 lane * ILP * 4 字节的数据
            konstexpr_for([&, addr = &simd_v]<auto Index> {
                auto &partial = partial_sum[Index];
                auto loadu = _mm256_loadu_si256(
                    (__m256i*)(addr + Index * lane));
                partial = _mm256_add_epi32(partial, loadu);
            });
        }
    };

    // SIMD+ILP
    auto bulk_simd_view = rng
                        | simdify<lane>
                        | simdify<ILP>; // 这里就知道定制 simdify 的妙处了
    process_simd(constexpr_for<0, ILP>, bulk_simd_view);

    // SIMD
    auto single_simd_view = rng
                          | std::views::drop(bulk * std::ranges::size(bulk_simd_view))
                          | simdify<lane>;
    process_simd(constexpr_for<0, 1>, single_simd_view);

    auto scalar_view = rng
                     | std::views::drop(bulk * std::ranges::size(bulk_simd_view))
                     | std::views::drop(lane * std::ranges::size(single_simd_view));
    int sum = std::ranges::fold_left(scalar_view, 0, std::plus());

    constexpr_for<1, ILP>([&]<size_t Index> {
        partial_sum[0] = _mm256_add_epi32(partial_sum[0], partial_sum[Index]);
    });
    int temp[lane];
    _mm256_storeu_si256((__m256i*)std::ranges::data(temp), partial_sum[0]);
    sum = std::ranges::fold_left(temp, sum, std::plus());

    return sum;
}

显式 ILP 的实现如上。在前面例子的基础上扩展，将输入 range 拆分为三个视图：bulk SIMD / SIMD / scalar，同样是只有无法完整覆盖的块才会被下一级的视图处理。

也就是说 ILP 版本的思路是扩展到流水线的宽度，而不是单个向量寄存器的宽度：处理器端口这么多，大于一总比等于一好。

动机和跑分

可是为什么需要特意写 constexpr_for，编译器不会优化吗？

与其理论，不如直接跑分。我们使用 Clang-20 编译器和 Google Benchmark，开启 -O3 优化和 -march=znver3 目标架构（面向我的 Zen 3 主机）。

代码折叠

#include <benchmark/benchmark.h>
#include <x86intrin.h>
#include <type_traits>
#include <algorithm>
#include <ranges>
#include <array>
#include <random>

// ...此前的代码

// ----------------------------------------------------------------------------
// Google Benchmark
// ----------------------------------------------------------------------------

template <size_t Size>
const auto& generate_data() {
    static const auto arr = [] {
        std::array<int, Size> result;
        std::mt19937 gen{std::random_device{}()};
        // No overflow.
        std::uniform_int_distribution dist(-10, 10);
        std::ranges::generate(result, [&] { return dist(gen); });
        return result;
    } ();
    return arr;
}

template <size_t Size>
void BM_run(benchmark::State& state, auto &&func) {
    const auto &rng = generate_data<Size>();

    for(auto _ : state) {
        auto res = func(rng);
        benchmark::DoNotOptimize(res);
    }

    state.SetBytesProcessed(int64_t(state.iterations()) * int64_t(Size) * sizeof(int));
}

template <size_t ...Is>
void register_tests(std::integer_sequence<size_t, Is...>,
                    auto name, auto func) {
    (benchmark::RegisterBenchmark(
        std::string(name) + "/" + std::to_string(Is),
        [func](benchmark::State& state) {
            BM_run<Is>(state, func);
        }
    ), ...);
}

int main(int argc, char** argv) {
    std::integer_sequence<size_t,
        35,
        350,
        3502,
        35023,
        350234
    > seq;

    auto register_test = [seq](auto name, auto func) {
        register_tests(seq, name, func);
    };

    register_test("BM_sum_avx2",
        [](const auto& r) { return sum_avx2(r); });

    register_test("BM_sum_avx2_ilp<4>",
        [](const auto& r) { return sum_avx2_ilp<4>(r); });

    benchmark::Initialize(&argc, argv);
    benchmark::RunSpecifiedBenchmarks();
    benchmark::Shutdown();
    return 0;
}

Run on (16 X 3193.93 MHz CPU s)
CPU Caches:
  L1 Data 32 KiB (x8)
  L1 Instruction 32 KiB (x8)
  L2 Unified 512 KiB (x8)
  L3 Unified 16384 KiB (x1)
Load Average: 0.29, 0.20, 0.11
------------------------------------------------------------------------------------
Benchmark                          Time             CPU   Iterations UserCounters...
------------------------------------------------------------------------------------
BM_sum_avx2/35                  1.00 ns         1.00 ns    690551079 bytes_per_second=130.031Gi/s
BM_sum_avx2/350                 53.4 ns         53.4 ns     12520116 bytes_per_second=24.3976Gi/s
BM_sum_avx2/3502                 621 ns          621 ns      1098424 bytes_per_second=21.007Gi/s
BM_sum_avx2/35023               6116 ns         6116 ns       109537 bytes_per_second=21.3326Gi/s
BM_sum_avx2/350234             63115 ns        63113 ns        10968 bytes_per_second=20.673Gi/s
BM_sum_avx2_ilp<4>/35          0.978 ns        0.978 ns    697817742 bytes_per_second=133.345Gi/s
BM_sum_avx2_ilp<4>/350          14.2 ns         14.2 ns     48383810 bytes_per_second=91.7483Gi/s
BM_sum_avx2_ilp<4>/3502          155 ns          155 ns      4560351 bytes_per_second=84.2754Gi/s
BM_sum_avx2_ilp<4>/35023        1457 ns         1457 ns       475022 bytes_per_second=89.5557Gi/s
BM_sum_avx2_ilp<4>/350234      17217 ns        17210 ns        41257 bytes_per_second=75.8117Gi/s

总之，我们确实得到了数倍的并行性能提升，也验证了 ILP 策略的有效性。

汇编结果

avx::test_medium(std::array<int, 3502ul>):
        // ...
        vpaddd  (%rsi), %ymm0, %ymm0
        // ...

avx_ilp::test_medium(std::array<int, 3502ul>):
        // ...
        vpaddd  (%rcx,%rdx), %ymm1, %ymm1
        vpaddd  32(%rcx,%rdx), %ymm2, %ymm2
        vpaddd  64(%rcx,%rdx), %ymm0, %ymm0
        vpaddd  96(%rcx,%rdx), %ymm3, %ymm3
        // ...

我们可以在 godbolt 快速查验编译器的生成质量。一个很明确的结论是：SIMD intrinsic 确实不会在指令并行这方面有优化，至少现在（2025 年）不行。

NOTES:

• intrinsic 有一定程度的优化，但是下面会指出它相对常规指令的优化会怂很多。
• 能不能用 for 替代 constexpr_for？结论是能但是不可靠，Clang 在 -O2 级别就没法优化到位。
• 关于 constexpr_for 的实现方式有一点话题无关的讨论，见下方折叠。

constexpr_for 实现方式

// 一种更直白的写法
template <auto First, auto Last>
constexpr void constexpr_for(auto &&f) {
    if constexpr (First < Last) {
        f.template operator()<First>();
        constexpr_for<First + 1, Last>(f);
    }
}

其实 constexpr 可以写得更加直白，但是这样的实现是函数而不是函数对象。

// 如果这么写，结合函数版本的 constexpr_for……
auto process_simd = [&]<size_t Width>(auto simd_view) {
    for(auto &&simd_v : simd_view) {
        constexpr_for<0, Width>([]<auto> {
            // ...
        });
    }
};

// 这个时候需要
process_simd.template operator()<ILP>(bulk_simd_view);

而又由于 process_simd 是一个 lambda，它的模板需要传入类型信息，如果不想调用用户不友好的 .template operator() 的话需要很曲折的做法，而改为传入一个带类型信息的对象很容易。

template <auto First, auto Last>
constexpr auto constexpr_for = [](auto &&f) {
    [&]<auto ...I>(std::index_sequence<I...>) {
        (f.template operator()<First + I>(), ...);
    }(std::make_index_sequence<Last - First>{});
};

另一个符合要求的做法也可以不使用 C++23 deducing this 特性。这种实现需要额外的头文件（现实中通常不需要）引入 std::index sequence，并且可读性我觉得不太好，怕被喷就没敢在文章用……当然经典的仿函数类就不用说了。

目前的实现缺陷是，要求的编译器版本比较高，比如 GCC 要开到 14 及以上版本。

也就是说，程序中的 __m256i partial_sum 是只能作为一个 ymm 寄存器去使用。而 constexpr_for 显式地以 ILP 的形式去展开，所以 ILP=4 时会明确使用 4 路并行的 ymm 寄存器。

基线差距

// 在 Google benchmark 注册一个新的测试
register_test("BM_sum_std_fold_left",
    [](const auto& r) { return std::ranges::fold_left(r, 0, std::plus()); });

我们再以 std::ranges::fold_left 作为基线对比，它的实现非常朴素：libstdc++/ranges_algo.h。

--------------------------------------------------------------------------------------
Benchmark                            Time             CPU   Iterations UserCounters...
--------------------------------------------------------------------------------------
BM_sum_avx2/35                   0.985 ns        0.985 ns    694303937 bytes_per_second=132.394Gi/s
BM_sum_avx2/350                   53.5 ns         53.5 ns     12573037 bytes_per_second=24.3865Gi/s
BM_sum_avx2/3502                   612 ns          612 ns      1133219 bytes_per_second=21.3033Gi/s
BM_sum_avx2/35023                 6022 ns         6022 ns       114100 bytes_per_second=21.6667Gi/s
BM_sum_avx2/350234               62198 ns        62195 ns        11085 bytes_per_second=20.9778Gi/s
BM_sum_avx2_ilp<4>/35            0.974 ns        0.974 ns    710236295 bytes_per_second=133.852Gi/s
BM_sum_avx2_ilp<4>/350            14.2 ns         14.2 ns     49325710 bytes_per_second=91.594Gi/s
BM_sum_avx2_ilp<4>/3502            142 ns          142 ns      4957734 bytes_per_second=92.1278Gi/s
BM_sum_avx2_ilp<4>/35023          1481 ns         1481 ns       467864 bytes_per_second=88.086Gi/s
BM_sum_avx2_ilp<4>/350234        17261 ns        17256 ns        40704 bytes_per_second=75.6079Gi/s
BM_sum_std_fold_left/35          0.988 ns        0.988 ns    691087391 bytes_per_second=131.98Gi/s
BM_sum_std_fold_left/350          8.63 ns         8.63 ns     78505005 bytes_per_second=151.164Gi/s
BM_sum_std_fold_left/3502         90.9 ns         90.9 ns      7649494 bytes_per_second=143.574Gi/s
BM_sum_std_fold_left/35023        1026 ns         1025 ns       679782 bytes_per_second=127.277Gi/s
BM_sum_std_fold_left/350234      15973 ns        15972 ns        48874 bytes_per_second=81.6897Gi/s

可以看出，fold_left 优化明显更强，我们有必要再看一眼汇编去确认它是怎么做的。

rng::test_medium(std::array<int, 3502ul>):
        // ...
.LBB10_1:
        vpaddd  -992(%rcx,%rax,4), %ymm0, %ymm0
        vpaddd  -960(%rcx,%rax,4), %ymm1, %ymm1
        vpaddd  -928(%rcx,%rax,4), %ymm2, %ymm2
        vpaddd  -896(%rcx,%rax,4), %ymm3, %ymm3
        vpaddd  -864(%rcx,%rax,4), %ymm0, %ymm0
        vpaddd  -832(%rcx,%rax,4), %ymm1, %ymm1
        vpaddd  -800(%rcx,%rax,4), %ymm2, %ymm2
        vpaddd  -768(%rcx,%rax,4), %ymm3, %ymm3
        vpaddd  -736(%rcx,%rax,4), %ymm0, %ymm0
        vpaddd  -704(%rcx,%rax,4), %ymm1, %ymm1
        vpaddd  -672(%rcx,%rax,4), %ymm2, %ymm2
        vpaddd  -640(%rcx,%rax,4), %ymm3, %ymm3
        vpaddd  -608(%rcx,%rax,4), %ymm0, %ymm0
        vpaddd  -544(%rcx,%rax,4), %ymm2, %ymm4
        vpaddd  -576(%rcx,%rax,4), %ymm1, %ymm1
        vpaddd  -512(%rcx,%rax,4), %ymm3, %ymm5
        vpaddd  -480(%rcx,%rax,4), %ymm0, %ymm2
        vpaddd  -448(%rcx,%rax,4), %ymm1, %ymm3
        vpaddd  -416(%rcx,%rax,4), %ymm4, %ymm0
        vpaddd  -384(%rcx,%rax,4), %ymm5, %ymm1
        cmpq    $3577, %rax
        je      .LBB10_3
        vpaddd  -352(%rcx,%rax,4), %ymm2, %ymm2
        vpaddd  -320(%rcx,%rax,4), %ymm3, %ymm3
        vpaddd  -288(%rcx,%rax,4), %ymm0, %ymm0
        vpaddd  -256(%rcx,%rax,4), %ymm1, %ymm1
        vpaddd  -224(%rcx,%rax,4), %ymm2, %ymm2
        vpaddd  -160(%rcx,%rax,4), %ymm0, %ymm4
        vpaddd  -192(%rcx,%rax,4), %ymm3, %ymm3
        vpaddd  -128(%rcx,%rax,4), %ymm1, %ymm5
        vpaddd  -96(%rcx,%rax,4), %ymm2, %ymm0
        vpaddd  -64(%rcx,%rax,4), %ymm3, %ymm1
        vpaddd  -32(%rcx,%rax,4), %ymm4, %ymm2
        vpaddd  (%rcx,%rax,4), %ymm5, %ymm3
        addq    $256, %rax
        jmp     .LBB10_1
        // ...

其实 fold_left 不仅会并行使用多路寄存器，而且还进一步使用（非常暴力的）循环展开。

NOTES:

• 为什么不用 std::reduce 而用 stdr::fold_left？因为前者的跑分表现非常菜。
• 如果是 std::reduce+std::execution::unseq 的话，表现和 fold_left 类似。
• 这里说的是 Clang 搭配 libstdc++ 的情况；而 libc++ 表现没有劣化，但是仍不强于 fold_left。
• 目前 libc++ 场景仍需要启用 -fexperimental-library 才能完成测试。
• 目前 libc++ 在 ranges 支持上不够完善，不足以完成所有测试。

进一步优化

现在以 fold_left 为目标改进原有的程序。

auto process_simd = [&](auto konstexpr_for, auto simd_view) {
    #pragma clang loop unroll_count(4)
    for(auto &&simd_v : simd_view) {
        // ...
    }
}

首先能想到的是编译器特定的 unroll 选项，结论是对 ranges 起负作用。感兴趣可以自行看下汇编，编译器会生成一堆非常混乱的指令。

// ILP 版本类似，但是问题更加严重
avx::test_medium(std::array<int, 3502ul>):
        // ...
.LBB1_1:
        vpaddd  (%rsi), %ymm0, %ymm0
        movq    %rax, %rdi
        subq    %rsi, %rdi
        leaq    32(%rsi), %r8
        sarq    $2, %rdi
        addq    $-9, %rdi
        cmpq    $-8, %rdi
        cmovaeq %rax, %r8
        cmpq    %rsi, %rax
        cmovneq %r8, %rsi
        cmpq    %rdx, %rsi
        jne     .LBB1_1

进一步调查 AVX/ILP 版本的汇编，发现另一个问题是：它总是在计算莫名其妙的边界，一个合理推测是 ranges 本身的编译器优化不到位。也就是说哪怕是 ILP 版本也是有大量无用指令的。一个简单的解决方案是使用最朴素的 for 循环，不交抽象税。

template <size_t Lane>
struct simdify_t: std::ranges::range_adaptor_closure<simdify_t<Lane>> {
    constexpr auto operator()(auto &&r) const noexcept {
        auto v = std::forward<decltype(r)>(r) | std::views::all;
        auto new_size = std::ranges::size(v) / Lane;
        return std::views::transform(
            std::views::iota(size_t{0}, new_size),
            [v](size_t i) -> decltype(auto) {
                return v[i * Lane];
            }
        );
    }
};
template <auto i>
constexpr simdify_t<i> simdify;

当然个人原因我不打算放弃 ranges，重写一份让编译器能更加优化友好的版本（DoD 风格？）。这样我们程序的核心算法完全不需要修改。

avx::test_medium(std::array<int, 3502ul>):
        movl    $704, %eax
        leaq    8(%rsp), %rcx
        vpxor   %xmm0, %xmm0, %xmm0
.LBB1_1:
        vpaddd  -704(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -672(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -640(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -608(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -576(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -544(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -512(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -480(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -448(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -416(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -384(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -352(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -320(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -288(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -256(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -224(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -192(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -160(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -128(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -96(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -64(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -32(%rcx,%rax), %ymm0, %ymm0
        vpaddd  (%rcx,%rax), %ymm0, %ymm0
        addq    $736, %rax
        cmpq    $14688, %rax
        jne     .LBB1_1
        // ...

avx_ilp::test_medium(std::array<int, 3502ul>):
        movl    $992, %eax
        leaq    8(%rsp), %rcx
        vpxor   %xmm0, %xmm0, %xmm0
        vpxor   %xmm1, %xmm1, %xmm1
        vpxor   %xmm2, %xmm2, %xmm2
        vpxor   %xmm3, %xmm3, %xmm3
.LBB5_1:
        vpaddd  -992(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -960(%rcx,%rax), %ymm1, %ymm1
        vpaddd  -928(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -896(%rcx,%rax), %ymm3, %ymm3
        vpaddd  -864(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -832(%rcx,%rax), %ymm1, %ymm1
        vpaddd  -800(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -768(%rcx,%rax), %ymm3, %ymm3
        vpaddd  -736(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -704(%rcx,%rax), %ymm1, %ymm1
        vpaddd  -672(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -640(%rcx,%rax), %ymm3, %ymm3
        vpaddd  -608(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -544(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -576(%rcx,%rax), %ymm1, %ymm1
        vpaddd  -512(%rcx,%rax), %ymm3, %ymm4
        vpaddd  -448(%rcx,%rax), %ymm1, %ymm3
        vpaddd  -480(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -416(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -384(%rcx,%rax), %ymm4, %ymm1
        cmpq    $14304, %rax
        je      .LBB5_3
        vpaddd  -352(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -320(%rcx,%rax), %ymm3, %ymm3
        vpaddd  -288(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -256(%rcx,%rax), %ymm1, %ymm1
        vpaddd  -224(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -160(%rcx,%rax), %ymm0, %ymm0
        vpaddd  -192(%rcx,%rax), %ymm3, %ymm3
        vpaddd  -128(%rcx,%rax), %ymm1, %ymm4
        vpaddd  -64(%rcx,%rax), %ymm3, %ymm1
        vpaddd  -96(%rcx,%rax), %ymm2, %ymm2
        vpaddd  -32(%rcx,%rax), %ymm0, %ymm0
        vpaddd  (%rcx,%rax), %ymm4, %ymm3
        addq    $1024, %rax
        jmp     .LBB5_1
        // ...

修改后可以确认循环展开和多路寄存器都是有的，而且不再有诡异的边界计算。一次性解决了两个问题。简单对比，可以认为 ILP 版本和 fold_left 基本一致，AVX 版本去掉了多路 YMM 并行（只依靠硬件调度）。

--------------------------------------------------------------------------------------
Benchmark                            Time             CPU   Iterations UserCounters...
--------------------------------------------------------------------------------------
BM_sum_avx2/35                   0.968 ns        0.968 ns    670405544 bytes_per_second=134.668Gi/s
BM_sum_avx2/350                   9.82 ns         9.82 ns     68731562 bytes_per_second=132.756Gi/s
BM_sum_avx2/3502                  96.4 ns         96.4 ns      7216600 bytes_per_second=135.349Gi/s
BM_sum_avx2/35023                 1055 ns         1055 ns       651462 bytes_per_second=123.717Gi/s
BM_sum_avx2/350234               15576 ns        15575 ns        49096 bytes_per_second=83.7717Gi/s
BM_sum_avx2_ilp<4>/35            0.952 ns        0.952 ns    715878822 bytes_per_second=136.893Gi/s
BM_sum_avx2_ilp<4>/350            9.79 ns         9.79 ns     71926717 bytes_per_second=133.126Gi/s
BM_sum_avx2_ilp<4>/3502           85.1 ns         85.1 ns      8083910 bytes_per_second=153.384Gi/s
BM_sum_avx2_ilp<4>/35023          1047 ns         1045 ns       650510 bytes_per_second=124.797Gi/s
BM_sum_avx2_ilp<4>/350234        15939 ns        15938 ns        48235 bytes_per_second=81.863Gi/s
BM_sum_std_fold_left/35           1.05 ns         1.05 ns    663313893 bytes_per_second=124.056Gi/s
BM_sum_std_fold_left/350          8.74 ns         8.74 ns     78781858 bytes_per_second=149.122Gi/s
BM_sum_std_fold_left/3502         91.0 ns         91.0 ns      7675929 bytes_per_second=143.435Gi/s
BM_sum_std_fold_left/35023        1115 ns         1113 ns       630869 bytes_per_second=117.201Gi/s
BM_sum_std_fold_left/350234      16219 ns        16218 ns        45655 bytes_per_second=80.449Gi/s

另外可以看出，此时 AVX 版本已经和 ILP 版本没有性能差异。也即是说，即使只用单个 YMM 寄存器，Zen 3 处理器的硬件调度仍然是非常可靠的，在编译器优化到位的前提下，不需要显式手写 ILP 设计。

浮点 benchmark

--------------------------------------------------------------------------------------
Benchmark                            Time             CPU   Iterations UserCounters...
--------------------------------------------------------------------------------------
BM_sum_avx2/35                    1.46 ns         1.46 ns    472044277 bytes_per_second=89.2223Gi/s
BM_sum_avx2/350                   12.2 ns         12.2 ns     55809571 bytes_per_second=106.701Gi/s
BM_sum_avx2/3502                   274 ns          274 ns      2538655 bytes_per_second=47.5458Gi/s
BM_sum_avx2/35023                 3090 ns         3089 ns       229860 bytes_per_second=42.2333Gi/s
BM_sum_avx2/350234               31352 ns        31351 ns        21813 bytes_per_second=41.6167Gi/s
BM_sum_avx2_ilp<4>/35             1.38 ns         1.38 ns    492076739 bytes_per_second=94.4973Gi/s
BM_sum_avx2_ilp<4>/350            12.0 ns         12.0 ns     57376116 bytes_per_second=108.563Gi/s
BM_sum_avx2_ilp<4>/3502           89.9 ns         89.9 ns      7725717 bytes_per_second=145.178Gi/s
BM_sum_avx2_ilp<4>/35023          1044 ns         1044 ns       651714 bytes_per_second=125.029Gi/s
BM_sum_avx2_ilp<4>/350234        15670 ns        15670 ns        48038 bytes_per_second=83.2645Gi/s
BM_sum_std_fold_left/35           1.39 ns         1.39 ns    503907127 bytes_per_second=93.8597Gi/s
BM_sum_std_fold_left/350          10.6 ns         10.6 ns     66140767 bytes_per_second=123.357Gi/s
BM_sum_std_fold_left/3502         93.7 ns         93.7 ns      7384452 bytes_per_second=139.258Gi/s
BM_sum_std_fold_left/35023        1053 ns         1053 ns       641127 bytes_per_second=123.88Gi/s
BM_sum_std_fold_left/350234      15627 ns        15627 ns        47696 bytes_per_second=83.4937Gi/s

// 需要明确启用 -ffast-math 编译，否则标准库吞吐量是个位数
#include <benchmark/benchmark.h>
#include <x86intrin.h>
#include <algorithm>
#include <ranges>
#include <array>
#include <vector>
#include <random>

template <auto First, auto Last>
constexpr auto constexpr_for = [](auto &&f) {
    [&]<auto I>(this auto &&self) {
        if constexpr (I < Last) {
            f.template operator()<I>();
            self.template operator()<I + 1>();
        }
    }.template operator()<First>();
};

template <size_t Lane>
struct simdify_t: std::ranges::range_adaptor_closure<simdify_t<Lane>> {
    constexpr auto operator()(auto &&r) const noexcept {
        auto v = std::forward<decltype(r)>(r) | std::views::all;
        auto new_size = std::ranges::size(v) / Lane;
        return std::views::transform(
            std::views::iota(size_t{0}, new_size),
            [v](size_t i) -> decltype(auto) {
                return v[i * Lane];
            }
        );
    }
};
template <auto i>
constexpr simdify_t<i> simdify;

float sum_avx2(std::ranges::range auto &&rng) {
    constexpr auto lane = sizeof(__m256) / sizeof(float);
    __m256 partial_sum {};

    auto simd_view = rng
                   | simdify<lane>;
    for(auto &&simd_v : simd_view) {
        auto addr = &simd_v;
        auto loadu = _mm256_loadu_ps(addr);
        partial_sum = _mm256_add_ps(partial_sum, loadu);
    }

    auto scalar_view = rng
                     | std::views::drop(lane * std::ranges::size(simd_view));
    float sum = std::ranges::fold_left(scalar_view, 0.0f, std::plus());

    float temp[lane];
    _mm256_storeu_ps(std::ranges::data(temp), partial_sum);
    sum = std::ranges::fold_left(temp, sum, std::plus());

    return sum;
}

template <size_t ILP = 4>
float sum_avx2_ilp(std::ranges::range auto &&rng) {
    constexpr auto lane = sizeof(__m256) / sizeof(float);
    constexpr auto bulk = lane * ILP;
    __m256 partial_sum[ILP] {};
    auto process_simd = [&](auto konstexpr_for, auto simd_view) {
        for(auto &&simd_v : simd_view) {
            konstexpr_for([&, addr = &simd_v]<auto Index> {
                auto &partial = partial_sum[Index];
                auto loadu = _mm256_loadu_ps(addr + Index * lane);
                partial = _mm256_add_ps(partial, loadu);
            });
        }
    };

    // SIMD+ILP
    auto bulk_simd_view = rng
                        | simdify<lane>
                        | simdify<ILP>;
    process_simd(constexpr_for<0, ILP>, bulk_simd_view);

    // SIMD
    auto single_simd_view = rng
                          | std::views::drop(bulk * std::ranges::size(bulk_simd_view))
                          | simdify<lane>;
    process_simd(constexpr_for<0, 1>, single_simd_view);

    auto scalar_view = rng
                     | std::views::drop(bulk * std::ranges::size(bulk_simd_view))
                     | std::views::drop(lane * std::ranges::size(single_simd_view));
    float sum = std::ranges::fold_left(scalar_view, 0.0f, std::plus());

    constexpr_for<1, ILP>([&]<size_t Index> {
        partial_sum[0] = _mm256_add_ps(partial_sum[0], partial_sum[Index]);
    });
    float temp[lane];
    _mm256_storeu_ps(std::ranges::data(temp), partial_sum[0]);
    sum = std::ranges::fold_left(temp, sum, std::plus());

    return sum;
}

// ----------------------------------------------------------------------------
// Google Benchmark
// ----------------------------------------------------------------------------

template <size_t Size>
const auto& generate_data() {
    static const auto arr = [] {
        std::array<float, Size> result;
        std::mt19937 gen{std::random_device{}()};
        // No overflow.
        std::uniform_real_distribution<float> dist(-10, 10);
        std::ranges::generate(result, [&] { return dist(gen); });
        return result;
    } ();
    return arr;
}

template <size_t Size>
void BM_run(benchmark::State& state, auto &&func) {
    const auto &rng = generate_data<Size>();

    for(auto _ : state) {
        auto res = func(rng);
        benchmark::DoNotOptimize(res);
    }

    state.SetBytesProcessed(int64_t(state.iterations()) * int64_t(Size) * sizeof(float));
}

template <size_t ...Is>
void register_tests(std::integer_sequence<size_t, Is...>,
                    auto name, auto func) {
    (benchmark::RegisterBenchmark(
        std::string(name) + "/" + std::to_string(Is),
        [func](benchmark::State& state) {
            BM_run<Is>(state, func);
        }
    ), ...);
}

int main(int argc, char** argv) {
    std::integer_sequence<size_t,
        35,
        350,
        3502,
        35023,
        350234
    > seq;

    auto register_test = [seq](auto name, auto func) {
        register_tests(seq, name, func);
    };

    register_test("BM_sum_avx2",
        [](const auto& r) { return sum_avx2(r); });

    register_test("BM_sum_avx2_ilp<4>",
        [](const auto& r) { return sum_avx2_ilp<4>(r); });

    register_test("BM_sum_std_fold_left",
        [](const auto& r) { return std::ranges::fold_left(r, 0.0f, std::plus()); });

    benchmark::Initialize(&argc, argv);
    benchmark::RunSpecifiedBenchmarks();
    benchmark::Shutdown();
    return 0;
}

NOTE: 如果改成延迟更高的浮点指令 benchmark 的话，显式 ILP 还是有价值的。

总结

基本上没别的事情了，因为向量求和算法很固定，所以自动向量化本来挺好的。

但是如果考虑手写/混用 SIMD intrinsic 的话，开发者还是要有手段去控制编译器的优化行为。constexpr_for 是一个极好的控制 ILP 的工具，它是现代且可移植的写法。

处理器的硬件调度能力也是个考虑点，Zen 3 的表现就是优化（暴力展开）到位了，显式 ILP 无所谓。但是前面在你不知道性能瓶颈的时候，简单的显式 ILP 设计就能把性能提升数倍，还不需要做剖析（本文有基线，省了 perf 环节）。考虑到其他更加复杂的算法和延迟更高的向量指令，显式 ILP 的做法仍然是非常的甜点。

std::ranges 的组合玩法很好，但是不看汇编的话你不知道有什么坑。本文中的 ILP 程序性能瓶颈就是 ranges 导致的。