chudnovskyPi.c
char *chudnovsky(unsigned long digits)
{
	mpf_t result, 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(result, 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(result, con, A, B, F, sum, NULL);
	mpz_clears(a, b, c, d, e, NULL);

	return output;
}
/**
 * Compute pi to the specified number of decimal digits using the
 * Chudnovsky Algorithm.
 *
 * http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
 *
 * NOTE: this function returns a malloc()'d string!
 *
 * @param digits number of decimal digits to compute
 *
 * @return a malloc'd string result (with no decimal marker)
 */
Compute pi to the specified number of decimal digits using the
Chudnovsky Algorithm.



this function returns a malloc()'d string!

@param digits number of decimal digits to compute

@return a malloc'd string result (with no decimal marker)
Compute pi to the specified number of decimal digits using the Chudnovsky Algorithm . this function returns a malloc () ' d string! @param digits number of decimal digits to compute @return a malloc ' d string result ( with no decimal marker )
char *chudnovsky(unsigned long digits)
{
	mpf_t result, 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;
	bits_per_digit = 3.32192809488736234789; 
	precision_bits = (digits * bits_per_digit) + 1;
	mpf_set_default_prec(precision_bits);
	mpf_inits(result, con, A, B, F, sum, NULL);
	mpz_inits(a, b, c, d, e, NULL);
	mpf_set_ui(sum, 0); 
	mpf_sqrt_ui(con, 10005);
	mpf_mul_ui(con, con, 426880);
	for (k = 0; k < iterations; k++) {
		threek = 3*k;
		mpz_fac_ui(a, 6*k);  
		mpz_set_ui(b, 545140134); 
		mpz_mul_ui(b, b, k);
		mpz_add_ui(b, b, 13591409);
		mpz_fac_ui(c, threek);  
		mpz_fac_ui(d, k);  
		mpz_pow_ui(d, d, 3);
		mpz_ui_pow_ui(e, 640320, threek); 
		if ((threek&1) == 1) { mpz_neg(e, e); }
		mpz_mul(a, a, b);
		mpf_set_z(A, a);
		mpz_mul(c, c, d);
		mpz_mul(c, c, e);
		mpf_set_z(B, c);
		mpf_div(F, A, B);
		mpf_add(sum, sum, F);
	}
	mpf_ui_div(sum, 1, sum); 
	mpf_mul(sum, sum, con); 
	output = mpf_get_str(NULL, &exp, 10, digits, sum); 
	mpf_clears(result, con, A, B, F, sum, NULL);
	mpz_clears(a, b, c, d, e, NULL);
	return output;
}
char * chudnovsky ( unsigned long digits ) { mpf_t result , 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 ; bits_per_digit = 3.32192809488736234789 ; precision_bits = ( digits * bits_per_digit ) + 1 ; mpf_set_default_prec ( precision_bits ) ; mpf_inits ( result , con , A , B , F , sum , NULL ) ; mpz_inits ( a , b , c , d , e , NULL ) ; mpf_set_ui ( sum , 0 ) ; mpf_sqrt_ui ( con , 10005 ) ; mpf_mul_ui ( con , con , 426880 ) ; for ( k = 0 ; k < iterations ; k ++ ) { threek = 3 * k ; mpz_fac_ui ( a , 6 * k ) ; mpz_set_ui ( b , 545140134 ) ; mpz_mul_ui ( b , b , k ) ; mpz_add_ui ( b , b , 13591409 ) ; mpz_fac_ui ( c , threek ) ; mpz_fac_ui ( d , k ) ; mpz_pow_ui ( d , d , 3 ) ; mpz_ui_pow_ui ( e , 640320 , threek ) ; if ( ( threek & 1 ) == 1 ) { mpz_neg ( e , e ) ; } mpz_mul ( a , a , b ) ; mpf_set_z ( A , a ) ; mpz_mul ( c , c , d ) ; mpz_mul ( c , c , e ) ; mpf_set_z ( B , c ) ; mpf_div ( F , A , B ) ; mpf_add ( sum , sum , F ) ; } mpf_ui_div ( sum , 1 , sum ) ; mpf_mul ( sum , sum , con ) ; output = mpf_get_str ( NULL , & exp , 10 , digits , sum ) ; mpf_clears ( result , con , A , B , F , sum , NULL ) ; mpz_clears ( a , b , c , d , e , NULL ) ; return output ; }
Compute pi to the specified number of decimal digits using the
 Chudnovsky Algorithm.
Compute pi to the specified number of decimal digits using the Chudnovsky Algorithm .
// roughly compute how many bits of precision we need for // this many digit: // log2(10) // allocate GMP variables // sum already zero at this point, so just FYI // first the constant sqrt part // now the fun bit // (6k)! // 13591409 + 545140134k // (3k)! // (k!)^3 // -640320^(3k) // numerator (in A) // denominator (in B) // result // add on to sum // final calculations (solve for pi) // invert result // multiply by constant sqrt part // get result base-10 in a string: // calls malloc() // free GMP variables