GNU MP Bignum 库的数值问题

Posted 2023-02-23

【中文标题】GNU MP Bignum 库的数值问题

【英文标题】：Numerical Problems with GNU MP Bignum Library

【发布时间】：2012-07-22 22:08:07

【问题描述】：

我正在做一些数字运算，这需要高精度算术。我正在使用 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;

【问题讨论】：

9.3 和 9.8 不能用二进制精确表示。

这似乎只是 15 的幂，而不是 100 的幂。

你是对的 DeadMG。我使用的是 100 的幂，但后来意识到 15 足以说明我的观点。已编辑。

@Mystical：感谢您的评论。老实说，我不明白你为什么这么说，但如果它有所不同，问题仍然存在于 93 和 98 上。这些只是整数 - 当然它们可以精确表示。

由于9.3 和9.8 在double 中无法表示，因此当它们被精确扩展到您想要的任何GMP 精度时，它们将保持+/- 10^-16 错误。这是主要问题。如果它不适用于偶数整数，那么还有一个单独的问题。我怀疑您没有正确使用该库，但我没有使用 GMP，所以我无法判断。

【参考方案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

【讨论】：

没注意到，但我明白你的意思——9.3 和 9.8 是双打，而不是 mpf_class。

@Oleg2718281828 你的意思是定义：mpf_class * x = new mpf_class; ?

@Oleg2718281828 感谢您的贡献。它确实有所作为 - 误差减少到 10^9 的数量级。不过我真的需要它为零。

那是因为它仍然使用9.3 和9.8。请记住，它们被视为double，因此它们在转换为mpf_class之前进行了相应的四舍五入。

@mga 我认为你需要有理数而不是浮点数

【参考方案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 位，所以这仅意味着两种方式都以相同的方式累积误差）。

【讨论】：

我相信a = 9.3; 仍然需要更改。也许它会接受a = "9.3";。虽然我在文档中找不到赋值运算符重载列表。

即使我在声明中有a(9.3,200), b(9.8,200)，输出保持不变。在构造函数中使用字符串实际上更糟糕。

hmm... 缺乏关于 GMP 精确行为的文档肯定无济于事。可以肯定的是，除非它们完全可表示，否则您不能使用浮点文字。

啊，呃，我忘了设置c的精度:(

想一想，这个特定的示例通过了，因为它使用相同的（由于double 舍入不正确）值9.3 和9.8 用于两个计算。如果将结果与不同的库进行比较，它们仍然只有大约 16 位数的准确度。如果 GMP 不接受字符串作为参数，那么您几乎必须手动以全精度除以 10。