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 位或更多位（或多或少对应于 post 中的大约 300 位），pow (10.0, digits) 溢出，给出无穷大，因此ceil(log2(pow(10.0, (double)digits))) 也给出无穷大，当转换为整数时，这具有未定义的行为。这里溢出是double类型的限制，与MPFR无关。为避免数量级过大（因此溢出），可以将数学表达式 log2(10^n) 替换为 n*log2(10)，这在实践中不会溢出。

注意：MPFR 的指数范围比 double 大得多，因此在 MPFR 中计算 log2(10^n) 而不是 double 也可以避免 n 的合理值溢出，但是 n*log2(10) 形式在任何情况下都更安全并且应该更快地评估。

MPFR 程序以高精度崩溃

MPFR Program Crashes With High Precision

c

arbitrary-precision

mpfr