0

私は最近、円周率を指定された桁数に計算するプログラムを書きました。桁数は、最初のコマンド ライン引数として渡す必要があります。

約 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
4

1 に答える 1

1

問題は、309 桁以上 (投稿の約 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) 形式の方が安全です。評価が速くなるはずです。

于 2015-08-31T07:34:28.523 に答える