计算大数值的 exp(x)/2 避免浮点数溢出

问题描述

在编写双曲正弦函数 sinh(x) 的实现时，遇到一个边界情况：当 x 处于 709.782 到 710.475 范围内时，直接计算 exp(x)/2 会导致浮点数溢出。

在 IEEE 754 双精度浮点数（double）中：

• DBL_MAX 约为 1.798 × 10³⁰⁸
• exp(x) 能够保持有限值的最大 x 约为 709.782
• sinh(x) 和 cosh(x) 能够保持有限值的最大 x 约为 710.475

当 x >= 20 左右时，sinh(x) 和 cosh(x) 都可以用 exp(x)/2 来近似。但当 x 接近上限时，直接计算 exp(x) 会溢出，后续的 /2 操作无法挽回。

解决方案

方法一：数学恒等式变换（基础版）

使用数学恒等式：

eˣ/2 = e^(x - v) × eᵛ/2

选择 v = ln(2)（约 0.693147），则：

double sinh_approx_v1(double x) {
    static const double v = 0.69314718055994530941; // ln(2)
    static const double scale = 2.0 / 2; // eᵛ/2
    return exp(x - v) * scale;
}

但这种方法存在精度问题，误差高达 495 ULP，因为减法 x - v 的误差会被 exp() 函数放大。

方法二：精确减法优化

选择特殊的 v 值，使得减法 x - v 是精确的：

double sinh_approx_v2(double x) {
    static const double v = 0x1.62e42fefa3c00p-1; // 0.69314718056000401
    static const double scale = 2.000000000000117415215 / 2; // eᵛ/2
    return exp(x - v) * scale;
}

这种方法将误差降低到约 1.792 ULP，因为减法操作现在是精确的。

方法三：优化缩放常数

进一步优化，选择使 eᵛ 的二进制表示中有更多尾随零的 v 值：

double sinh_approx_v3(double x) {
    static const double v = 0x1.62e4307c58800p-1; // 0.69314719694034466
    static const double scale = 0x1.000000465a709000p+1 / 2; // eᵛ/2
    return exp(x - v) * scale;
}

这种方法将误差进一步降低到约 1.008 ULP。

方法四：使用 FMA 优化

将缩放因子分解为 1 + scale_minus_1 的形式，利用 FMA 指令：

double sinh_approx_v4(double x) {
    static const double v = 0x1.62e42fefa3c00p-1;
    static const double scale_minus_1 = 5.870760753453405e-14;
    double y = exp(x - v);
    return fma(y, scale_minus_1, y); // 或者 y * scale_minus_1 + y
}

这种方法误差约为 1.006 ULP，且计算效率较高。

方法五：幂运算分解

使用数学变换避免大数运算：

double sinh_approx_v5(double x) {
    double s = exp(x / 2);
    return 0.5 * s * s;
}

这种方法简单直接，误差约为 1.918 ULP，适合对精度要求不极高的场景。

方法六：双精度分解法（Mark Dickinson 方法）

对于 x 在 [1024×ln(2), 1025×ln(2)] 范围内的情况：

double sinh_approx_v6(double x) {
    static const double c1 = 0x1.62e42fefa39efp+9;
    static const double c2 = 0x1.acp-46;
    return exp(x - c1 - c2) * 0x1p+1023;
}

这种方法误差极低（约 0.612 ULP），且所有运算都是精确的。

性能与精度比较

根据测试结果（使用 100,000,007 个样本点）：

方法	最大误差 (ULP)	平均时间 (μs)
sinh_approx_v6	0.600	0.0162
sinh_approx_v4	1.006	0.0253
sinh_approx_v5	1.918	0.0169
标准库 sinh	1.918	0.0234

实现建议

高精度需求场景

对于需要最高精度的应用，推荐使用 Mark Dickinson 的方法：

double precise_sinh(double x) {
    if (x < 710.12928648366403) {
        static const double c1 = 0x1.62e42fefa39efp+9;
        static const double c2 = 0x1.abc9e3b39803fp-46;
        return exp(x - c1 - c2) * 0x1p+1023;
    } else {
        static const double c1 = 0x1.633ce8fb9f87ep+9;
        static const double c2 = -0x1.3ae594e9bd8bp-45;
        return exp(x - c1 - c2) * 0x1p+1023 * 2;
    }
}

平衡精度与复杂度

对于大多数应用，简单的幂运算分解方法已足够：

double simple_sinh(double x) {
    double s = exp(x / 2);
    return 0.5 * s * s;
}

注意事项

1. 避免使用 0.5 / exp(-x)，在处理亚正规数时可能导致精度损失
2. exp(x - log(2)) 由于减法误差会被放大，精度较差
3. 所有方法都应进行边界检查，确保输入在有效范围内

总结

计算大数值的 exp(x)/2 时需要特别注意浮点数溢出问题。通过数学变换和精心选择的常数，可以实现在整个 double 范围内的精确计算。选择哪种方法取决于具体应用对精度和性能的要求。

对于大多数情况，简单的 0.5 * exp(x/2) * exp(x/2) 已经足够；对于需要最高精度的科学计算，推荐使用基于双精度分解的 Mark Dickinson 方法。