如何在C++中高效地计算2的分数次幂？

Question

我想高效地计算一个有理数的二次幂。这将成为我所有其他数学函数实现的基础（如log、log2、exp、sin、cos），因为我使用牛顿法和大量的迭代方案涉及指数运算。

我想要比cmath库更快的实现方法。由于我们使用的是二进制计算机，与数字2相关的许多操作可以进行位操作。

单精度浮点格式使用32位表示小数部分，第一位是符号位，接下来的8位是带有偏置值127的指数，接下来的24位编码小数部分。双精度浮点格式中，指数用11位表示，偏置值为1023，尾数为52位。

我的想法很简单，我使用乘法的幂规则，即2的a次方乘以2的b次方等于2的(a + b)次方。对于分数指数e，其小数部分d为e % 1，然后其整数部分为i = e - d。由于浮点数的编码方式，我们可以通过组合指数和尾数两部分来组装最终的浮点表示。i加上偏置值直接成为指数部分。

尾数部分有点棘手，但有一些技巧可以解决。尾数被编码为一个带有隐含整数部分为1的二进制小数，我们只需将1加到d上，就可以从浮点位得到二进制小数。注意，如果d是负数，则应该加到2而不是1。

接下来的步骤很简单，第一个位对应于2的0.5次方，第二个位对应于2的0.25次方，第三个位对应于2的0.125次方，依此类推。重要的是它们是常数，因此可以存储在查找表中。因此，对于分数幂，找到尾数的设置位，并将一个初始化为1的数字与查找表中的相应值相乘。

这给出了最终值的尾数。减去1，并将差乘以1 << 23（对于float）或1 << 52（对于double）。结果由位转换组装而成（e << p）| m。对于float，最坏情况下需要进行23次乘法，对于double，需要进行52次乘法。

我的代码：

#include <array>
#include <bitset>
#include <chrono>
#include <cmath>
#include <iostream>
#include <utility>
#include <vector>

using std::array;
using std::chrono::steady_clock;
using std::chrono::duration;
using std::cout;
using std::vector;
float r = 0.0;

const array<double, 52> EXP2_BITS = [] {
    array<double, 52> exp2a;
    for (int i = 1; i < 53; i++) {
        exp2a[i - 1] = exp2(1.0 / (int64_t(1) << i));
    }
    return exp2a;
}();
constexpr int F52 = 1 << 23;
inline float fast_exp2_f(float f) {
    float frac = fmodf(f, 1.0f);
    int i, e;
    e = f - frac + 127;
    frac += frac >= 0 ? 1 : 2;
    uint32_t bits = std::bit_cast<uint32_t>(frac);
    uint32_t m = bits & 0x7fffff;
    i = 22;
    double p = 1.0;
    while (m) {
        if (m & 1) {
            p *= EXP2_BITS[i];
        }
        m >>= 1;
        i--;
    }
    uint32_t m1 = (p - 1.0) * F52;
    return std::bit_cast<float>((e << 23) | m1);
}

template <size_t N, size_t... I>
inline double exp2_helper_impl(std::bitset<N> m, std::index_sequence<I...>) {
    double p = 1.0;
    ((p *= std::get<I>(m)), ...);
    return p;
}

template <std::size_t N>
inline double exp2_helper(const std::bitset<N>& m) {
    return exp2_helper_impl(m, std::make_index_sequence<N>{});
}

inline float fast_exp2_f1(float f) {
    float frac = fmodf(f, 1.0f);
    int i, e;
    e = f - frac + 127;
    frac += frac >= 0 ? 1 : 2;
    uint32_t bits = std::bit_cast<uint32_t>(frac);
    std::bitset<23> m(bits & 0x7fffff);
    double p = exp2_helper(m);
    uint32_t m1 = (p - 1.0) * F52;
    return std::bit_cast<float>((e << 23) | m1);
}

