Программа MPFR падает с высокой точностью

Question

Недавно я написал программу для расчета числа Пи до заданного количества цифр. Количество цифр должно быть передано в качестве первого аргумента командной строки.

При запуске с числовым значением менее 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

Answer 1

Проблема связана с тем, что для 309 цифр или более (что более или менее соответствует примерно 300 в вашем посте), pow (10.0, digits) переполняет, давая бесконечность, так что ceil(log2(pow(10.0, (double)digits))) также дает бесконечность, и при преобразовании в целое число это имеет неопределенное поведение. Переполнение здесь является ограничением double тип, не относящийся к MPFR. Чтобы избежать огромных чисел по величине (таким образом, переполнения), математическое выражение log2(10^n) можно заменить на n*log2(10), который на практике не будет переполнен.

Примечание: MPFR имеет гораздо больший диапазон экспонент, чем double, так что вычисление log2(10^n) в MPFR вместо double также позволит избежать переполнения при разумных значениях n, но форма n * log2 (10) в любом случае безопаснее и должна быть быстрее оценена.