/*
 * The Chudnovsky formula that was used is from:
 * http://www.craig-wood.com/nick/articles/pi-chudnovsky/
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <pthread.h>


int main(int argc, char *argv[]){
	setpriority(0, 0, 19);
	mpf_t  a_k, a_sum, b_sum, C, C3;
	mpf_t k, tmp1, tmp2, tmp3, n1, total, pi;
	mpf_t one;
	int i;

	//get enough bits to put 1million digits about 4mil bits
	mpf_set_default_prec(3330000);

	//inii all the gmp vars
	mpf_init(k);
	mpf_init(a_k);
	mpf_init(a_sum);
	mpf_init(b_sum);
	mpf_init(C);
	mpf_init(C3);
	mpf_init(tmp1);
	mpf_init(tmp2);
	mpf_init(tmp3);
	mpf_init(total);
	mpf_init(pi);
	mpf_init(n1); //negative 1
	mpf_init(one);


	//declare the values for the vars (first iteration of chudnovsky) k = 0
	mpf_set_str(k, "1", 10);
	mpf_set_str(a_k, "1000", 10);
	mpf_set_str(a_sum, "1000", 10);
	mpf_set_str(b_sum, "0", 10);
	mpf_set_str(C, "640320", 10);
	mpf_set_str(n1, "-1", 10);

	//C3 = C*C*C /24
	mpf_pow_ui(C3, C, (unsigned long int) 3);
	mpf_div_ui(C3, C3, (unsigned long int) 24);
	
	i = 0;
	while(1){
		i++;
		//tmp1 = (6*k-5)*(2*k)*(6*k-1)
		mpf_mul_ui(tmp1, k, (unsigned long int) 6); 
		mpf_mul_ui(tmp2, k, (unsigned long int) 2); 
		mpf_sub_ui(tmp3, tmp1, (unsigned long int) 1);
		mpf_sub_ui(tmp2, tmp2, (unsigned long int) 1);
		mpf_sub_ui(tmp1, tmp1, (unsigned long int) 5);
		mpf_mul(tmp1, tmp1, tmp2);
		mpf_mul(tmp1, tmp1, tmp3);
		mpf_mul(tmp1, tmp1, n1);

		//a_k = tmp1
		mpf_mul(a_k, a_k, tmp1);

		//tmp1 = k*k*k*c3
		mpf_pow_ui(tmp1, k, (unsigned long int) 3);
		mpf_mul(tmp1, tmp1, C3);

		//a_k = a_k / tmp1
		mpf_div(a_k, a_k, tmp1);
		
		//a_sum = a_sum + a_k
		mpf_add(a_sum, a_sum, a_k);

		//b_sum = b_sum + (k*a_k);
		mpf_mul(tmp1, k, a_k);
		mpf_add(b_sum, b_sum, tmp1);

		mpf_add_ui(k, k, (unsigned long int) 1);
		
		//Iterate over Chudnovsky fomrula upto k = i = 70514
		if(i == 70514) break;
	}

	//total = 13591409*a_sum + 545140134*b_sum;
	mpf_mul_ui(a_sum, a_sum, (unsigned long int) 13591409);
	mpf_mul_ui(b_sum, b_sum, (unsigned long int) 545140134);
	mpf_add(total, a_sum, b_sum);

	// pi = (426880*sqrt(10005) / total;
	mpf_set_str(tmp1, "10005", 10);
	mpf_sqrt(tmp1, tmp1);
	mpf_mul_ui(tmp1, tmp1, (unsigned long int) 426880);
	mpf_div(pi, tmp1, total);

	//print pi
	gmp_printf("%.Ff\n", pi);

	//clear the vars used
	mpf_clear(k);
	mpf_clear(a_k);
	mpf_clear(a_sum);
	mpf_clear(C);
	mpf_clear(tmp1);
	mpf_clear(tmp2);
	mpf_clear(tmp3);
	mpf_clear(total);
	mpf_clear(pi);
	mpf_clear(n1);
	mpf_clear(one);
}