Calculating Pi using GMP [closed]

Question

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.

Closed 9 years ago.

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.

Answer 2

Your problem is most likely on this line:

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

On this line you're essentially passing 0 to the mpz_set_d function, which I assume sets your number.

That's because you are evaluating 4 / (a * b * c) using integer arithmetic. (a * b * c) is always greater than 4, so this expression always evaluates to 0. This is not what you want.

You should probably put 4 and (a * b * c) into two separate GMP floating point numbers and use the GMP functions to do the divison, the treat the result as your next variable.

EDIT: Looking at GMP documentation, the mpz family of functions deals with integers. You probably need to use the mpf functions, that deal with floats, or the mpq functions, that deal with rationals.

Answer 3

4 / (a * b * c) is always zero.