int main()
{
    cout << std::setprecision(7);
    vector<float> bases(256);
    vector<int> ns(256);
    float r256 = 1.0 / 256;
    for (int i = 0; i < 256; i++) {
        bases[i] = 1.0 + rand() % 120 + (rand() % 256) * r256;
        ns[i] = 2 + rand() % 30;
    }
    auto start = steady_clock::now();
    for (int64_t i = 0; i < 1048576; i++) {
        r += fast_exp2_f(bases[i % 256]);
    }
    auto end = steady_clock::now();
    duration<double, std::nano> time = end - start;
    cout << "fast_exp2_f: " << time.count() / 1048576 << " nanoseconds\n";
    start = steady_clock::now();
    for (int64_t i = 0; i < 1048576; i++) {
        r += exp2f(bases[i % 256]);
    }
    end = steady_clock::now();
    time = end - start;
    cout << "exp2f: " << time.count() / 1048576 << " nanoseconds\n";
        start = steady_clock::now();
    for (int64_t i = 0; i < 1048576; i++) {
        r += fast_exp2_f1(bases[i % 256]);
    }
    end = steady_clock::now();
    time = end - start;
    cout << "fast_exp2_f1: " << time.count() / 1048576 << " nanoseconds\n";
    float n;
    for (int i = 0; i < 256; i++) {
        n = bases[i];
        cout << "n: " << n << ", fast_exp2_f: " << fast_exp2_f(n) << ", fast_exp2_f1: " << fast_exp2_f1(n) << ", exp2f: " << exp2f(n) << '\n';
    }
}

我的方法很有效且非常准确，但速度较慢。

fast_exp2_f: 57.74498 nanoseconds
exp2f: 36.58924 nanoseconds
n: 42.13672, fast_exp2_f: 4.83522e+12, exp2f: 4.83522e+12
n: 101.8789, fast_exp2_f: 4.66237e+30, exp2f: 4.66237e+30
n: 79.67969, fast_exp2_f: 9.682243e+23, exp2f: 9.682243e+23
n: 105.2852, fast_exp2_f: 4.942994e+31, exp2f: 4.942994e+31
n: 2.730469, fast_exp2_f: 6.636712, exp2f: 6.636713
n: 12.69922, fast_exp2_f: 6650.369, exp2f: 6650.369
n: 28.23438, fast_exp2_f: 3.157867e+08, exp2f: 3.157867e+08
n: 85.24219, fast_exp2_f: 4.575677e+25, exp2f: 4.575677e+25
n: 53.36719, fast_exp2_f: 1.161781e+16, exp2f: 1.161781e+16
n: 117.0234, fast_exp2_f: 1.688748e+35, exp2f: 1.688748e+35
n: 48.86719, fast_exp2_f: 5.134394e+14, exp2f: 5.134394e+14
n: 19.30078, fast_exp2_f: 645823.8, exp2f: 645823.9
n: 108.7305, fast_exp2_f: 5.384341e+32, exp2f: 5.384341e+32
n: 55.12109, fast_exp2_f: 3.918344e+16, exp2f: 3.918345e+16
n: 3.488281, fast_exp2_f: 11.22218, exp2f: 11.22218
n: 105.1445, fast_exp2_f: 4.483919e+31, exp2f: 4.483919e+31
n: 54.82812, fast_exp2_f: 3.198234e+16, exp2f: 3.198234e+16
n: 45.58594, fast_exp2_f: 5.281223e+13, exp2f: 5.281224e+13
n: 118.2305, fast_exp2_f: 3.898679e+35, exp2f: 3.898679e+35
n: 22.53516, fast_exp2_f: 6077962, exp2f: 6077962
n: 77.85547, fast_exp2_f: 2.734207e+23, exp2f: 2.734207e+23
n: 43.125, fast_exp2_f: 9.592207e+12, exp2f: 9.592208e+12
n: 41.92969, fast_exp2_f: 4.188839e+12, exp2f: 4.188839e+12
n: 89.21094, fast_exp2_f: 7.164207e+26, exp2f: 7.164207e+26
n: 51.28516, fast_exp2_f: 2.743913e+15, exp2f: 2.743913e+15
n: 111.6172, fast_exp2_f: 3.982185e+33, exp2f: 3.982185e+33
n: 34.85938, fast_exp2_f: 3.116861e+10, exp2f: 3.116861e+10
n: 24.07812, fast_exp2_f: 1.771079e+07, exp2f: 1.771079e+07
n: 37.25, fast_exp2_f: 1.634434e+11, exp2f: 1.634434e+11
n: 57.41797, fast_exp2_f: 1.925444e+17, exp2f: 1.925444e+17
n: 25.71484, fast_exp2_f: 5.507307e+07, exp2f: 5.507307e+07
n: 44.62891, fast_exp2_f: 2.720442e+13, exp2f: 2.720442e+13
n: 39.13281, fast_exp2_f: 6.027682e+11, exp2f: 6.027683e+11
n: 102.8789, fast_exp2_f: 9.324739e+30, exp2f: 9.324739e+30
n: 80.85156, fast_exp2_f: 2.181451e+24, exp2f: 2.181451e+24
n: 91.59766, fast_exp2_f: 3.746641e+27, exp2f: 3.746641e+27
n: 114.4453, fast_exp2_f: 2.827951e+34, exp2f: 2.827951e+34
n: 66.17188, fast_exp2_f: 8.312262e+19, exp2f: 8.312262e+19
n: 31.76953, fast_exp2_f: 3.660849e+09, exp2f: 3.660849e+09
n: 94.91016, fast_exp2_f: 3.722236e+28, exp2f: 3.722236e+28
n: 107.9141, fast_exp2_f: 3.057523e+32, exp2f: 3.057523e+32
n: 37.32422, fast_exp2_f: 1.720717e+11, exp2f: 1.720717e+11
n: 16.83594, fast_exp2_f: 116982.8, exp2f: 116982.9

