如何在 C++ 中高效应用多项式而无需循环？

Question

我想获得一些复杂函数的准确近似值（pow, exp, log,log2比 C++ 标准库中 cmath 提供的函数更快地

为此，我想利用浮点编码方式并使用位操作获取指数和尾数，然后进行多项式近似。尾数在 1 和 2 之间，因此我使用 n 阶多项式来近似 [1, 2] 中域 x 中的目标函数，并对浮点表达式进行位操作和简单数学运算以使计算有效。

我用来np.polyfit生成多项式。作为示例，以下是我用来在 1 <= x <= 2 上近似 log2 的 7 阶多项式：

P = np.array(
    [
        0.01459855,
        -0.17811046,
        0.95074541,
        -2.91450247,
        5.67353733,
        -7.39616658,
        7.08511059,
        -3.23521156,
    ],
    dtype=float,
)

要应用多项式，请将第一项乘以 x 的 7 次方，第二项乘以 x 的 6 次方，依此类推...

在代码中：

P[0] * x**7 + P[1] * x**6 + P[2] * x**5 + P[3] * x**4 + P[4] * x**3 + P[5] * x**2 + P[6] * x + P[7]

当然，这是非常低效的，首先计算较大的幂，因此存在大量重复计算，如果我们颠倒顺序，我们可以根据之前的幂计算出当前的幂，如下所示：

PR = P[::-1]
s = 0
c = 1
for i in PR:
    s += i * c
    c *= x

这正是我在 C++ 中所做的：

PR = P[::-1]
s = 0
c = 1
for i in PR:
    s += i * c
    c *= x

它比 cmath 的 log2 快得多，同时获得相同的精度：

log2(3.1415927f) = 1.651496171951294 : 42.68856048583984 nanoseconds
fast_log2_accurate(3.1415927f) = 1.651496171951294 : 9.967899322509766 nanoseconds

我使用 Visual Studio 2022 进行编译，编译器标志：

/permissive- /ifcOutput "x64\Release\" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /Ob1 /sdl /Fd"x64\Release\vc143.pdb" /Zc:inline /fp:fast /D "NDEBUG" /D "_CONSOLE" /D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope /std:c17 /Gd /Oi /MD /std:c++20 /FC /Fa"x64\Release\" /EHsc /nologo /Fo"x64\Release\" /Ot /Fp"x64\Release\exponentiation.pch" /diagnostics:column /arch:AVX2

不过我认为这样可以更有效。有一个循环开销，如果我可以优化循环，它应该会更快。

如何在没有循环的情况下应用多项式？

---

如果循环已经展开，那么是否可以使用 SIMD 指令进行计算以使其更快？

---

我对下面提供的解决方案以及我之前编写的其他一些函数进行了基准测试：

#include <vector>
#include <numbers>
using std::vector;

using std::numbers;
using numbers::ln2;
using numbers::pi;

constexpr double LOG2_POLY7[8] = {
    -3.23521156,
    7.08511059,
    -7.39616658,
    5.67353733,
    -2.91450247,
    0.95074541,
    -0.17811046,
    0.01459855,
};
constexpr float FU = 1.0 / (1 << 23);

