c - 使用 GMP 从公式计算非常大的整数

Question

我很难弄清楚如何处理我遇到的问题。 作为复杂公式的一部分，我需要计算一个快速溢出的部分，即结果达到 ~ 1.59*10^(1331)(用 mathematica 计算)。 当然这不在double的范围内。然后我在考虑使用 long double，在我的 gcc 4.6.3 的 linux 系统上它是 16byte。

1) double (8byte) 的可能范围最大为 10^(308)。我说的对吗，那个 long double 然后会增加实际精度，但不会增加可能的数值范围？记得听说可以是or，取决于系统和编译器。真的吗 ？至少当我尝试用 long double 计算我的值时，我仍然得到 NaN。

2.) 然后我在寻找一种方法来实际计算这些结果，我找到了 GNU gmp。我听说您可以表示非常大的整数，我认为这可能会有所帮助。但是，阅读文档，似乎与

 mpz_t x;
 mpz_init(x);
 mpz_set_*(x,#);

我可以为 gmp 整数数据类型赋值，但为了做到这一点，我“只能”选择赋值可以由内置数据类型表示的值，如 double 或 ( u/s)int 等。我发现如何分配非常大的数字是使用 mpz_set_str() 从字符串中分配数字。 我将如何分配一个复杂计算结果的数字？ 简单来说，公式如下所示:

long double res1,res2=0.0;
int a,b;
a=780;
b=741;
float d,d1,o,s; // can be values in [0.01,100]

res1=(2*(pow(b,2)*pow(E,b*(o + s))*(pow(d1,2) + pow(E,a*s)*(-1 + pow(E,a*o)) + pow(d,2)*(-1 + pow(E,a*s))) + pow(a,2)*pow(E,a*(o + s))*(pow(E,b*s)*(pow(E,b*o) + (-1 + d)*(1 + d + b*o)) + (-d + d1)*(d + d1 + b*o + b*d*s)) - a*b*(pow(d1,2)*(pow(E,a*(o + s)) + pow(E,b*(o + s))) + pow(E,a*o + (a + b)*s)*(-2 + 2*pow(E,b*o) - b*o) + d1*pow(E,a*s)*(-pow(E,b*o) + pow(E,a*o)*(1 + b*o)) + pow(d,2)*pow(E,a*s)*(-pow(E,b*o) + pow(E,b*(o + s)) + pow(E,a*o)*(-1 + pow(E,b*s) - b*s)) +  d*(-(d1*pow(E,b*(o + s))) + (1 + d1)*pow(E,b*o + a*s) - pow(E,a*s + b*(o + s)) + pow(E,b*s + a*(o + s))*(1 + b*o) + pow(E,a*(o + s))*(-1 - b*o + b*d1*s)))))......;

res2也会是这种，最后我要计算res1/res2，通常会变成一个很小的数。

我曾考虑拆分公式并向 mpg_z 添加项，以免每个项超出双倍范围，但由于公式又长又复杂，这几乎是不可能的。

所以总而言之，问题是，我的中间结果可能会变得如此之大，以至于没有任何数据类型能够存储它们，所以我无法将它分配给 mpz 并摆脱这个问题。

我知道我想计算一个 double 值并实际使用 mpz_t 作为整数。据我所知，这是存储如此大数据的唯一方法，因为 mpf_t 只能处理 float 类型。老实说，我对 gmp 中的表示仍然感到困惑。

有什么想法可以解决这个问题吗？

Answer 1

问题 1 Long double 允许处理比 double 更大的数(指数和尾数精度)。但是你必须认为，如果你的目标是存储大整数，你的 1e308 数量级没有任何意义；您只需关心有效数字精度的大小，52/53 位( double )或 64 位(x86 扩展精度)。如果您尝试将它与更大的整数一起使用，您将获得正确的数量级，但精确值将会丢失(当使用整数计算时，人们通常比使用近似数字时更关心这一点)。

问题 2 使用 GMP 是一个不错的选择。其他图书馆也存在；对于较小的值，我经常使用 libqd它具有扩展的固定精度并且非常快 但这不足以解决您自己的问题。现在你的问题是关于设置值:

• 使用数字的字符串版本通常不是一个好主意(您应该仅将其保留用于输入/输出目的)；这是一个涉及基础转换和逐位处理的缓慢操作
• 尽可能多地使用 GMP 类型是可行的方法(除非您真的关心速度并且完全控制要使用 native 类型计算的值的预期范围
• 如果您的公式太长，我也无能为力。但是使用GMP并没有那么难，你真的不能转换你的配方吗？您的公式中是否有一些可以嵌入循环中的逻辑？或许您可以编写一个快速而复杂的 Python 脚本，使用 GMP 将您的公式转换为一段 C 代码？

现在我不完全明白为什么要使用mpz_t 而不是mpf_t。这种类型实现了任意长 float ；您是否注意到可以使用 mpf_set_default_prec 设置精度？