GNU MP Bignum 库的数值问题

Question

我正在做一些数字运算，这需要高精度算术。我正在使用 GNU MP 库，according to the GMP manual:

“浮点数或简称Float，是具有有限精度指数的任意精度尾数。”

虽然尾数应该具有任意精度，但我仍然遇到精度问题。与其用我的实际代码让你厌烦，不如用一个近乎最小的工作示例来说明我的问题。代码计算 9.3^15、9.8^15 和 (9.3*9.8)^15。在我的机器上，(9.3^15)*(9.8^15) 和 (9.3*9.8)^15 的值从第 16 位开始开始不同，在这种情况下导致（大约）4.94* 的错误10^13。

任何帮助将不胜感激。代码如下。

#include <gmp.h>
#include <gmpxx.h>

#include <iostream>
#include <iomanip>

int main()
{
    mpf_class x, y, z;
    x = y = z = 1.0;

    for (int i = 0; i < 15; i++)
    {
        x *= 9.3;
        y *= 9.8;
        z *= 9.3*9.8;
    }

    std::cout << z - x*y << std::endl;

    return 0;
}

Answer 1

您看到的问题是由于 9.3 * 9.8 的计算结果近似。请将文字更改为 mpf_class 的实例：

mpf_class a, b;
a = 9.3;
b = 9.8;

// ...

x *= a;
y *= b;
z *= a * b;

如果您需要无限精度，请考虑使用有理数：

#include <gmp.h>
#include <gmpxx.h>

#include <iostream>
#include <iomanip>

int main()
{
    mpq_class x(1), y(1), z(1), a(93, 10), b(98, 10);

    for (int i = 0; i < 15; i++)
    {
        x *= a;
        y *= b;
        z *= (a * b);
    }

    std::cout << z - x*y << std::endl << z << std::endl;

    return 0;
}

打印

0
7589015305950762920038660273144124106674963183136666693/30517578125000000000000000

Answer 2

问题是你没有明确设置精度，所以你得到默认精度，通常（我认为）64位，因此由于不同的舍入方式不同，最后一位的结果不同计算。这使得大约 20 位数字成为公共前缀（随着计算的增加，差异可能会变得更大）。如果设置更高的精度，

#include <gmp.h>
#include <gmpxx.h>

#include <iostream>
#include <iomanip>

int main()
{
    mpf_class x(1.0,200), y(1.0,200), z(1.0,200), a("9.3",200), b("9.8",200), c(0,200); 
    c = a*b;

    for (int i = 0; i < 15; i++)
    {
        x *= a;
        y *= b;
        z *= c;
    }

    std::cout << z << "\n" << (x*y) << std::endl;
    std::cout << z - x*y << std::endl;

    return 0;
}

这里200位，你会得到更准确的结果：

$ ./a.out 
2.48677e+29
2.48677e+29
-4.80637e-49

因此是大约 80 个十进制数字或近 256 位（大于 199 的 64 的最小倍数）的公共前缀。

精度为 2000，与字符串构造函数的差异为 -2.78942e-588，如果从 double 初始化则为 0（但当然，初始精度限制为 53 位，所以这仅意味着两种方式都以相同的方式累积误差）。