GNU MP Bignum库的数值问题

Numerical Problems with GNU MP Bignum Library

本文关键字：问题 MP Bignum GNU      更新时间：2023-10-16

我正在做一些数字运算，这需要高精度的算术运算。我使用GNU MP库，根据GMP手册：

浮点数，简称Float，是一个具有有限精度指数的任意精度尾数

虽然尾数的精度应该是任意的，但我仍然会遇到精度问题。与其说我的实际代码让你厌烦，不如说这里有一个几乎最小的工作示例来说明我的问题。代码计算9.3^15、9.8^15和（9.3*9.8）^15。在我的机器上，从第16位开始，（9.3^15）*（9.8^15）和（9.3*9.8）^15的值开始不同，在这种情况下会导致（大约）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;
}

您看到的问题是由于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

问题是，您没有显式设置精度，因此您获得默认精度，通常（我认为）为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位（64的最小倍数大于199）的公共前缀。

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

相关文章：