如何使用Boost Multiprecision与Boost Math实现大数值指数运算？

我在运行有偏蒙特卡洛（Monte Carlo）模拟时，需要处理能量输出值，核心计算是exp(-beta*energy)——其中energy为负数，beta取值约为100。当energy = -90时（此时指数项为900），std::exp开始输出inf。以下是我的测试代码和结果：

测试代码：

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <iostream>
#include <cmath>

int main()
{
    // using namespace boost::multiprecision;

    double u = std::exp (900.0);

    boost::multiprecision::cpp_dec_float_50 v = 2.71;

    // loop for e^900
    for(unsigned i = 1; i < 901; ++i){
        v *= 2.71;
    }
    
    boost::multiprecision::cpp_dec_float_50 x = boost::math::expm1 (900.0);

    std::cout << "u = " << u << std::endl; 
    std::cout << "v = " << v << std::endl; 
    std::cout << "x = " << x << std::endl;  

    return 0;
}

测试结果：

u = inf
v = 1.27447e+390
x = inf

解决方法 #

问题根源：

1. std::exp基于双精度浮点数，其最大可处理的指数约为709（std::exp(709.78)接近双精度最大值1.79769e+308），900远超这个范围，因此返回inf。
2. boost::math::expm1是为计算exp(x)-1设计的，仅针对x接近0的场景做精度优化，不适合处理大指数输入，所以也返回inf。

正确实现方式是直接使用Boost多精度类型对应的exp函数，替代低效且精度差的手动循环乘近似值：

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <iostream>

int main()
{
    using namespace boost::multiprecision;

    // 计算e^900
    cpp_dec_float_50 exp_result = exp(cpp_dec_float_50(900.0));
    std::cout << "e^900 = " << exp_result << std::endl;  

    // 对应业务场景：exp(-beta*energy)，示例参数beta=100，energy=-90
    cpp_dec_float_50 beta = 100.0;
    cpp_dec_float_50 energy = -90.0;
    cpp_dec_float_50 business_result = exp(-beta * energy);
    std::cout << "exp(-beta*energy) = " << business_result << std::endl;

    return 0;
}

说明 #

• cpp_dec_float_50支持更大范围的数值表示，能轻松处理900级别的大指数计算，且精度远高于手动循环乘近似值的方式。
• 必须确保调用Boost多精度库的exp函数，可通过using namespace boost::multiprecision;或明确写boost::multiprecision::exp实现，避免调用标准库的std::exp。

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