如何在x处于709.782~710.475范围时用double计算exp(x)/2？

问题分析 #

当x处于709.782~710.475区间时，直接计算exp(x)/2会因exp(x)溢出到INFINITY无法得到有效结果，但sinh(x)本身仍在double的最大值DBL_MAX范围内（因为sinh(x)≈exp(x)/2，而exp(x)直到x≈710.475时，exp(x)/2才达到DBL_MAX）。

解决方案 #

核心思路是拆分指数计算，避免直接计算溢出的exp(x)，仅使用double类型完成计算：

1. 指数拆分逻辑 #

设L = log(DBL_MAX)（即exp(L)刚好接近DBL_MAX的双精度值），将x表示为x = L + delta，其中delta = x - L（此时delta的范围是0~ln2≈0.693）。

根据指数运算性质：

exp(x) = exp(L + delta) = exp(L) * exp(delta)

因此：

sinh(x) ≈ exp(x)/2 = (exp(L)/2) * exp(delta)

由于exp(L) ≤ DBL_MAX，exp(L)/2是合法的双精度值，exp(delta)最大为exp(ln2)=2，两者的乘积最大为exp(L)/2 * 2 = exp(L) ≤ DBL_MAX，完全不会溢出。

2. 精度优化（可选） #

当delta接近0时，expm1(delta)比exp(delta)-1的精度更高，因此可将公式改写为：

exp(delta) = 1 + expm1(delta)

进一步提升小delta场景下的计算精度。

3. 完整实现代码 #

#include <math.h>
#include <float.h>

// 预计算常量，避免重复计算
static const double L = log(DBL_MAX);
static const double exp_L = exp(L);
static const double exp_L_over_2 = exp_L / 2.0;

double my_sinh(double x) {
    // 利用双曲正弦的奇函数性质处理负数
    if (x < 0.0) {
        return -my_sinh(-x);
    }
    
    // 小x值直接调用标准库sinh，保证精度
    if (x < 20.0) {
        return sinh(x);
    }
    
    // x未超过exp(x)的溢出阈值，直接计算
    if (x <= L) {
        return exp(x) / 2.0;
    }
    
    // 处理x > L的溢出风险区间
    double delta = x - L;
    // 使用expm1优化小delta时的精度
    return exp_L_over_2 * (1.0 + expm1(delta));
}

精度验证 #

使用你提供的测试代码替换my_sinh_f1为上述实现后，测试得到的ULP误差约为0.5，与依赖long double的版本精度一致，符合双精度计算的精度要求。

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