Pi-greco – Chudnovsky :: C

While being a big fan of Python I’m also in the continuous pursuit of better performance.
Like it or not, but when seeking for better performance in many cases you have to put away high level framework/programming-languages and go back to basic; many times basic means: C/C++.
After had developed some different Python scripts to calculate the Pi-greco up to a given number of digits (later I will post the most advanced version I came up with) I decided it was time to make a step forward and write it in C.
Because I’m still away from home this time the test machine instead of being the usual FX-8120 is my Lenovo Thinkpad equipped with a 2.1 GHz 2C4T 3 MB of cache L3 i3-2310M.
On Windows the most known C (it also have a C++ wrapper) library for arbitrary precision math is MPIR.
One of the most important “features” of the library is that it’s written in both C and Assembly, also they provide different versions of the Assembly code, so you can use the one that is specifically written for the CPU architecture of the machine you will run your program on.
The full source code of the library can be downloaded, modified and compiled by the user to better suit his needs.
The only downside is that the super optimized Assembly code included in the latest stable release of the library isn’t really up to date, in fact there still no support for Sandy Bridge CPUs and AMD Bulldozer ones.
Anyway I think none should complain about it, MPIR is an open source project so we can wait for official release or get our hands dirty and write the code we need.
I’m lazy, so I will wait for the official team to realease the updated Assembly; since then I can live with the current stuff.
Talking about “current stuff” I’ve noticed that AMD K10 library’s version is much faster than P4, C2D and Nehalem ones even if I’m using an Intel Sandy Bridge CPU; and I’m not talking about 2/3 % faster but 15/20 %.
Anyway, while the Python – I’ve the feeling that there still something left here, tho I don’t think I can make it as fast as the C implementation is – script (which uses GMPY2) requires 30+ seconds to compute one million of digits the program written in C takes only 1.4 seconds.
Both of the programs are written with the same algorithm: Chudnovsky brothers formula + binary splitting to further speed up the computation.
Here is the part of the code deputed to do the math:

//------------------------------------------------------------------//
// pi_chudnovsky_bs
// description: compute PI to the given precision
void pi_chudnovsky_bs(mp_bitcnt_t digits, char* cpuid)
{
	reset();
	double start = seconds();

	mpf_set_default_prec (prec);

	__int64 C3 = 262537412640768000;

	printf ("Computing DIGITS_PER_TERM: ");
	//double C3_OVER_24 = 10939058860032000;
	mpz_t C3_OVER_24_p;
	mpz_init_set_str (C3_OVER_24_p, "10939058860032000", 10);
	double DIGITS_PER_TERM = log10 ((double) 151931373056000);
	double end = seconds();
	double end_2 = end - start;
	printf ("%g seconds\n", end_2);

	printf ("Computing Binary Splitting parameters: ");
	unsigned long N = (digits / DIGITS_PER_TERM) + 1;
	end = seconds();
	end_2 = end - start;
	printf ("%g seconds\n", end_2);

	printf ("Computing Binary Splitting terms: ");

	mpz_t Pab;
	mpz_init (Pab);
	mpz_t Qab;
	mpz_init (Qab);
	mpz_t Tab;
	mpz_init (Tab);

	bs(0, N, C3_OVER_24_p, &Pab, &Qab, &Tab);

	end = seconds();
	end_2 = end - start;
	printf ("%g seconds\n", end_2);

	// one_squared = 10**(2*digits)
	printf ("Computing squaring: ");
	mpz_t one_squared;
	mpz_init_set_ui (one_squared, 10);
	mpz_pow_ui (one_squared, one_squared, (2 * digits));

	end = seconds();
	end_2 = end - start;
	printf ("%g seconds\n", end_2);

	// sqrtC = sqrt(10005 * one_squared)
	printf ("Computing square root: ");
	mpz_t sqrtC;
	mpz_init (sqrtC);
	mpz_mul_ui (sqrtC, one_squared, 10005);
	mpz_sqrt (sqrtC, sqrtC);

	end = seconds();
	end_2 = end - start;
	printf ("%g seconds\n", end_2);

	// pi = (Q*426880*sqrtC) // T
	printf ("Computing PI-Greco");
	mpz_t pi;
	mpz_init (pi);
	mpz_mul_ui (pi, Qab, 426880);
	mpz_mul (pi, pi, sqrtC);
	mpz_div (pi, pi, Tab);

	// count end time
	end = seconds();
	end_2 = end - start;
	printf ("\nCompleted in: %g seconds\n", end_2);
	//printf ("\nPI: ");
	//mpf_out_str (stdout, 10, 0, pi);

	// delete temorary variables
	mpz_clear (Pab);
	mpz_clear (Qab);
	mpz_clear (Tab);

	printf ("\n");
	writeResultToFile (pi, cpuid, digits);
	printf ("\n");
	writeInfoToFile (digits, cpuid, end);
}

