如何提高乘法的精度？答案

【问题标题】：How to improve precision in multiplication?如何提高乘法的精度？
【发布时间】：2014-03-11 15:48:06
【问题描述】：

我正在使用 C++ 实现超限插值算法 (http://en.wikipedia.org/wiki/Transfinite_interpolation)。一切看起来都很好，直到我尝试测试一些小数字时，结果看起来很奇怪且不正确。我想这一定与精度损失有关。代码是

for (int j = 0; j < Ny+1; ++j)
{
    for (int i = 0; i < Nx+1; ++i)
    {
        int imax = Nx;
        int jmax = Ny;
        double CA1 = (double)i/imax;  // s
        double CB1 = (double)j/jmax;  // t
        double CA2 = 1.0-CA1;   // 1-s
        double CB2 = 1.0-CB1;   // 1-t

        point U = BD41[j]*CA2 + BD23[j]*CA1;
        point V = BD12[i]*CB2 + BD34[i]*CB1;
        point UV = 
            (BD12[0]*CB2 + BD41[jmax]*CB1)*CA2
          + (BD12[imax]*CB2 + BD23[jmax]*CB1)*CA1;

        tfiGrid[k][j][i] = U + V - UV;
    }
}

我猜当BD12[i]（或BD34[i]或BD41[j]或BD23[j]）非常小时，舍入误差或某些东西会累积并变得可以忽略不计。任何想法如何处理这种情况？

PS：尽管类似的问题已经被问了数百万次。我仍然无法弄清楚这与我的乘法或除法或减法有关吗？

【问题讨论】：

标签： c++ precision

【解决方案1】：

除了安托万提出的几点（都非常好）：可能值得记住的是，将两个值相加不同的数量级会造成很大的损失精确。例如，如果CA1 小于大约1E-16， 1.0 - CA1 可能仍然是 1.0，即使它只是稍大一点，你会失去相当多的精度。如果这是问题所在，您应该能够通过以下方式将其隔离在内部循环中放置一些打印语句，然后查看您正在添加的值（或者甚至可能使用调试器）。

如何处理是另一个问题。可能有一些您正在尝试做的事情的数值稳定算法；我不知道。否则，您可能必须检测动态问题，并重写方程以避免它，如果它发生。（例如，检测CA1 是否太小另外，您可以检查1.0 / CA1 是否超过几千，或几百万，或者无论你有多精确可以松动。）

【讨论】：

【解决方案2】：

C/C++ 中内置的算法的准确性是有限的。当然，错误会在您的情况下累积。

您是否考虑过使用提供更高精度的库？也许看看https://gmplib.org/

一个说明更高准确性的简短示例：

double d, e, f;
d = 1000000000000000000000000.0;
e = d + 1.0;
f = e - d;
printf( "the difference of %f and %f is %f\n", e, d, f);

这不会打印 1 而是 0。使用 gmplib 的代码如下所示：

#include "gmp.h"

mpz_t m, n, r;

mpz_init( m);
mpz_init( n);
mpz_init( r);

mpz_set_str( m, "1 000 000 000 000 000 000 000 000", 10);
mpz_add_ui( n, m, 1);
mpz_sub( r, n, m);
printf( "the difference of %s and %s is %s (using gmp)\n",
    mpz_get_str( NULL, 10, n),
    mpz_get_str( NULL, 10, m),
    mpz_get_str( NULL, 10, r));

mpz_clear( m);
mpz_clear( n);
mpz_clear( r);

这将返回 1。

【讨论】：

嗯，我应该大量更改我的代码吗？或者我可以简单地使用gmp_double 或类似的东西吗？
谢谢你的例子。 gmp 和 mpfr 听了很多，但一直没有找到机会使用它。 :)
我认为交换数据类型和对 gmplib 等价物的操作就足够了。所以你可以使用 mpz_t 代替 double 和 mpz_mul(r, n, m) 代替 r = n * m。不要忘记使用 mpz_init() 初始化变量。
如果你看他的代码，没有什么会导致错误累积。另一方面，在没有任何累加的情况下，加减不同幅度的值可能会导致精度损失很大。

【解决方案3】：

您的算法似乎没有通过在每次迭代中重复使用以前的计算来累积错误，因此如果不查看您的数据就很难回答您的问题。

基本上，您有 3 个选项：

提高数字的精度：x86 CPU 可以处理float（32 位）、double（64 位），通常还有long double（80 位）。除此之外，您必须使用“软”浮点，其中所有操作都是在软件而不是硬件中实现的。有一个很好的 C 库可以做到这一点：MPFR 基于GMP，GNU 建议使用 MPFR。我强烈建议使用更易于使用的 C++ 包装器，例如 boost multiprocesion。预计您的计算速度会慢几个数量级。
通过在计算中使用比单个标量数更多的信息来分析精度损失的来源。查看基于 MPFR 的 interval arithmetic 和 MPFI 库。 CADENA 是另一种非常有前途的解决方案，它基于随机改变硬件的舍入模式，并且运行时成本低。
执行静态分析，这甚至不需要运行您的代码并通过分析您的算法来工作。不幸的是，我没有使用此类工具的经验，因此我无法推荐任何东西。

我认为最好的方法是在开发算法时运行静态或动态分析，确定精度问题，通过更改算法或对最不稳定的变量使用更高的精度来解决这些问题 - 而不是其他的，以避免过多的性能运行时的影响。

【讨论】：

【解决方案4】：

数值（即浮点）计算很难精确地进行。您必须特别注意减法，这主要是在您失去精度的地方。在这种情况下，1.0 - CA1 等是可疑的（如果CA1 非常小，您将得到 1）。重新组织你的表达式，维基百科的文章（一个存根！）可能是为了理解（显示对称性，...）和美学而不是数值稳健性而编写的。

搜索关于数值计算的课程/讲义，你应该包括一个介绍性的章节。并查看戈德堡的经典作品What every computer scientist should know about floating point arithmetic。

【讨论】：