如何在 C++ 中有效地将 2 提高到小数次幂？-解网

问：

我想有效地将两个提升到理性的力量。这将是我所有其他数学函数实现（）的基础，因为我使用牛顿方法，并且迭代方案的负载涉及幂。log, log2, exp, sin, cos

我想实现比 cmath 更快的方法。因为我们使用二进制计算机，所以很多与数字二相关的东西都可能被黑客入侵。

单精度浮点格式使用 32 位表示分数，第一位是符号位，接下来的 8 位是偏置为 127 的指数，接下来的 24 位对分数进行编码。对于双精度，指数编码为 11 位，偏差为 1023，尾数为 52 位。

我的想法很简单，我使用乘法的幂法则，在这种情况下是 2 a * 2 b = 2^{a + b}。对于分数指数 e，其小数部分 d 为 e % 1，然后其积分部分是 i = e - d。因为浮点的编码方式，我们可以通过组合指数和尾数两部分来组装最终的浮点表示。我带着它的偏差直接进入指数部分。

尾数部分有点棘手，但有一些技巧。尾数被编码为二进制小数，隐式整数部分是 1，我们只需将 1 添加到 d 即可从浮点位中获取二进制分数。请注意，如果 d 为负数，则应将其添加到 2 而不是 1。

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

这给出了最终值的尾数。减去 1 并将差值乘以 1 << 23 表示浮点数，1 << 52 表示双倍数。结果由钻头铸造组装而成（e << p） |对于浮点数，最坏的情况是 23 次乘法，对于双精度，最坏情况是 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

但仍然存在条件开销。

C++ 算法性能浮点 C++20

2赞 njuffa 9/27/2023

您可能希望更详细地指定“将 2 提高到有理幂”。所有数值二进制浮点操作数都是有理数（分数），因此您可能会询问类似的东西，此处提供了 ISO-C99 中的示例实现。exp2f()

1赞 Bob__ 9/27/2023

这似乎也相关。

2赞 Jan Schultke 9/27/2023

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

1赞 Simon Goater 9/27/2023

有一个 FPU （x87）指令可以执行此操作（减 1）。我怀疑在 x86 CPU 上的准确性和性能很难超越它。请参阅 en.wikibooks.org/wiki/X86_Assembly/Floating_Point 上的 f2xm1

1赞 Jan Schultke 9/27/2023

还有一个 AVX-512 指令 felixcloutier.com/x86/vexp2ps 它计算 2^X 而不使用 -1。glibc 实现已经使用了 x87 指令，因此无论您编写什么代码，都不太可能击败它。如果 glibc 方法不够准确，您可能应该使用牛顿方法或其他方法对其进行改进，使其稍微准确一些，否则请使用标准库。

答：

6赞 Jérôme Richard 9/28/2023 #1

TL;DR：实现受 GCC 上依赖条件跳转链的约束。Clang 生成更好的汇编代码（感谢 SIMD 指令和无分支代码）。该基准显然是有偏见的。与 SIMD 实现相比，标量实现效率低下，因此如果可能，应使用后者。exp2f

分析

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

[...]
.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 ns 仅适用于一系列条件。如果不进行跳转，则 FMA 延迟是瓶颈，因为在大多数主流机器上，双精度乘法通常需要大约 2-8 个周期。在我的 i5-9600KF 处理器上，它需要 4 个周期。因此，如果未执行所有条件跳转，则代码至少需要周期，因为所有乘法运算都是相互依赖的。请注意，由于双精度数字的非关联性，处理器无法对它们进行重新排序，虽然编译器理论上可以对它们进行重新排序，但由于条件的原因，这对他们来说太复杂了。在实践中，有条件跳跃和未有条件跳跃的混合，因此，只有在未进行多次条件跳跃时，FMA 延迟才是一个问题。事实证明，大多数条件跳跃都是在基准测试中进行的，但不是全部。事实上，平均而言，没有进行 4 次有条件的跳跃，而所有其他的都有。这意味着任何有条件的跳跃都有 83% 的机会被采取。23*4=92

