I recently wrote a program to calculate pi to a specified number of digits. The number of digits must be passed as the first command line argument.
When run with a digit value under about 300 it works just fine. However, when run with a larger digit value it crashes with the following exception.
../../src/init2.c:52: MPFR assertion failed: p >= 2 && p <= ((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1))
Aborted (core dumped)
I believe that this is due to a maximum precision in the MPFR library, however, I do not understand how to change that maximum or work around it.
Here is my code 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;
}
I compile the program using gcc with the following command:
gcc -o pi pi.c -lgmp -lmpfr -lm
The problem is due to the fact that for 309 digits or more (which more or less corresponds to the about 300 in your post), pow (10.0, digits) overflows, giving an infinity, so that ceil(log2(pow(10.0, (double)digits))) also gives an infinity, and when converted to an integer, this has an undefined behavior. The overflow is here a limitation of the double type, not related to MPFR. To avoid huge numbers in magnitude (thus the overflow), the mathematical expression log2(10^n) can be replaced by n*log2(10), which will not overflow in practice.
Note: MPFR has a much larger exponent range than double, so that computing log2(10^n) in MPFR instead of double would also avoid the overflow for reasonable values of n, but the n*log2(10) form is safer in any case and should be faster to evaluate.