display | more...
	c = sqrt (0.125);
	a = 1.0 + 3.0 * c;
	b = sqrt (a);
	e = b - 0.625;
	b *= 2.0;
	c = e - c;
	a += e;
	npow = 4.0;

	do
	{
		npow *= 2.0;
		e = (a + b) / 2.0;
		b = sqrt (a * b);
		e -= b;
		b *= 2.0;
		c -= e;
		a = e + b;
	}
	while (e > SQRT_SQRT_EPSILON);

	e = e * e / 4.0;
	a += b;
	pi = (a * a - e - e / 2.0) / (a * c - e) / npow;

This is a modified version of the Gauss-Legendre formula by Takuya Ooura.

a.k.a. the AGM

The main advantage of using this algorithm (in addition to being simple) is that it has squared convergence. This means the number of accurate digits of pi doubles every iteration.

This algorithm was one of the two algorithms used to compute (and check) the world record of 206,158,430,000 digits of pi (π) by Yasumasa Kanada and Daisuke Takahashi on September 20th, 1999. The other one was Borweins' 4-th order convergent algorithm.

Reference:
http://www.lupi.ch/PiSites/Pi-Rekord.html