inline float fast_log2_accurate(float f) {
    uint32_t bits = std::bit_cast<uint32_t>(f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = 0;
    double m = 1;
    for (const double& p : LOG2_POLY7) {
        s += p * m;
        m *= m1;
    }
    return e + s;
}

template <int N> inline double poly(const double* a, const float x) {
    return (a[0] + x * poly<N - 1>(a + 1, x));
}

template <> inline double poly<0>(const double* a, const float x) {
    return x * a[0];
}

inline float fast_log2_accurate2(float f) {
    uint32_t bits = std::bit_cast<uint32_t>(f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    return e + poly<8>(LOG2_POLY7, m1);
}


inline float fast_log2_accurate3(float f) {
    uint32_t bits = std::bit_cast<uint32_t>(f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = 0;
    double m = 1;
    for (int i = 0; i < 8; i++) {
        s += LOG2_POLY7[i] * m;
        m *= m1;
    }
    return e + s;
}

vector<float> log2_float() {
    int lim = 1 << 24;
    vector<float> table(lim);
    for (int i = 0; i < lim; i++) {
        table[i] = float(log(i) / ln2) - 150;
    }
    return table;
}
const vector<float> LOG2_FLOAT = log2_float();

inline float fast_log2(float f) {
    uint32_t bits = std::bit_cast<uint32_t>(f);
    int e = (bits >> 23) & 0xff;
    int m = bits & 0x7fffff;
    return (e == 0 ? LOG2_FLOAT[m << 1] : e + LOG2_FLOAT[m | 0x00800000]);
}

inline float fast_log(float f) {
    return fast_log2(f) * ln2;
}

vector<double> log2_double() {
    int lim = 1 << 24;
    vector<double> table(lim);
    for (uint64_t i = 0; i < lim; i++) {
        table[i] = log(i << 29) / ln2 - 1075;
    }
    return table;
}
const vector<double> LOG2_DOUBLE = log2_double();

inline double fast_log2_double(double d) {
    uint64_t bits = std::bit_cast<uint64_t>(d);
    uint64_t e = (bits >> 52) & 0x7ff;
    uint64_t m = bits & 0xfffffffffffff;
    return (e == 0 ? LOG2_DOUBLE[m >> 28] : e + LOG2_DOUBLE[(m | 0x10000000000000) >> 29]);
}

fast_log2(3.1415927f) = 1.651496887207031 : 0.2610206604003906 nanoseconds
log2f(3.1415927f) = 1.651496171951294 : 33.27693939208984 nanoseconds
fast_log2_double(pi) = 1.651496060131421 : 0.3225326538085938 nanoseconds
fast_log2_accurate(3.1415927f) = 1.651496171951294 : 8.907032012939453 nanoseconds
fast_log2_accurate3(3.1415927f) = 1.651496171951294 : 7.831001281738281 nanoseconds
fast_log2_accurate2(3.1415927f) = 1.651496171951294 : 13.57889175415039 nanoseconds

虽然使用查找表的两个函数是无与伦比的，但它们相当不准确。我log2f在我的基准测试中明确使用了。正如您所看到的，在 MSVC 中它非常慢。

正如预期的那样，递归函数会显着减慢代码速度。使用旧式循环可使代码运行速度加快 2 纳秒。然而，我无法对使用的那个进行基准测试std::index_sequence，它导致编译器错误，并且我无法解决它。

---

我的基准测试代码中有一个错误，导致递归版本计时不准确，它使测量的时间更长，我修复了这个问题。

最新答案的解决方案：

constexpr double LOG2_POLY7[8] = {
    -3.23521156,
    7.08511059,
    -7.39616658,
    5.67353733,
    -2.91450247,
    0.95074541,
    -0.17811046,
    0.01459855,
};
constexpr float FU = 1.0 / (1 << 23);

inline float fast_log2_accurate(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = 0;
    double m = 1;
    for (const double& p : LOG2_POLY7) {
        s += p * m;
        m *= m1;
    }
    return e + s;
}

fast_log2_accurate4(3.1415927f) = 1.651496887207031 : 17.01173782348633 nanoseconds

它不像我的代码那么准确并且需要更长的时间，因为每次迭代都更昂贵。

---

索引序列版本之前无法编译，因为我使用double[8]而不是std::array<double, 8>，我以为它们是同一个东西！经指出后，我修复了该问题，并且编译成功。

基准：

ln(256) = 5.545613288879395 : 3.985881805419922 nanoseconds
log(256) = 5.545177459716797 : 7.047939300537109 nanoseconds
fast_log2(3.1415927f) = 1.651496887207031 : 0.25787353515625 nanoseconds
log2f(3.1415927f) = 1.651496171951294 : 35.03541946411133 nanoseconds
fast_log2_double(pi) = 1.651496060131421 : 0.3331184387207031 nanoseconds
fast_log2_accurate(3.1415927f) = 1.651496171951294 : 9.366512298583984 nanoseconds
fast_log2_accurate3(3.1415927f) = 1.651496171951294 : 7.454872131347656 nanoseconds
fast_log2_accurate2(3.1415927f) = 1.651496171951294 : 14.07079696655273 nanoseconds
fast_log2_accurate4(3.1415927f) = 1.651496887207031 : 16.6351318359375 nanoseconds
fast_log2_accurate5(3.1415927f) = 1.651496171951294 : 7.868862152099609 nanoseconds

事实证明它ln非常快，唯一击败它的方法是使用查找表，但它只给出 3 个正确的十进制数字，np.log(256) 给出 5.545177444479562。相比之下，我最快的函数可以给出 6 个正确的十进制数字，并且速度快十倍。我只需要乘以它就ln2可以得到ln(x)，而且它会更准确。

我对解决方案进行了多次基准测试，并且fast_log2_accurate5是索引序列版本。它和旧式循环版本始终比我基于范围的 for 循环版本更快。有时 for 循环版本更快，有时索引序列版本更快。在这个级别上，测量值波动很大，并且我同时运行许多其他程序。

但看起来索引序列版本的性能比for循环版本稳定得多，所以我会接受它。

---

更新：

我重新审视了代码，我对索引序列版本做了一点更改，我只是inline在函数前面添加了do_loop这一点，这个小小的更改使代码在不到一纳秒的时间内运行，我可以添加更多术语在不降低代码速度太多的情况下，它仍然比log2获得非常准确的结果要快得多。

相比之下，std::apply即使使用内联，该版本也很慢：

constexpr std::array<double, 8> LOG2_POLY7A = {
    -3.23521156,
    7.08511059,
    -7.39616658,
    5.67353733,
    -2.91450247,
    0.95074541,
    -0.17811046,
    0.01459855,
};

template <std::size_t... I>
inline double do_loop(double m1, std::index_sequence<I...>) {
    double s = 0;
    double m = 1;
    ((s += std::get<I>(LOG2_POLY7A) * m, m *= m1), ...);
    return s;
}

inline float fast_log2_accurate5(float f) {
    uint32_t bits = std::bit_cast<uint32_t>(f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = do_loop(m1, std::make_index_sequence<8>{});
    return e + s;
}

inline double do_loop1(double m1) {
    double s = 0;
    double m = 1;
    auto worker = [&](auto&...term) {
        ((s += term * m, m *= m1), ...);
        };
    std::apply(worker, LOG2_POLY7A);
    return s;
}

inline float fast_log2_accurate6(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = do_loop1(m1);
    return e + s;
}

fast_log2_accurate5(3.1415927f) = 1.651496171951294 : 0.9766578674316406 nanoseconds
fast_log2_accurate6(3.1415927f) = 1.651496171951294 : 7.168102264404297 nanoseconds

Answer 1

您可以使用 展开循环std::index_sequence，如下所示：

#include <array>
#include <cstdint>
#include <utility>

constexpr std::size_t size_log2 = 8;

constexpr std::array<double, size_log2> LOG2_POLY7 = {
    -3.23521156,
    7.08511059,
    -7.39616658,
    5.67353733,
    -2.91450247,
    0.95074541,
    -0.17811046,
    0.01459855,
};
constexpr float FU = 1.0 / (1 << 23);

template <std::size_t... I>
double do_loop(double m1, std::index_sequence<I...>) {
    double s = 0;
    double m = 1;
    ((s += std::get<I>(LOG2_POLY7) * m, m *= m1),...);
    return s;
}

inline float fast_log2_accurate(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = do_loop(m1, std::make_index_sequence<size_log2>{});
    return e + s;
}

在godbolt上测试代码。

请注意，我还将变量移到了函数m内部do_loop，因为只有在那里才需要它。并且，根据评论中的建议，返回存储在问题do_loop变量中的结果。s与您的初始版本相比，循环的展开是在编译时完成的，并且避免了p与LOG2_POLY7.end()。与往常一样，应该对实际增益进行基准测试。

另外，正如注释中所指出的，可以通过使用std::apply将数组转换LOG2_POLY7为可变参数来简化循环；请参阅下面或其他答案。

正如评论中所要求的，人们可以概括该do_loop函数以适用于通用std::array. 然而，这意味着额外的例程。这是因为，在do_loop上面，的类型index_sequence是由非类型模板参数决定的std::size_t... I；这些参数，以及第二个参数的类型都do_loop被推导出来。现在，给定一个泛型std::array，您可以推断出大小N，但为了完成do_loop工作，您需要将其转换N为索引序列0...N-1：这正是index_sequence。

因此，对于更通用的例程，可以将上面的代码替换为：

template <std::size_t N, std::size_t... I>
double do_loop_impl(double m1, const std::array<double, N>& data, std::index_sequence<I...>) {
    double s = 0;
    double m = 1;
    ((s += std::get<I>(data) * m, m *= m1),...);
    return s;
}

template <std::size_t N>
double do_loop(double m1, const std::array<double, N>& data) {
    return do_loop_impl(m1, data, std::make_index_sequence<N>{});
}

inline float fast_log2_accurate(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = do_loop(m1, LOG2_POLY7);
    return e + s;
}

或者，可以使用std::apply，稍微概括另一个答案：

template <std::size_t N>
double do_loop(double m1, const std::array<double, N>& data) {
    return std::apply([m1](auto... p) {
                      double s = 0;
                      double m = 1;
                      ((s += p * m, m *= m1), ...);
                      return s; }, data);
}

inline float fast_log2_accurate(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = ((bits >> 23) & 0xff) - 127;
    double m1 = 1 + (bits & 0x7fffff) * FU;
    double s = do_loop(m1, LOG2_POLY7);
    return e + s;
}

include <tuple>（如果您使用，请记住std::apply）。

最后，如果您对速度感兴趣，您可以考虑 SIMD 矢量化，请参阅openMP SIMD 指令、矢量类库、xsimd 库等。这可能需要重新考虑循环。