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

#define PI_DEBUG
/*
*/

int main() {
	int n, iterations=10;
	mpf_t y, a, ynew, anew, t1, t2, one;
	// mpf_set_default_prec(8*1000000+30);
	//mpf_set_default_prec(594100);
	mpf_set_default_prec(4*2376240+100);
#ifdef PI_DEBUG
	fprintf(stderr,"Set default precision\n");
#endif
	mpf_init(y);
	mpf_init(a);
	mpf_init(ynew);
	mpf_init(anew);
	mpf_init(t1);
	mpf_init(t2);

#ifdef PI_DEBUG
	fprintf(stderr,"Initialized variables\n");
#endif
	mpf_sqrt_ui(t1,2);
	mpf_sub_ui(y,t1,1);
	mpf_mul_ui(t2,t1,4);
	mpf_ui_sub(a,6,t2);
	
	mpf_mul(t1,y,y);
#ifdef PI_DEBUG
	fprintf(stderr,"Initial values set.\n");
#endif
	
	for (n=0;n<iterations;n++) {
		fprintf(stderr, "Starting iteration %d\n", n+1);
#ifdef PI_DEBUG
		//mpf_ui_div(t2,1,a);
		//fprintf(stderr,"pi_%d = ",n);
		//mpf_out_str(stderr, 10, 0, t2);
		//putc('\n',stderr);
#endif
		mpf_mul(t2,t1,t1);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished yn^4\n");
#endif
		mpf_ui_sub(t1,1,t2);
		/* Fourth root of 1-y_n^4*/
		mpf_sqrt(t2,t1);
		mpf_sqrt(t1,t2);
		
#ifdef PI_DEBUG
		fprintf(stderr,"Finished square roots\n");
#endif
		/* use anew for the moment as a temporary variable */
		mpf_ui_sub(anew,1,t1);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished 1-t1 \n");
#endif
		mpf_add_ui(t2,t1,1);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished 1+t1 \n");
#endif
		mpf_div(ynew,anew,t2);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished division\n");
#endif

		/* Work on anew */
		/* use y and anew as a temporary variable */
		mpf_add_ui(anew,ynew,1);
		mpf_pow_ui(y,anew,4);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished first power\n");
#endif
		mpf_mul(t2,y,a);
		/* t1=ynew^2 which will be used at the start of the next iteration */
		mpf_mul(t1,ynew,ynew);
		/* now use a as a temporary variable */
		mpf_add(a,t1,anew);
		mpf_mul(y,ynew,a);
		mpf_set_ui(a,2);
		mpf_pow_ui(anew,a,2*n+3);
#ifdef PI_DEBUG
		fprintf(stderr,"Finished second power\n");
#endif
		mpf_mul(a, anew, y);
		mpf_sub(anew, t2, a);

		if (n+1 < iterations) {
			mpf_set(a,anew);
			mpf_set(y,ynew);
		}
	}
	fprintf(stderr,"Finished with the loop.\n");
	mpf_ui_div(t2,1,anew);
	fprintf(stderr,"Finished inverting.\n");
	mpf_out_str(stdout, 10, 0, t2);
	putc('\n',stdout);
	fprintf(stderr,"Output finished.\n");
}