所以我想要进行优化。我认为性能瓶颈在于while循环，我看到使用内联的索引序列来展开循环可以使代码在不到2纳秒的时间内运行，而使用while循环将整数转换为二进制是一种非常低效的方式，因为整数已经以二进制形式存储在内部，所以我们应该直接访问二进制。

我尝试使用位集和索引序列，但没有成功。结果导致编译失败。错误信息：

D:\MyScript\CodeBlocks\testapp\main.cpp:103:23: error: no matching function for call to 'get<0>(std::bitset<23>&)'

我的第一种方法有效，我使用以下方式进行编译：

g++.exe -Wall -fexceptions -fomit-frame-pointer -fexpensive-optimizations -flto -O3 -m64 --std=c++20 -march=native -ffast-math  -c D:\MyScript\CodeBlocks\testapp\main.cpp -o obj\Release\main.o
g++.exe  -o bin\Release\testapp.exe obj\Release\main.o  -O3 -flto -s -static-libstdc++ -static-libgcc -static -m64  

如何优化这个？

---

就算不值得一提，我已经手动展开了循环，代码确实变快了，正如预期的那样，但是它很丑陋，也不够高效。

inline float fast_exp2_f1(float f) {
    int e = int(f);
    float frac = f - e;
    e += 127;
    frac += frac >= 0 ? 1 : 2;
    uint32_t bits = std::bit_cast<uint32_t>(frac);
    std::bitset<23> m(bits & 0x7fffff);
    double p = 1.0;
    if (m[0]) {p *= EXP2_BITS[22]; }
    if (m[1]) {p *= EXP2_BITS[21]; }
    if (m[2]) {p *= EXP2_BITS[20]; }
    if (m[3]) {p *= EXP2_BITS[19]; }
    if (m[4]) {p *= EXP2_BITS[18]; }
    if (m[5]) {p *= EXP2_BITS[17]; }
    if (m[6]) {p *= EXP2_BITS[16]; }
    if (m[7]) {p *= EXP2_BITS[15]; }
    if (m[8]) {p *= EXP2_BITS[14]; }
    if (m[9]) {p *= EXP2_BITS[13]; }
    if (m[10]) {p *= EXP2_BITS[12]; }
    if (m[11]) {p *= EXP2_BITS[11]; }
    if (m[12]) {p *= EXP2_BITS[10]; }
    if (m[13]) {p *= EXP2_BITS[9]; }
    if (m[14]) {p *= EXP2_BITS[8]; }
    if (m[15]) {p *= EXP2_BITS[7]; }
    if (m[16]) {p *= EXP2_BITS[6]; }
    if (m[17]) {p *= EXP2_BITS[5]; }
    if (m[18]) {p *= EXP2_BITS[4]; }
    if (m[19]) {p *= EXP2_BITS[3]; }
    if (m[20]) {p *= EXP2_BITS[2]; }
    if (m[21]) {p *= EXP2_BITS[1]; }
    if (m[22]) {p *= EXP2_BITS[0]; }
    uint32_t m1 = (p - 1.0) * F23;
    return std::bit_cast<float>((e << 23) | m1);
}

