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]= 2X[k]-aX[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.
");
