C++计算test_func返回inf而非预期有限值问题求助

问题根源是双精度浮点数溢出，出在I0函数的原始表达式写法：

• 调用test_func(3000,10)时，计算得到的k0≈3000.02，传入I0的两个参数为y0=3000.02、y=3000
• 计算exp(y0+y)时，指数值达到6000.02，而双精度浮点数（double）能容纳的exp最大正指数约为709，超过阈值后会直接返回inf（无穷大）
• 后续inf参与除法、对数运算后结果仍然是inf，最终导致输出异常。

对I0的表达式做等价数学变形，把正指数转换为负指数即可避免溢出：
原表达式：

log( (exp(y0+y)-1)/(exp(y0-y)-1) ) - y

对分子分母分别提取exp公因子后拆分对数，可简化得到完全等价的表达式：

y + log( (1 - exp(-y0 - y)) / (1 - exp(y - y0)) )

变形后所有exp的输入都是负数，计算结果都在0~1区间内，完全不会溢出，精度也符合要求。

仅需要替换I0函数即可：

double I0 (double y0, double y)
{
    double numerator = 1.0 - exp(-y0 - y);
    // 这里用expm1计算exp(x)-1，y0和y差值很小时也能保证精度
    double denominator = -expm1(y - y0);
    return y + log(numerator / denominator);
}

修改后运行即可得到你预期的计算结果。

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