fast_exp2_f: 50.11549 nanoseconds
exp2f: 36.29322 nanoseconds
fast_exp2_f1: 11.36074 nanoseconds

---

我现在已经将我的代码转换为无分支形式了。

const array<double, 104> EXP2_BITS = [] {
    array<double, 104> exp2a;
    for (int i = 1; i < 53; i++) {
        exp2a[i * 2 - 2] = 1.0;
        exp2a[i * 2 - 1] = exp2(1.0 / (int64_t(1) << i));
    }
    return exp2a;
}();

constexpr int F23 = 1 << 23;
inline float fast_exp2_f1(float f) {
    int e = int(f);
    float frac = f - e;
    e += 127;
    frac += frac >= 0 ? 1 : 2;
    uint32_t bits = std::bit_cast<uint32_t>(frac);
    std::bitset<23> m(bits & 0x7fffff);
    double p = 1.0;
    p *= EXP2_BITS[44+m[0]];
    p *= EXP2_BITS[42+m[1]];
    p *= EXP2_BITS[40+m[2]];
    p *= EXP2_BITS[38+m[3]];
    p *= EXP2_BITS[36+m[4]];
    p *= EXP2_BITS[34+m[5]];
    p *= EXP2_BITS[32+m[6]];
    p *= EXP2_BITS[30+m[7]];
    p *= EXP2_BITS[28+m[8]];
    p *= EXP2_BITS[26+m[9]];
    p *= EXP2_BITS[24+m[10]];
    p *= EXP2_BITS[22+m[11]];
    p *= EXP2_BITS[20+m[12]];
    p *= EXP2_BITS[18+m[13]];
    p *= EXP2_BITS[16+m[14]];
    p *= EXP2_BITS[14+m[15]];
    p *= EXP2_BITS[12+m[16]];
    p *= EXP2_BITS[10+m[17]];
    p *= EXP2_BITS[8+m[18]];
    p *= EXP2_BITS[6+m[19]];
    p *= EXP2_BITS[4+m[20]];
    p *= EXP2_BITS[2+m[21]];
    p *= EXP2_BITS[0+m[22]];
    uint32_t m1 = (p - 1.0) * F23;
    return std::bit_cast<float>((e << 23) | m1);
}

它应该更稳定，条件开销被消除了。但是每个位都有一个乘法赋值开销，所以代码执行时间较长。

fast_exp2_f1: 20.96491 nanoseconds

---

我做到了。我彻底击败了exp2f。

const array<double, 52> EXP2_BITS = [] {
    array<double, 52> exp2a;
    for (int i = 1; i < 53; i++) {
        exp2a[i - 1] = exp2(1.0 / (int64_t(1) << i));
    }
    return exp2a;
}();

constexpr int F23 = 1 << 23;
inline float fast_exp2_f(float f) {
    int e = int(f);
    float frac = f - e;
    e += 127;
    frac += frac >= 0 ? 1 : 2;
    uint32_t bits = std::bit_cast<uint32_t>(frac);
    std::bitset<23> m(bits & 0x7fffff);
    double p = 1.0;
    p *= m[0] ? EXP2_BITS[22] : 1;
    p *= m[1] ? EXP2_BITS[21] : 1;
    p *= m[2] ? EXP2_BITS[20] : 1;
    p *= m[3] ? EXP2_BITS[19] : 1;
    p *= m[4] ? EXP2_BITS[18] : 1;
    p *= m[5] ? EXP2_BITS[17] : 1;
    p *= m[6] ? EXP2_BITS[16] : 1;
    p *= m[7] ? EXP2_BITS[15] : 1;
    p *= m[8] ? EXP2_BITS[14] : 1;
    p *= m[9] ? EXP2_BITS[13] : 1;
    p *= m[10] ? EXP2_BITS[12] : 1;
    p *= m[11] ? EXP2_BITS[11] : 1;
    p *= m[12] ? EXP2_BITS[10] : 1;
    p *= m[13] ? EXP2_BITS[9] : 1;
    p *= m[14] ? EXP2_BITS[8] : 1;
    p *= m[15] ? EXP2_BITS[7] : 1;
    p *= m[16] ? EXP2_BITS[6] : 1;
    p *= m[17] ? EXP2_BITS[5] : 1;
    p *= m[18] ? EXP2_BITS[4] : 1;
    p *= m[19] ? EXP2_BITS[3] : 1;
    p *= m[20] ? EXP2_BITS[2] : 1;
    p *= m[21] ? EXP2_BITS[1] : 1;
    p *= m[22] ? EXP2_BITS[0] : 1;
    uint32_t m1 = (p - 1.0) * F23;
    return std::bit_cast<float>((e << 23) | m1);
}

