use strict;
use warnings;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Inline::C
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

use Inline 'C' => Config => ccflagsex => '-O2', inc => '-I/usr/local/include', libs => '-L/usr/local/lib -lgmp -lm', clean_after_build => 0;

use Inline 'C' => <<'END_C';

/* http://beej.us/blog/data/pi-chudnovsky-gmp/
 *
 * Copyright (c) 2012 Brian "Beej Jorgensen" Hall <beej@beej.us>
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

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

// how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555

char *chudnovsky (unsigned long digits)
{
    mpf_t con, A, B, F, sum;
    mpz_t a, b, c, d, e;
    char *output;
    mp_exp_t exp;
    double bits_per_digit;

    unsigned long int k, threek;
    unsigned long iterations = (digits/DIGITS_PER_ITERATION)+1;
    unsigned long precision_bits;

    // roughly compute how many bits of precision we need for
    // this many digit:
    bits_per_digit = 3.32192809488736234789; // log2(10)
    precision_bits = (digits * bits_per_digit) + 1;

    mpf_set_default_prec(precision_bits);

    // allocate GMP variables
    mpf_inits(con, A, B, F, sum, NULL);
    mpz_inits(a, b, c, d, e, NULL);

    mpf_set_ui(sum, 0); // sum already zero at this point, so just FYI

    // first the constant sqrt part
    mpf_sqrt_ui(con, 10005);
    mpf_mul_ui(con, con, 426880);

    // now the fun bit
    for (k = 0; k < iterations; k++) {
        threek = 3*k;

        mpz_fac_ui(a, 6*k);  // (6k)!

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

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

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

        mpz_ui_pow_ui(e, 640320, threek); // -640320^(3k)

        if ((threek&1) == 1) { mpz_neg(e, e); }

        // numerator (in A)
        mpz_mul(a, a, b);
        mpf_set_z(A, a);

        // denominator (in B)
        mpz_mul(c, c, d);
        mpz_mul(c, c, e);
        mpf_set_z(B, c);

        // result
        mpf_div(F, A, B);

        // add on to sum
        mpf_add(sum, sum, F);
    }

    // final calculations (solve for pi)
    mpf_ui_div(sum, 1, sum); // invert result
    mpf_mul(sum, sum, con); // multiply by constant sqrt part

    // get result base-10 in a string:
    output = mpf_get_str(NULL, &exp, 10, digits, sum); // calls malloc()

    // free GMP variables
    mpf_clears(con, A, B, F, sum, NULL);
    mpz_clears(a, b, c, d, e, NULL);

    return output;
}

END_C

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Perl
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

my $pi = chudnovsky(shift || 15);

print "3.", substr($pi, 1), "\n";