如果您正在寻找舍入误差的良好界限，则不一定需要任意精度库。您可以改用运行错误分析。

我无法找到一个好的在线参考资料，但这一切都在 Nick Higham 的书“数值算法的准确性和稳定性”的第 3.3 节中进行了描述。这个想法很简单：

1. 重构您的代码，以便您在每行上都有一个算术运算的单一分配。
2. 对于每个变量，例如x，创建一个变量，当分配一个常量时，该变量x_err初始化为零。x
3. 对于每个操作，例如，使用浮点算术的标准模型z = x * y更新变量以及结果和运行错误和。z_errzx_erry_err
4. 然后，您的函数的返回值也应该_err附加一个相应的值。这是您的总舍入误差的数据相关界限。

棘手的部分是第 3 步。对于最简单的算术运算，您可以使用以下规则：

• z = x + y->z_err = u*abs(z) + x_err + y_err
• z = x - y->z_err = u*abs(z) + x_err + y_err
• z = x * y->z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
• z = x / y->z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
• z = sqrt(x)->z_err = u*abs(z) + x_err/(2*abs(z))

哪里u = eps/2是单位舍入。+是的，和的规则-是相同的。op(x)使用应用于 的结果的泰勒级数展开，可以轻松提取任何其他操作的规则op(x + x_err)。或者你可以尝试谷歌搜索。或者使用尼克海厄姆的书。

例如，考虑以下 Matlab/Octave 代码，该代码使用 Horner 方案a在某个点处评估系数中的多项式：x

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        s = a(k) + x*s;
    end

第一步，我们将两个操作拆分为s = a(k) + x*s：

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        z = x*s;
        s = a(k) + z;
    end

然后我们引入_err变量。请注意，输入a和x假定是准确的，但我们也可以要求用户为a_err和传递相应的值x_err：

function [ s , s_err ] = horner ( a , x )
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = ...;
        s = a(k) + z;
        s_err = ...;
    end

最后，我们应用上述规则来获取错误项：

function [ s , s_err ] = horner ( a , x )
    u = eps/2;
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = u*abs(z) + s_err*abs(x);
        s = a(k) + z;
        s_err = u*abs(s) + z_err;
    end

请注意，由于我们没有a_error x_err，例如它们被假定为零，因此在错误表达式中将简单地忽略相应的项。

瞧！我们现在有一个霍纳方案，它在结果旁边返回一个与数据相关的误差估计（注意：这是误差的上限）。

附带说明一下，由于您使用的是 C++，因此您可能会考虑为浮点值创建自己的类，该类包含该_err术语并重载所有算术运算以更新这些值，如上所述。对于大型代码，这可能是更容易的路线，尽管计算效率较低。话虽如此，您也许可以在网上找到这样的课程。一个快速的谷歌搜索给了我这个链接。

PS 请注意，这一切都只适用于严格遵守 IEEE-754 的机器，即所有算术运算都精确到 $\pm u$。这种分析还给出了比使用区间算术更严格、更现实的界限，因为根据定义，您不能用浮点数表示数字 $x(1 \pm u)$，即您的区间只会舍入到数字本身。$\pm u$. This analysis also gives a tighter, more realistic bound than using interval arithmetic since, by definition, you can not represent a number $x(1 \pm u)$ in floating point, i.e. your interval would just round to the number itself.

Victor Shoup 的 NTL是一个很好的可移植和开源库，用于任意精度浮点运算（以及其他许多），它以 C++ 源代码形式提供。

较低级别是GNU 多精度 (GMP) Bignum 库，它也是一个开源包。

NTL 可以与 GMP 一起使用，需要更快的性能，但 NTL 提供了自己的基本例程，如果您“不关心速度”，这些例程肯定是可用的。GMP 声称是“最快的 bignum 库”。GMP 主要是用 C 编写的，但有一个 C++ 接口。

补充：虽然区间算术可以自动给出准确答案的上限和下限，但这并不能准确测量“标准”精度计算中的误差，因为区间大小通常会随着每次操作而增长（相对或绝对误差意义）。

对于舍入误差或离散化误差等，查找误差大小的典型方法是计算额外的精度值并将其与“标准”精度值进行比较。仅需要适度的额外精度即可将误差大小本身确定为合理的精度，因为“标准”精度中的舍入误差本身比额外精度计算中的舍入误差要大得多。

可以通过比较单精度和双精度计算来说明这一点。请注意，在 C++ 中，中间表达式总是（至少）以双精度计算，所以如果我们想说明“纯”单精度计算是什么样的，我们需要以单精度存储中间值。

C 代码片段

    float fa,fb;
    double da,db,err;
    fa = 4.0;
    fb = 3.0;
    fa = fa/fb;
    fa -= 1.0;

    da = 4.0;
    db = 3.0;
    da = da/db;
    da -= 1.0;

    err = fa - da;
    printf("Single precision error wrt double precision value\n");
    printf("Error in getting 1/3rd is %e\n",err);
    return 0;

上面的输出（Cygwin/MinGW32 GCC 工具链）：

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08

因此，误差是关于将 1/3 舍入为单精度所期望的。一个人不会（我怀疑）关心在错误中得到超过几个小数位是正确的，因为错误的测量是为了幅度而不是精确度。

GMP（即 GNU 多精度库）是我所知道的最好的任意精度库。

我不知道有任何编程方法可以测量任意浮点函数结果中的误差。您可以尝试的一件事是使用区间算术计算函数的区间扩展。在 C++ 中，您必须使用某种库来计算间隔扩展；一个这样的库是Boost Interval Arithmetic Library. 基本上，为了测量误差，您将作为参数提供给您的函数间隔，其宽度为单位舍入的 2 倍（大致），以感兴趣的值为中心，然后您的输出将是间隔的集合，宽度为这会给你一些保守的误差估计。这种方法的一个困难在于，以这种方式使用的区间算术可能会高估大量误差，但这种方法是我能想到的最“程序化”的方法。

通过区间分析可以实现严格和自动的误差估计。您使用间隔而不是数字。例如加法：

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]

舍入也可以严格处理，请参阅舍入区间算术。

只要您的输入包含较窄的间隔，估计值就可以并且计算起来非常便宜。不幸的是，错误经常被高估，请参阅依赖问题。

我不知道任何任意精度区间算术库。

区间算术能否满足您的需求取决于您手头的问题。