即使数组中有32768个随机值，上述方法仍然更快。

fast_exp2_f: 11.4665 nanoseconds
exp2f: 36.11965 nanoseconds

但仍然存在条件开销。

- Ξένη Γήινος

---

你可能想要更详细地说明“将二提升为有理数的幂”。所有的数值二进制浮点操作数都是有理数（分数），所以你可能是在询问类似于exp2f()的函数，这是一个在ISO-C99中提供的示例实现，可以在这里找到。 - undefined

这个问题似乎也与这个问题有关。 - undefined

这个问题应该明确期望的精度级别。您是否需要完美的精度，直到IEEE浮点数的最后一位数字？您是否需要考虑浮点环境中的舍入模式？您愿意依赖于哪个级别的CPU扩展？如果有的话，是否可以使用cpuid进行运行时切换？您是否希望完全避免内联汇编和/或内部函数？有这么多需要回答的问题... - undefined

有一条单独的FPU（x87）指令可以实现这个功能（减去1）。我怀疑在x86 CPU上很难找到比这更准确和高效的方法。请参考https://en.wikibooks.org/wiki/X86_Assembly/Floating_Point上的f2xm1。 - undefined

如果你能在这方面击败FPU，我会感到惊讶。这个游戏真的值得吗？ - undefined

Answer 1

TL;DR：实现受GCC上的依赖条件跳转链限制。Clang生成更好的汇编代码（得益于SIMD指令和无分支代码）。基准测试明显存在偏见。与SIMD相比，标量exp2f实现效率低下，因此如果可能的话应该使用后者。

---

分析

fast_exp2_f1的主要问题是长的依赖条件指令链。实际上，这是GCC生成的代码的一部分：

[...]
.L33:
        test    ah, 16
        je      .L34
        vmulsd  xmm0, xmm0, QWORD PTR .LC21[rip]
.L34:
        test    ah, 32
        je      .L35
        vmulsd  xmm0, xmm0, QWORD PTR .LC22[rip]
.L35:
        test    ah, 64
        je      .L36
        vmulsd  xmm0, xmm0, QWORD PTR .LC23[rip]
.L36:
        test    ah, -128
        je      .L37
        vmulsd  xmm0, xmm0, QWORD PTR .LC24[rip]
.L37:
        test    eax, 65536
        je      .L38
        vmulsd  xmm0, xmm0, QWORD PTR .LC25[rip]
.L38:
        test    eax, 131072
        je      .L39
        vmulsd  xmm0, xmm0, QWORD PTR .LC26[rip]
.L39:
        test    eax, 262144
        je      .L40
        vmulsd  xmm0, xmm0, QWORD PTR .LC27[rip]
[...]

处理器需要评估条件跳转是否被执行，并根据结果进行跳转。当条件跳转被执行时，处理器无法在1个周期内完成。问题是，你有23个条件跳转，如果全部执行，最快也需要23个周期。在我的i5-9600KF处理器上，这意味着至少需要5.11纳秒来处理这些条件。如果条件跳转未执行，则FMA延迟成为瓶颈，因为在大多数主流机器上，双精度乘法通常需要2-8个周期。在我的i5-9600KF处理器上，需要4个周期。因此，如果所有条件跳转都未执行，代码至少需要23*4=92个周期，因为所有乘法操作都是依赖的。请注意，处理器无法重新排序它们，因为双精度数的非结合性，虽然编译器理论上可以重新排序，但由于条件的复杂性，这对它们来说太复杂了。实际上，在基准测试中，条件跳转的执行情况是混合的，有些被执行，有些未执行，所以只有当多个条件跳转未执行时，FMA延迟才会成为问题。事实上，平均而言，有4个条件跳转未执行，其他都执行了。这意味着任何一个条件跳转有83%的概率被执行。

