Calculating Pi using GMP

Question

I'm trying to calculate Pi using GMP(5.1.0) with this series:

What i have done is:

#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);
}

I'm compiling on 64 bit Linux:

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

Output:

3
3
3
and so on

Expected Output:

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

---

Any help will be greatly appreciated.

Answer 1

You do integer division:

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

Dividing two integers will round them towards zero, which is 0 in this case, since a * b * c > 4 in every iteration.

You could fix this by writing

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

However, you should perform the division using GMP since the code above suffers from the limits of the native number types. Also, the result of this division shouldn't be stored in a GMP integer, but in a GMP float:

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

Note that also a, b, c have to be GMP floats in order to make this work.