MPFR Program Crashes With High Precision

Question

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

If you haven't done so yet, the please take some time to read [the help pages](http://stackoverflow.com/help), including [about how to ask good questions](http://stackoverflow.com/help/how-to-ask). But more importantly you should learn how to create a [Minimal, Complete, and Verifiable Example](http://stackoverflow.com/help/mcve). — Some programmer dude, Aug 27 '15 at 17:02
What is `MPFR_PREC_MAX`? *"Warning! MPFR needs to increase the precision internally, in order to provide accurate results (and in particular, correct rounding). Do not attempt to set the precision to any value near MPFR_PREC_MAX, otherwise MPFR will abort due to an assertion"* failure.http://www.mpfr.org/mpfr-current/mpfr.html — Weather Vane, Aug 27 '15 at 17:10
while improving the program. 1) check that the command line parameter actually exists before using it. 2) suggest a) rewrite the code to not use multiprocessor features while testing b) step through the code with a debugger to see what actual values are being generated. — user3629249, Aug 27 '15 at 18:40
I suspect that you have an overflow in `pow (10.0, digits)`, thus with undefined behavior. log2(10^n) can be changed to n*log2(10) without an overflow in practice. — vinc17, Aug 27 '15 at 21:12
More precisely, the `pow()` will give you an overflow for 309 decimal digits or more, which more or less corresponds to the about 300 in your post. This is a limitation of `double`, not related to MPFR. Do the change I suggested. — vinc17, Aug 27 '15 at 21:40
@JoachimPileborg Thank your for pointing me to these links I will take these into account when writing any further questions on this site. — Emmanuel Mathi-Amorim, Aug 30 '15 at 00:03
@vinc17 Thank you for your advice regarding the power function. I have tested your suggested edit and it works fine. Please write it as an answer and I will accept it. — Emmanuel Mathi-Amorim, Aug 30 '15 at 00:10

score 1 · Accepted Answer · answered Aug 31 '15 at 07:34

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.

MPFR Program Crashes With High Precision

1 Answers1