条件跳转的成本是可变的，很难预测。当进行条件跳转时，需要首先从内存中的不连续位置获取和解码要执行的指令。问题在于结果取决于流水线中其他部分产生的数据。这可能会导致显著的开销，因为流水线停顿时间较长（通常超过10个周期）。现代主流处理器非常复杂，经过几十年的智慧工程师和研究人员的设计，他们试图尽量减少这种开销。主流处理器现在可以预测条件跳转，以减少流水线停顿时间。当条件预测良好时，在最新的主流处理器上通常（几乎）没有流水线停顿。当出现错误预测时，开销可能会很大，尤其是由于条件跳转链的长度。

事实证明，我的处理器（像大多数最近的主流处理器一样）能够预测许多条件跳转。实际上，它能够根据正在计算的实际数字来判断跳转是否发生。嗯，这是一个巨大的简化，因为实际情况要复杂得多（分支预测是一个非常复杂的主题）。事实上，输入数组中只有256个双精度数字，所以处理器可以记住大多数条件跳转的结果，以便成功预测它们（顺便说一句，这相当疯狂）。简而言之，这意味着基准测试由于处理器的优化而存在偏差。通过将输入数组中的项目数量从256增加到32768（在Windows上使用Clang进行快速检查，不使用`-ffast-math`生成接近GCC的代码），可以在我的机器上看到这一点。

256 items:
 - fast_exp2_f: 44.950962 nanoseconds
 - exp2f: 40.409470 nanoseconds
 - fast_exp2_f1: 13.479328 nanoseconds

32768 items:
 - fast_exp2_f: 68.788147 nanoseconds
 - exp2f: 40.300083 nanoseconds
 - fast_exp2_f1: 45.244217 nanoseconds

可以看到，exp2f 函数不受输入数组中项目数量的影响。这是因为它内部要么没有条件跳转，要么只有很少的条件跳转，并且这些跳转很容易预测（即不太依赖于基准测试中实际使用的输入数字）。事实上，在更实际的使用情况下，fast_exp2_f1 实际上比 exp2f 函数更慢。实际上，输入数字往往彼此接近，所以 fast_exp2_f1 相对于 exp2f 有优势。如果使用范围更广的随机数作为输入，fast_exp2_f1 的速度会更慢。在我的机器上，它应该需要比 60 纳秒更长的时间。

---

如何加快代码速度

一般来说，避免分支预测错误的简单解决方案就是编写一个无分支的代码。在x86-64处理器上，由于SIMD（SSE）指令如blendpd的存在，这是可能的。然而，在这里这并不是一个好的解决方案，因为大多数条件跳转都会被执行，所以大多数乘法都不会被执行。使用一个简单的无分支代码会导致FMA延迟成为一个问题（至少在我的机器上，这部分代码的延迟为20纳秒）。