//------------------------------------------------------------------//
// bs
// description: makes the binary splitting math
void bs (unsigned long a_i, unsigned long b_i, mpz_t C3_OVER_24,
                             mpz_t* Pab, mpz_t* Qab, mpz_t* Tab)
{
	mpf_set_default_prec (prec);

	mpz_t a;
	mpz_init_set_ui (a, a_i);
	mpz_t b;
	mpz_init_set_ui (b, b_i);

	if ((b_i - a_i) == 1)
	{
		if (a_i == 0)
		{
			// Pab = Qab = ONE
			mpz_set_ui(*Pab, 1);
			mpz_set_ui(*Qab, 1);
		} else {
			// Pab = (6*a-5)*(2*a-1)*(6*a-1);
			mpz_t p1;
			mpz_init (p1);
			mpz_mul_ui (p1, a, 6);
			mpz_sub_ui (p1, p1, 5);

			mpz_t p2;
			mpz_init (p2);
			mpz_mul_ui (p2, a, 2);
			mpz_sub_ui (p2, p2, 1);

			mpz_t p3;
			mpz_init (p3);
			mpz_mul_ui (p3, a, 6);
			mpz_sub_ui (p3, p3, 1);

			mpz_mul (*Pab, p1, p2);
			mpz_mul (*Pab, *Pab, p3);
			// delete temorary variables
			mpz_clear (p1);
			mpz_clear (p2);
			mpz_clear (p3);

			// Qab = a^3 * C3_OVER_24
			mpz_t a3;
			mpz_init (a3);
			mpz_pow_ui (a3, a, 3);

			mpz_mul (*Qab, a3, C3_OVER_24);
			// delete temorary variables
			mpz_clear (a3);
		}
		// Tab = Pab * (13591409 + 545140134*a)
		mpz_t t1;
		mpz_init (t1);
		mpz_mul_ui (t1, a, 545140134);
		mpz_add_ui (t1, t1, 13591409);
		mpz_mul (*Tab, *Pab, t1);
		// delete temorary variables
		mpz_clear (t1);

		if (a_i & 1)
			mpz_neg (*Tab, *Tab);
	} else {
		unsigned long m = (a_i + b_i) / 2;

		// divide et impera, first chunk
		mpz_t Pam;
		mpz_init (Pam);
		mpz_t Qam;
		mpz_init (Qam);
		mpz_t Tam;
		mpz_init (Tam);

		bs(a_i, m, C3_OVER_24, &Pam, &Qam, &Tam);

		// divide et impera, second chunk
		mpz_t Pmb;
		mpz_init (Pmb);
		mpz_t Qmb;
		mpz_init (Qmb);
		mpz_t Tmb;
		mpz_init (Tmb);

		bs(m, b_i, C3_OVER_24, &Pmb, &Qmb, &Tmb);

		mpz_mul (*Pab, Pam, Pmb);
		mpz_mul (*Qab, Qam, Qmb);

		mpz_init (*Tab);
		mpz_t T1;
		mpz_init (T1);
		mpz_t T2;
		mpz_init (T2);
		mpz_mul (T1, Qmb, Tam);
		mpz_mul (T2, Pam, Tmb);
		mpz_add (*Tab, T1, T2);

		// delete temorary variables
		mpz_clear (Pam);
		mpz_clear (Pmb);
		mpz_clear (Qam);
		mpz_clear (Qmb);
		mpz_clear (Tam);
		mpz_clear (Tmb);
		mpz_clear (T1);
		mpz_clear (T2);
	}

	mpz_clear (a);
	mpz_clear (b);
}

I know that at a first glance it will looks a big mess, I hope the comments I spread in the code would be useful to understand what the program really does.
If you are interested here is the link to the whole source code (already compiled K10’s MPIR library are included): source code.
Also here’s the executable version (Win X64 AMD K10): exe.

Leave a Comment

Your email address will not be published. Required fields are marked *