c - GMPを使用した円周率の計算

Question

このシリーズでGMP（5.1.0）を 使用して円周率を計算しようとしています。

私がしたことは：

#include <stdio.h>
#include <gmp.h>

mpz_t pi;
mpz_t next; // Every number in the serie that comes after 3

int a = 2, b = 3, c = 4; // Divisors. Everytime add 2 to all of them at the end of the loop
int i; // Loop counter

int main (void) {
    mpz_init (pi); mpz_init (next);
    mpz_set_str (pi, "3", 10);

    for (i = 2 ;; ++i) {
        printf ("%s\n", mpz_get_str (NULL, 10, pi));

        mpz_set_d (next, 4 / (a * b * c));

        if (i % 2 == 0)
            mpz_add (pi, pi, next);
        else
            mpz_sub (pi, pi, next);

        a += 2; b += 2; c += 2;
    }

    mpz_clear (next);
    mpz_clear (pi);
}

私は64ビットLinuxでコンパイルしています：

gcc -Wall -pedantic -o foo foo.c -lgmp

出力：

3
3
3
and so on

期待される出力：

3
3.1666...
3.1333...
3.1452...
3.1396...
and so on

---

どんな助けでも大歓迎です。

Answer 1

あなたの問題はこの行にある可能性が最も高いです：

mpz_set_d (next, 4 / (a * b * c));

この行では、基本的にmpz_set_d関数に0を渡しています。これにより、数値が設定されると思います。

4 / (a * b * c)これは、整数演算を使用して評価しているためです。(a * b * c)は常により大きいので4、この式は常に0と評価されます。これはあなたが望むものではありません。

おそらく、2つの別々のGMP浮動小数点数に入れ4、(a * b * c)GMP関数を使用して除算を実行し、結果をnext変数として扱う必要があります。

編集： GMPドキュメントを見るとmpz、関数ファミリーは整数を扱います。おそらくmpf、浮動小数点数を処理するmpq関数、または有理数を処理する関数を使用する必要があります。

Answer 2

整数除算を行います。

mpz_set_d (next, 4 / (a * b * c));
//               ^   ^^^^^^^^^^^
//              int      int

2つの整数を除算すると、ゼロに向かって丸められます。0この場合は、a * b * c > 4すべての反復で発生するためです。

あなたは書くことによってこれを修正することができます

mpz_set_d (next, 4.0 / (a * b * c));
//               ^^^   ^^^^^^^^^^^
//             double      int

ただし、上記のコードにはネイティブの数値タイプの制限があるため、GMPを使用して除算を実行する必要があります。また、この除算の結果はGMP整数ではなく、GMPフロートに保存する必要があります。

mpf_t next;
//...
mpf_set_d(next, 4.0);
mpf_div(next, a);
mpf_div(next, b);
mpf_div(next, c);

これを機能させるには、a、b、cもGMPフロートである必要があることに注意してください。

Answer 3

4 / (a * b * c)常にゼロです。