一个更聪明的解决方案是将乘法进行二进制规约并与混合操作结合起来。自己手动实现这个过程是繁琐、不可移植且低级的。事实上，GCC倾向于不生成无分支代码，据我所知，没有办法通过可移植的代码来强制它生成无分支代码。因此，SIMD（SSE）内联函数是进行这种优化的方法。实际上，人们可以使用SIMD指令同时操作多个数字。事实证明，Clang在使用-O3 -ffast-math时做得非常好。以下是Clang生成的代码的一部分：

        vpand   ymm6, ymm0, ymmword ptr [rip + .LCPI1_5]
        vpand   ymm7, ymm0, ymmword ptr [rip + .LCPI1_6]
        vpand   ymm8, ymm0, ymmword ptr [rip + .LCPI1_7]
        vpand   ymm9, ymm0, ymmword ptr [rip + .LCPI1_8]
        vpcmpeqq        ymm9, ymm9, ymm3
        vpcmpeqq        ymm7, ymm7, ymm3
        vpcmpeqq        ymm6, ymm6, ymm3
        vmovupd xmm10, xmmword ptr [rip + EXP2_BITS+168]
        vunpcklpd       xmm5, xmm10, xmm5       # xmm5 = xmm10[0],xmm5[0]
        vpermpd ymm10, ymmword ptr [rip + EXP2_BITS+152], 20 # ymm10 = mem[0,1,1,0]
        vblendpd        ymm5, ymm5, ymm10, 12           # ymm5 = ymm5[0,1],ymm10[2,3]
        vblendvpd       ymm5, ymm5, ymmword ptr [rip + .LCPI1_9], ymm7
        vpermpd ymm7, ymmword ptr [rip + EXP2_BITS+120], 27 # ymm7 = mem[3,2,1,0]
        vpermpd ymm10, ymmword ptr [rip + EXP2_BITS+88], 27 # ymm10 = mem[3,2,1,0]
        vblendvpd       ymm6, ymm10, ymm4, ymm6
        vpermpd ymm10, ymmword ptr [rip + EXP2_BITS+56], 27 # ymm10 = mem[3,2,1,0]
        vpcmpeqq        ymm3, ymm8, ymm3
        vblendvpd       ymm3, ymm10, ymm4, ymm3
        vblendvpd       ymm4, ymm7, ymm4, ymm9
        vmulpd  ymm3, ymm4, ymm3
        vmulpd  ymm3, ymm6, ymm3
        vmulpd  ymm3, ymm5, ymm3
        vextractf128    xmm4, ymm3, 1
        vmulpd  xmm3, xmm3, xmm4
        vshufpd xmm4, xmm3, xmm3, 1             # xmm4 = xmm3[1,0]
        vmulsd  xmm3, xmm3, xmm4
        vextractf128    xmm4, ymm2, 1
        vmulpd  xmm2, xmm2, xmm4
        vshufpd xmm4, xmm2, xmm2, 1             # xmm4 = xmm2[1,0]
        vmulsd  xmm2, xmm2, xmm4
        vpand   xmm0, xmm0, xmmword ptr [rip + .LCPI1_10]
        vxorpd  xmm4, xmm4, xmm4
        vpcmpeqq        xmm0, xmm0, xmm4
        vmovupd xmm4, xmmword ptr [rip + EXP2_BITS+8]
        vblendvpd       xmm0, xmm4, xmmword ptr [rip + .LCPI1_11], xmm0

可以注意到，vmulpd一次操作多个值，vpand，vpcmpeqq和vblendvpd以无分支的方式执行条件。这段代码相对快速，并且与输入值无关。

256 items:
 - fast_exp2_f: 43.608665 nanoseconds
 - exp2f: 40.797138 nanoseconds
 - fast_exp2_f1: 15.594292 nanoseconds

32768 items:
 - fast_exp2_f: 70.216560 nanoseconds
 - exp2f: 41.187191 nanoseconds
 - fast_exp2_f1: 16.847897 nanoseconds

可以看到，fast_exp2_f1的时间不受输入数组大小的影响。事实上，fast_exp2_f1比exp2f更快。

简而言之：如果可以的话，请在这里使用Clang。

---

注意事项

在这种情况下，所有的实现都是低效的。你可以使用SIMD指令来同时处理多个输入值。这比Clang对fast_exp2_f1使用-ffast-math更高效（即对函数进行向量化，但仍然一次处理一个输入值）。大多数优化的数学库都这样做（如Intel SVML、ISPC，甚至是索尼PS4/PS5上的标准数学库）。我预计这些优化库在我的机器上的速度比每个项目10纳秒要快得多。然而，它们通常会引入几个ULP的小误差。在你的情况下，这不应该是一个问题，因为你无论如何都使用了-ffast-math。

只是为了好玩，我尝试在你的基准测试中使用ISPC。看起来ISPC没有实现exp2f，但有pow，可能比直接计算exp2f要慢（因为pow通常计算一个log和一个exp，但我不确定基数是否为2）。如预期的那样，ISPC比你的基准测试中的所有实现都要快得多（精度达到>=5位数）。在我的机器上，它比其他实现要快30倍以上（得益于AVX-2 SIMD指令集和积极优化的无分支计算）。

ISPC pow(2.0f,x): 0.417328 nanoseconds

我建议你不要重复造轮子，如果可能的话，使用SIMD实现，并阅读相关的研究论文。实现比优化的数学库更快的指数函数，同时提供至少相同的精度，对于即使是高技能的开发者来说也非常困难。

- Jérôme Richard