# The problem with Pi (ii)

To print out precise numbers, we need another approach. One such approach to π was offered in 1990, when  Rabinowitz and Wagon discovered a “spigot” algorithm for π [1]. This algorithm allows for successive digits of π to be calculated using a simple recursive algorithm.

Here is the Pascal program from the original paper:

```program pi_spigot;
const n = 1000;
len = 10*n div 3;
var i, j, k, q, x, nines, predigit : integer;
a : array[1..len] of longint;
begin
for j := 1 to len do a[j] := 2;
nines := 0;
predigit := 0;
for j := 1 to n do
begin
q := 0;
for i := len downto 1 do
begin
x := 10*a[i] + q*i;
a[i] := x mod (2*i-1);
q := x div (2*i-1);
end;
a[1] := q mod 10;
q := q div 10;
if q = 9 then nines := nines + 1
else if q = 10 then
begin
write(predigit+1);
for k := 1 to nines do write(0);
predigit := 0;
nines := 0;
end
else begin
write(predigit);
predigit := q;
if nines <> 0 then
begin
for k := 1 to nines do write(9);
nines := 0;
end
end
end;
writeln(predigit);
end.
```

The reader is left to browse the original paper for algorithm details. When running the code (I used the Free Pascal Compiler), it gives the following output:

03141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923459099207118856806516745194252854340964271825877594867237227818052636476240001253966769588037865014608282022810293731532093125719608182508805545474738567199841858637460939653266515205973230200470327367748014549418069976867613243271975438847522856914469354541235658607836752109748167637744042444357591344102954987267348532155345226074874607261272294756681440117171146164262329464280245194165531036682264959997624747949649719425809626723943536390912536760246664166117527686634968595119533790744923019073345333321686533168484290724097901531077274052175927055788110268222469183375946146833893295768948223454910618813165693638675020776147309271039650546608291414446931319129743860775232869818000001373465107828774083538328579350606183962012466

Only the first 265 digits are correct. This is probably due to overflow errors due to datatypes being too small. Convert the program to C, and it works as it should:

```
#include <stdio.h>
#include <stdlib.h>

void spigot();

int main(void) {
spigot();
return 0;
}

void spigot() {
int i, j, k, q, x;
int len, nines=0, predigit=0;
int N=1000;

len = (10*N/3)+1;
int a[len];

// Initialize A to (2,2,2,2,2,...,2)
for (i=0; i<len; i=i+1)
a[i] = 2;

// Repeat n times
for (j=1; j<=1000; j=j+1) {
q = 0;
for (i=len; i>0; i=i-1) {
x = 10 * a[i-1] + q*i;
a[i-1] = x % (2*i-1);
q = x / (2*i-1);
}
a[0] = q % 10;
q = q / 10;
if (q == 9)
nines = nines + 1;
else if (q == 10) {
printf("%d", predigit+1);
for (k=0; k<nines; k=k+1)
printf("%d",0);
predigit = 0;
nines = 0;
}
else {
printf("%d", predigit);
predigit = q;
if (nines != 0) {
for (k=0; k<nines; k=k+1)
printf("%d",9);
nines = 0;
}
}
}
printf("%d\n", predigit);
}
```

[1] ] Rabinowitz, S.D., Wagon, S., ” A spigot algorithm for pi”. American Mathematical Monthly, 102, pp.195–203 (1995)

This site uses Akismet to reduce spam. Learn how your comment data is processed.