高精度MPFR程序崩溃
高精度MPFR程序崩溃
提问于 2015-08-27 16:59:48
最近，我编写了一个计算pi到指定数字的程序。数字数必须作为第一个命令行参数传递。
当运行的数字值低于300，它的工作很好。但是，当以更大的数字值运行时，它会崩溃，但会出现以下异常。
../../src/init2.c:52: MPFR assertion failed: p >= 2 && p <= ((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1))
Aborted (core dumped)我相信这是由于MPFR库中的最大精度所致，但是，我不明白如何更改该最大值或绕过它。
这是我的代码pi.c
#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>
#include <math.h>
void CalculatePi(mpfr_t* pi, int precision, unsigned long long digits, unsigned long long* iterations)
{
//this function uses the gauss-legendre algorithm
unsigned long long loopCount = 0;
mpfr_set_default_prec(precision);
mpfr_t epsilon;
mpfr_t oldA;
mpfr_t oldB;
mpfr_t oldT;
mpfr_t oldP;
mpfr_t A;
mpfr_t B;
mpfr_t T;
mpfr_t P;
mpfr_init2(epsilon, precision);
mpfr_init2(oldA, precision);
mpfr_init2(oldB, precision);
mpfr_init2(oldT, precision);
mpfr_init2(oldP, precision);
mpfr_init2(A, precision);
mpfr_init2(B, precision);
mpfr_init2(T, precision);
mpfr_init2(P, precision);
//epsilon = pow((1 / 10), digits);
mpfr_set_ui(epsilon, 1, MPFR_RNDD);
mpfr_div_ui(epsilon, epsilon, 10, MPFR_RNDD);
mpfr_pow_ui(epsilon, epsilon, digits, MPFR_RNDD);
//oldA = 1;
mpfr_set_ui(oldA, 1, MPFR_RNDD);
//loldB = (1 / sqrt(2));
mpfr_sqrt_ui(oldB, 2, MPFR_RNDD);
mpfr_ui_div(oldB, 1, oldB, MPFR_RNDD);
//oldT = (1 / 4);
mpfr_set_ui(oldT, 1, MPFR_RNDD);
mpfr_div_ui(oldT, oldT, 4, MPFR_RNDD);
//oldP = 1;
mpfr_set_ui(oldP, 1, MPFR_RNDD);
while (1)
{
// A = ((oldA + oldB) / 2)
mpfr_add(A, oldA, oldB, MPFR_RNDD);
mpfr_div_ui(A, A, 2, MPFR_RNDD);
// B = sqrt(oldA * oldB)
mpfr_mul(B, oldA, oldB, MPFR_RNDD);
mpfr_sqrt(B, B, MPFR_RNDD);
// T = (oldT - (oldP * pow((oldA - A), 2)))
mpfr_sub(T, oldA, A, MPFR_RNDD);
mpfr_pow_ui(T, T, 2, MPFR_RNDD);
mpfr_mul(T, oldP, T, MPFR_RNDD);
mpfr_sub(T, oldT, T, MPFR_RNDD);
// P = (oldP * 2)
mpfr_mul_ui(P, oldP, 2, MPFR_RNDD);
// oldA = A
mpfr_set(oldA, A, MPFR_RNDD);
// oldB = B
mpfr_set(oldB, B, MPFR_RNDD);
// oldT = T
mpfr_set(oldT, T, MPFR_RNDD);
// oldP = P
mpfr_set(oldP, P, MPFR_RNDD);
loopCount++;
mpfr_add(*pi, A, B, MPFR_RNDD);
mpfr_pow_ui(*pi, *pi, 2, MPFR_RNDD);
mpfr_mul_ui(T, T, 4, MPFR_RNDD);
mpfr_div(*pi, *pi, T, MPFR_RNDD);
printf("Pi is ");
mpfr_out_str (stdout, 10, digits, *pi, MPFR_RNDD);
putchar ('\n');
//if (abs(A - B) < epsilon)
//break;
mpfr_sub(P, A, B, MPFR_RNDN);
mpfr_abs(P, P, MPFR_RNDN);
if(mpfr_lessequal_p(P, epsilon) != 0)
break;
}
//pi = (pow((A + B), 2) / (T * 4));
//mpfr_add(*pi, A, B, MPFR_RNDD);
//mpfr_pow_ui(*pi, *pi, 2, MPFR_RNDD);
//mpfr_mul_ui(T, T, 4, MPFR_RNDD);
//mpfr_div(*pi, *pi, T, MPFR_RNDD);
iterations = &loopCount;
mpfr_clear(epsilon);
mpfr_clear(oldA);
mpfr_clear(oldB);
mpfr_clear(oldT);
mpfr_clear(oldP);
mpfr_clear(A);
mpfr_clear(B);
mpfr_clear(T);
mpfr_clear(P);
}
int main(int argc, char* argv[])
{
unsigned long long digits = strtoull(argv[1], NULL, 10);
unsigned long long iterations = 0;
//precision = log(10 ^ digits)
int precision = ceil(log2(pow(10.0, (double)digits)));
mpfr_t pi;
mpfr_init2(pi, precision);
CalculatePi(&pi, precision, digits, &iterations);
printf("-----FINAL RESULT REACHED-----\n");
printf("Pi is ");
mpfr_out_str (stdout, 10, digits, pi, MPFR_RNDD);
putchar ('\n');
mpfr_clear(pi);
return 0;
}我用以下命令使用gcc编写了这个程序：
gcc -o pi pi.c -lgmp -lmpfr -lm回答 1
Stack Overflow用户
回答已采纳
发布于 2015-08-31 07:34:28
这个问题是由于这样一个事实:对于309位或更多的数字(这或多或少对应于您的帖子中的大约300位)，pow (10.0, digits)溢出，给出了无穷大，因此ceil(log2(pow(10.0, (double)digits)))也给出了无穷大，当转换为整数时，这是一种未定义的行为。溢出是double类型的限制，与MPFR无关。为了避免大量的数值(例如溢出)，可以用n*log2 2(10)代替数学表达式log2(10^n)，这在实际中是不会溢出的。
注意: MPFR的指数范围比double大得多，因此在MPFR中而不是double中计算log2(10^ n )也将避免n的合理值的溢出，但n*log2 2(10)形式在任何情况下都更安全，应该更快地计算。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/32255374
复制相关文章
相似问题
腾讯云开发者
