pi.c

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

char* calculatePi(unsigned long digits);

int main()
{
    const unsigned long piDigits = 500;
    char* pi = calculatePi(piDigits);
    printf("%.1s.%s\n", pi, pi + 1); // Since there is no decimal marker in our string, we'll print one
    free(pi);

    return 0;
}

/**
* Calculates pi to the specified number of digits
* using the Chudnovsky brothers' algorithm
*
* Function returns a malloc()'ed string
*
* @param digits Number of pi digits to compute
* @return       malloc()'ed string result with no decimal marker
**/

char* calculatePi(unsigned long digits)
{
    mpz_t a, b, c, d, e;
    mpf_t numerator, denominator, fraction, sum, A, result;

    mp_exp_t exp;

    const double bitsPerDigit = 3.32192809488736234789;    // precision = digits * log2(10)
    unsigned long precision = (digits * bitsPerDigit) + 1; // +1 for safety
    mpf_set_default_prec(precision);

    const double digitsPerIteration = 14.1816474627254776555;
    unsigned long iterations = (digits / digitsPerIteration) + 1; // +1 for safety

    mpz_inits(a, b, c, d, e, NULL);
    mpf_inits(numerator, denominator, fraction, sum, A, result, NULL);

    for(unsigned long k = 0; k < iterations; k++)
    {
        mpz_fac_ui(a, 6*k); // (6k)!

        mpz_set_ui(b, 545140134);   //
        mpz_mul_ui(b, b, k);        // 13591409 + 545140134k
        mpz_add_ui(b, b, 13591409); //

        mpz_mul(a, a, b);        // (6k!)(13591409 + 545140134k)
        mpf_set_z(numerator, a); //

        mpz_fac_ui(c, 3*k); // (3k)!

        mpz_fac_ui(d, k); // (k!)^3
        mpz_pow_ui(d, d, 3); //

        mpz_ui_pow_ui(e, 640320, 3*k); //
        if (3*k % 2 != 0)              // (-640320)^(3k)
            mpz_neg(e, e);             //

        mpz_mul(c, c, d);          //
        mpz_mul(c, c, e);          // (3k!)(k!)^3(-640320)^3k
        mpf_set_z(denominator, c); //

        mpf_div(fraction, numerator, denominator);
        mpf_add(sum, sum, fraction);
    }

    mpf_ui_div(sum, 1, sum);  // inverted sum

    mpf_sqrt_ui(A, 10005);    // 426880 * sqrt(10005)
    mpf_mul_ui(A, A, 426880); //

    mpf_mul(result, sum, A);

    const int base = 10;
    char* output = mpf_get_str(NULL, &exp, base, digits, result);

    mpf_clears(numerator, denominator, fraction, sum, A, result, NULL);
    mpz_clears(a, b, c, d, e, NULL);

    return output;
}