条件跳跃的成本是可变的，难以预测。当进行条件跳转时，首先需要从内存中不同的非连续位置获取和解码要执行的指令。问题是结果取决于管道其他部分产生的数据。由于管道停滞（通常为 >10 个周期），这可能会带来显着的开销。现代主流处理器非常复杂，它们是由聪明的工程师和研究人员设计了几十年，他们试图尽可能地减少这种开销。主流处理器现在预测条件跳跃，以减少管道停顿。当一个条件被很好地预测时，在最近的主流处理器上通常（几乎）没有管道停顿。当出现预测失误时，开销可能会很大，尤其是由于条件跳跃的长链。

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

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

可以看出，这不受输入数组中项目数量的影响。这是因为其中要么没有条件跳跃，要么很少有被很好地预测的跳跃（即与基准测试中使用的实际输入数字没有太大关系）。最后，在更现实的用例中，fast_exp2_f1实际上比 exp2f 慢。事实上，输入数字往往彼此接近，因此比具有优势。使用范围更广的输入随机数应该会减慢速度。在我的机器上，它应该花费超过 60 ns 的时间。exp2ffast_exp2_f1exp2ffast_exp2_f1

如何使代码更快

通常，避免分支错误预测问题的简单解决方案是编写无分支代码。这在 x86-64 处理器上是可能的，这要归功于 SIMD （SSE）指令，例如 .然而，这并不是一个好的解决方案，因为大多数条件跳跃都是采用的，因此大多数乘法都不会执行。使用朴素的无分支代码会导致 FMA 延迟成为一个问题（这部分代码在我的机器上至少需要 20 ns）。blendpd

更聪明的解决方案是将乘法与混合相结合进行二进制约简。自己做这件事很乏味，不便携和低级。事实上，GCC 倾向于不生成无分支代码，而 AFAIK 没有办法使用可移植代码强制它。因此，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

可以注意到一次对多个值的操作，以无分支的方式执行条件。此代码相对较快，并且与输入值无关：vmulpdvpandvpcmpeqqvblendvpd

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_f1fast_exp2_f1exp2f

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

笔记

在这种情况下，所有实现都是低效的。您可以使用 SIMD 指令一次对多个输入值进行操作。这比 Clang 所做的使用更有效（即矢量化函数，该函数一次仍对 1 个输入值进行操作）。大多数优化的数学库都可以做到这一点（如英特尔 SVML 或 ISPC，甚至是 PS4/PS5 上索尼的标准数学库）。我希望这种优化的库在我的机器上明显快于 10 ns/项。但是，它们通常会引入很少的 ULP 的小误差。这在您的情况下应该不是问题，因为您无论如何都会使用。fast_exp2_f1-ffast-math-ffast-math

只是为了好玩，我尝试在您的基准测试中使用 ISPC。看起来 ISPC 没有实现，但可能比直接计算慢（因为经常计算 a 和 an，但我不确定常数基数是否为 2）。正如预期的那样，ISPC 比基准测试中的所有实现都要快得多（精度为 >=5 位）。它比我机器上的其他实现快 >30 倍（这要归功于 AVX-2 SIMD 指令集和积极优化的无分支计算）：exp2fpowexp2fpowlogexp

ISPC pow(2.0f,x): 0.417328 nanoseconds

我建议你不要重新发明轮子，使用 SIMD 实现是可能的，并阅读有关它的研究论文。实现比优化的数学库更快的指数函数，同时提供至少相同的精度是非常困难的，即使对于技术娴熟的开发人员也是如此。

上一个：电话号码格式，在国家/地区代码后加点

下一个：使用 %f 但获取“format %f expects type double”

如何在 C++ 中有效地将 2 提高到小数次幂？

How to efficiently raise 2 to a fractional power in C++?

评论

分析

如何使代码更快

笔记