To the spigot programs, I have several comments .
At spigotpi.c
The program is a 1:1 port of the PASCAL program of the original paper of
Rabinowitz and Wagon.
The ported program has herited several bugs from the
original program and this in turn suffers from the erreanous algorithm:
1. The program fails when used to calculate N=1 or N=32 decimals
of pi. In these cases is the chain length LEN = [10N/3] by 1 too small.
Ironically enough, Rabinowitz and Wagon "proof" in their paper
the expression [10N/3] to be "correct" (Lemma 5).
2. The program fails when used to calculate N=85,167,245 and
many more values of N.
The reason is that the number of iterations of the outer loop is too
small. This loop must be iterated "as long as the number of true
decimals is not greater N".
3. The program computes in many cases less decimals than it should.
Such cases are e.g. N= 6,13,15,31 and many more. An extreme is N=768 where
only 761 digits will be computed.
This occurs for the same reason as point 2.
4. Here is a version of the program which seems to be correct in these
respects and hopefully without new bugs. (It is also faster.)
------------------------------------------------------------------
/* A spigot algorithm for the digits of pi, Stanley Rabinowitz
and Stan Wagon, Amer.Math.Monthly, March 1995, 195-203
Bug fixes by C. Haenel.
The program can be compiled with an 16- or 32-Bit compiler.
The upper limit for N is 9826 on 16 Bit and (at least) 50000
on 32 Bit compilers.
*/
#define N 1000 // Decimals to compute.
#define LEN (10L * N) / 3 + 1 // Chain length
unsigned j, predigit, nines, a[LEN];
long x, q, k, len, i;
void main()
{
for(j=N; j; )
{
q = 0;
k = LEN+LEN-1;
for(i=LEN; i; --i)
{
x = (j == N ? 20 : 10L*a[i-1]) + q*i;
q = x / k;
a[i-1] = (unsigned)(x-q*k);
k -= 2;
}
k = x % 10;
if (k==9) ++nines;
else
{
if (j)
{
--j;
printf("%ld", predigit+x/10);
}
for(; nines; --nines)
{
if (j) --j, printf(x >= 10 ? "0" : "9");
}
predigit = (unsigned)k;
}
}
}
------------------------------------------------------------------
Ad "shortest program for pi you have ever seen":
1. The program fails when used to calculate 360 decimals. The
size of array f is in that case by one too small.
2. Here is a version which seems to be correct
in this respect and which is hopefully without new bugs:
(It is also faster and shorter.)
Version "Fast": (approx. 25% faster)
/* Calculation of pi to 32372 decimal digits */
/* Size of program: 152 characters */
/* After Dik T. Winter, CWI Amsterdam */
unsigned a=1e4,b,c=113316,d,e,f[113316],g,h,i;
main(){for(;b=c,c-=14;i=printf("%04d",e+d/a),e=d%a)
while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d-g*h;}
I wonder whether such artists as the one who made "xmas.c" can
help to squeeze Winter's program even more.
Version "Short":
/* Calculation of pi to 16276 decimal digits */
/* Size of program: 143 characters */
/* After Dik T. Winter, CWI Amsterdam */
int a=1e4,b,c=56980,d,e,f[56980],g,h,i;
main(){for(;b=c,c-=14;i=printf("%04d",e+d/a),e=d%a)
while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d%g;}
---------------------------------------------------------
Christoph Haenel
chrhaenel@aol.com