是否可以使此表达式更准确并且不返回负值？

用

4*z - 3 - z^2 = (3-z)*(z-1)

和使用

e^(-aT)-1 = expm1(-aT)

它以完全相对精度计算此因子为 。-aT*(1+O(aT))

使用 gnuplot 可视化差异：

原始函数出现问题的根本原因是减法抵消，即对两个大小相似的量进行有效减法。关于如何处理这种情况，有一些一般的启发式方法：（1）对于涉及接近零的计算，改用。（2） 将减法转换为除法或乘法 （3） 在发生减法抵消的地方使用融合乘加法，前提是所涉及的两个操作数之一是一个乘积。这之所以有效，是因为 FMA 将双倍宽度未舍入乘积合并到总和中。exp()expm1()

在提出合适的缓解重新排列后，我注意到精度损失持续存在，小于 ^2-6 小时。对于非常接近于零的参数，此类问题的标准方法是尝试泰勒级数展开。我使用配置为 130 位十进制数字的多精度库，以数值方式计算了泰勒展开的前六项。更快的替代方案是使用 Wolfram Alpha。a * T

我应该指出，我以行人的方式处理了减少错误的替换。Lutz Lehmann 找到了一种更优雅、更高效且不需要 FMA 的重新排列;看他的回答。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define USE_NJ (1)
#define FUNC   accurate_func

double accurate_func (double a, double T)
{
    double r, u = a * T;
    // use Taylor series expansion near zero
    if (u <= 0.015625) {  // 1/64
        double u3 = u * u * u;
        double p =  -1.0 / 160.0;
        p = p * u + 31.0 / 1260.0; 
        p = p * u -  1.0 / 12.0;
        p = p * u +  7.0 / 30.0;
        p = p * u -  1.0 / 2.0;
        p = p * u +  2.0 / 3.0;
        r = p * u3;
    } else {
#if USE_NJ // use exmp1() and fma() to mitigate subtractive cancellation
        double v = expm1 (u);
        double w = exp (-2 * u);
        r = fma (-3 * v * v, w, 2 * fma (-v, w, u));
#else // Lutz Lehmann's elegant arrangement
        double s = expm1 (-u);
        r = (2.0 - s) * s + 2 * u;
#endif
    }
    return sqrt (r);
}

double original_func (double a, double T)
{
    double u = a * T;
    return sqrt (4* exp (-u) - 3 - exp (-2 * u) + 2*u);
}

int main (void)
{
    double mpref [13] = {
        5.7981245368580744e-1, // 2^0
        2.4133627509989329e-1, // 2^-1
        9.3126111123497739e-2, // 2^-2
        3.4450359460776940e-2, // 2^-3
        1.2463894628395316e-2, // 2^-4
        4.4581489320727356e-3, // 2^-5
        1.5854164404162396e-3, // 2^-6
        5.6217040623145288e-4, // 2^-7
        1.9904830179668834e-4, // 2^-8
        7.0425736778377448e-5, // 2^-9
        2.4908375624050728e-5, // 2^-10
        8.8080530945851612e-6, // 2^-11
        3.1144021359546032e-6, // 2^-12
    };

    double mpref2 [13] = {
        5.7981245368580744e-1, // 2^0
        1.2341091904923634e+0, // 2^1
        2.2523159398554711e+0, // 2^2
        3.6057373362429543e+0, // 2^3
        5.3851648489290173e+0, // 2^4
        7.8102496759066576e+0, // 2^5
        1.1180339887498948e+1, // 2^6
        1.5905973720586866e+1, // 2^7
    };

    int i = 0;
    double a = 1.0;
    do {
        double T = 1.0 / (1 << i);
        printf ("i=%2d a*T=%23.16e  func=%23.16e  relerrfunc=% 15.8e\n", 
                i, a*T, FUNC(a,T), (FUNC(a,T) - mpref[i]) / mpref[i]);
        i++;
    } while (i < 13);
    printf ("\n");
    i = 0;
    do {
        double T = 1.0 * (1 << i);
        printf ("i=%2d a*T=%23.16e  func=%23.16e  relerrfunc=% 15.8e\n", 
                i, a*T, FUNC(a,T), (FUNC(a,T) - mpref2[i]) / mpref2[i]);
        i++;
    } while (i < 8);

    return EXIT_SUCCESS;
}

我只是一个随机的 Julia 人，不是数字学家，但我想贡献 Herbie，一个自动提高浮点精度的工具，给你的东西：

function f(a, T)
    t1 = Float64((T ^ 4.0) * -0.5)
    t2 = a ^ 4.0
    t3 = Float64(Float64(T * -2.0) + Float64(T * 2.0))
    t4 = a ^ 3.0
    t5 = Float64((T ^ 3.0) * 0.6666666666666666)
    t6 = log(exp(Float64(0.23333333333333334 * (Float64(T * a) ^ 5.0))))

    u1 = fma(t4, t5, t6)
    u2 = fma(t3, a, u1)
    return fma(t1, t2, u2)
end

（赫比可以直接生成茱莉亚;我只是将输出从单行重新格式化为更可分析的内容。

它甚至可以为你生成一个证明，证明为什么转换是正确的。