Following up on **Orbifold**, I can confirm that Borweins algorithm works quite fast in practice. Years ago when I still had mathematical aspirations, I implemented it in a small C program. If you’re interested I can mail it (or even post it: it is not that large). Here is the ‘descriptive’ part of the program. Even though I didn’t use Fast Fourier Transformations for multiplications it was still reasonably quick on an old XT (in the old days of 5 1/4 disks, when 8 MHz was considered a fast PC). The nice thing was that you could see the number of correct digits double with each iteration.

I’m afraid I don’t even understand my own program anymore.

printf("

This program calculates PI. It uses Borweins Algorithm, and works as follows:

");

printf("Consider three series of numbers, A, B and P.

");

printf("The series are recursive defined as follows:

");

printf("A[0]=û2, B[0]=0, P[0]=2+û2

");

printf("A[k+1]=(ûA[k]+(1/ûA[k]))/2

");

printf("B[k+1]=( ((ûA[k])*(1+B[k])) / (A[k]+B[k]) )

");

printf("P[k+1]=P[k]*B[k+1]*(1+A[k+1])/(1+B[k+1])

");

printf("Now P[ì] = PI. This algorithm is quadratically convergent, which means that

");

printf("the number of correct digits rougly doubles with each iteration.

");

printf("For the calculation of the reciprocal 1/a and the square root ûa two

");

printf("selfcorrecting, quadratically convergent iterations are used:

");

printf("X[k+1]= 2*X[k]-a*X[k]*X[k], converges to 1/a.

");

printf("X[k+1]= (X[k]+(a/X[k]))/2, converges to ûa.

");

printf("For more details on the algorithm, see Math. Comp January 1988, vol.50,

no. 188, p. 283-286, article by D.H.Bailey.

");

printf("The iterations can be ‘seen’ by choosing mode 2. rN is the Nth iteration of

");

printf("the reciprocal, sqN is the Nth iteration of the squareroot. In mode 1 only

");

printf("only A and B are printed besides PI.

");