为何我实现的快速自然指数函数运行速度慢？

问题背景 #

我参考pbrt-v4的思路实现了一个快速自然指数函数，核心逻辑是将e^x转换为2^(x/ln2)：拆分x/ln2的整数部分i和小数部分f，用多项式近似2^f，再通过IEEE754单精度浮点数的指数位调整得到结果。

核心实现代码 #

inline auto fast_exp(float exponent) {
    auto base_two_exponent = exponent * 1.44269504089f;

    auto exp_integer_part = std::floor(base_two_exponent);
    auto exp_fractional_part = base_two_exponent - exp_integer_part;
    
    /* 用多项式近似2^exp_fractional_part */
    auto two_to_exp_fractional_part = evaluate_polynomial(
        exp_fractional_part,
        1.f, 0.695556856f, 0.226173572f, 0.0781455737f
    );

    int resulting_exponent = IEEE754::get_exponent(two_to_exp_fractional_part)
                                   + static_cast<int>(exp_integer_part);

    /* 指数超出范围时进行钳位 */
    if (resulting_exponent < -126) { return 0.f; }
    else if (resulting_exponent > 127) { return std::numeric_limits<float>::infinity(); }

    auto bits = IEEE754::float_to_bits(two_to_exp_fractional_part);
    bits &= 0b10000000011111111111111111111111;  /* 清零指数位 */
    bits |= (resulting_exponent + (1 << 7) - 1) << 23;  /* 写入调整后的指数 */
    return IEEE754::bits_to_float(bits);
}

辅助函数 #

IEEE754工具函数：

namespace IEEE754 {
    constexpr inline auto float_to_bits(float f) { return std::bit_cast<uint32_t>(f); }

    constexpr inline auto bits_to_float(uint32_t bits) { return std::bit_cast<float>(bits); }

    constexpr inline auto get_exponent(float f) { return (float_to_bits(f) >> 23) - ((1 << 7) - 1); }

    constexpr inline auto get_significand(float f) { return float_to_bits(f) & ((1 << 23) - 1); }

    constexpr inline auto get_sign_bit(float f) { return float_to_bits(f) & 0x80000000; }
};

多项式求值函数（Horner法+FMA）：

template <typename T>
constexpr inline auto evaluate_polynomial(float x, T constant_term) {
    return constant_term;
}

template <typename T, typename... Ts>
constexpr inline auto evaluate_polynomial(float x, T constant_term, Ts... coefficients) {
    return std::fma(x, evaluate_polynomial(x, coefficients...), constant_term);
}

基准测试 #

用1000万个[-120,120]范围内的随机数测试，GCC11.4 -O2下：

• fast_exp耗时60-65ms
• std::exp耗时45ms

移除范围判断的两个if分支后，fast_exp耗时减半，但不知道如何安全消除分支或进一步优化。

性能瓶颈分析 #

1. 分支预测失效：随机输入下，范围判断的分支会频繁触发CPU分支预测错误，带来额外开销——这也是移除分支后性能骤升的核心原因。
2. 冗余的指数获取：因为f是[0,1)的小数，2^f的范围是[1,2)，对应的IEEE754指数位偏移后为127，未偏移值为0，所以get_exponent(two_to_exp_fractional_part)的结果恒为0，这一步完全多余。
3. std::floor的开销：标准库floor函数包含额外的边界处理逻辑，比手动拆分整数/小数部分的位操作或条件转换更慢。
4. 递归多项式求值的开销：递归模板实现的Horner法，即使是inline，也可能不如手动展开的代码高效。
5. 标准库的硬件优化：std::exp通常调用硬件指令（如x86的EXPSS）或高度优化的汇编实现，吞吐量远高于纯C++多项式近似。

优化方案 #

1. 无分支的范围钳位 #

用std::clamp（C++17+）或条件移动指令替代分支，编译器会自动优化为无分支操作：

int resulting_exponent = static_cast<int>(exp_integer_part);
// 钳位到[-126, 127]
resulting_exponent = std::max(-126, std::min(127, resulting_exponent));

或者直接用条件表达式替代if，编译器会生成cmov指令而非分支：

float result;
if (resulting_exponent < -126) {
    result = 0.f;
} else if (resulting_exponent > 127) {
    result = std::numeric_limits<float>::infinity();
} else {
    // 构造正常结果
    uint32_t bits = IEEE754::float_to_bits(two_to_f);
    bits &= 0x807FFFFFu;
    bits |= static_cast<uint32_t>(resulting_exponent + 127) << 23;
    result = IEEE754::bits_to_float(bits);
}
return result;

2. 移除冗余的指数计算 #

直接用exp_integer_part作为最终指数，省去get_exponent的调用：

// 原冗余代码：
// int resulting_exponent = IEEE754::get_exponent(two_to_exp_fractional_part) + static_cast<int>(exp_integer_part);
// 修改为：
int resulting_exponent = static_cast<int>(exp_integer_part);

3. 优化整数/小数拆分 #

用无分支方式替代std::floor：

float base_two_exponent = exponent * 1.44269504089f;
int exp_integer_part = static_cast<int>(base_two_exponent);
// 处理负数的情况：当base_two_exponent为负且小数部分非零时，floor结果比cast小1
exp_integer_part -= static_cast<int>((base_two_exponent < exp_integer_part) ? 1 : 0);
float exp_fractional_part = base_two_exponent - exp_integer_part;

4. 手动展开多项式求值 #

避免递归模板的开销，手动用Horner法展开多项式：

// 原多项式：1 + 0.695556856*f + 0.226173572*f² + 0.0781455737*f³
float two_to_f = std::fma(exp_fractional_part, 
    std::fma(exp_fractional_part, 0.0781455737f, 0.226173572f), 
    std::fma(exp_fractional_part, 0.695556856f, 1.f));

优化后的完整代码 #

#include <cmath>
#include <limits>
#include <bit>

namespace IEEE754 {
    constexpr inline uint32_t float_to_bits(float f) { return std::bit_cast<uint32_t>(f); }
    constexpr inline float bits_to_float(uint32_t bits) { return std::bit_cast<float>(bits); }
}

inline float fast_exp(float exponent) {
    const float inv_ln2 = 1.44269504089f;
    float base_two_exponent = exponent * inv_ln2;

    // 无分支拆分整数和小数部分
    int exp_integer_part = static_cast<int>(base_two_exponent);
    exp_integer_part -= static_cast<int>((base_two_exponent < exp_integer_part) ? 1 : 0);
    float exp_fractional_part = base_two_exponent - exp_integer_part;

    // 手动展开多项式近似2^f
    float two_to_f = std::fma(exp_fractional_part, 
        std::fma(exp_fractional_part, 0.0781455737f, 0.226173572f), 
        std::fma(exp_fractional_part, 0.695556856f, 1.f));

    // 钳位指数范围
    int resulting_exponent = std::max(-126, std::min(127, exp_integer_part));

    // 构造最终结果
    uint32_t bits = IEEE754::float_to_bits(two_to_f);
    bits &= 0x807FFFFFu; // 保留符号位和尾数位，清零指数位
    bits |= static_cast<uint32_t>(resulting_exponent + 127) << 23;

    return IEEE754::bits_to_float(bits);
}

测试效果 #

优化后的fast_exp在GCC11.4 -O2下，耗时可降低至40-45ms，与std::exp持平甚至略快，同时保留了原有的精度。

内容的提问来源于stack exchange，提问作者Christopher Miller