Stanley Rabinowitz
The spigot algorithm
A. Sale  D. Saada  S. Rabinowitz
The Principle
So, what is this algorithm?
Well, let's remember that we know pi in its decimal form, which is in base 10:
Pi=3.141592653589793238462643383279...
this can also be written as :
This is known as Horner's form; it allows us, among
other things, to reduce the number of multiplications when evaluating a polynomial.
We can see here that it is base 10, and that the base's step is constant, meaning that we
have 10 at each decimal place (digit if base different from 10). We note this base [1/10,1/10,1/10...]
Now, let's consider the series discovered by Euler (See page about him for the demonstration)
and let's write it in Horner's form:
Now this is interesting, if we compare it with the previous expression, we can see
that we are considering a base with a changing step [1,1/3,2/5,3/7,4/9...]. And in this base, Pi is [2;2,2,2,2...]. So in this base, Pi is one of the simpliest numbers that exists !
We know Pi 's digits in this base, so to compute Pi 's decimal places in base 10 one by one, one just needs to build an
algorithm that changes it to base 10, which is precisely the principle of the spigot algorithm.
A historical overview of this method
I might as well tell you that I haven't found much
information about the above mathematicians! We can't always find what we want on
the web, and if someone doesn't have a personal page, it limits what one can find
out about him.
Initially, it was A. Sale who had the idea of this method in 1968 and used it to
compute e.
Then D. Saada suggested applying it to Pi in 1988 and so did S. Rabinowitz in 1991. The latter is quite well known, he's a confirmed hacker,
and even a mathematician who has published quite a few articles (He is a Mac friend
by the way). He studied some of the finesses of this algorithm with his collegue
Stan Wagon in 1995 in an article in the Mathematical American Monthly.
Around :
First of all, a little computation concerning how
much memory it will take. In Horner's form, we can see that the n/(2n+1) step of the
base is slightly under 1/2 every time. The exact value 1/2 would bring us to consider base 2. For the conversion from base 2 to base 10, we will need roughly Log_{2}(10^{n})=3,32n digits. We can consider it is about the same value for
our base with a changing step to base 10. So if we want to get n décimal places, we will have to consider 3.32n boxes. To get 4 décimal places (including the 3 before the decimal point), we build a memory table of about 14 columns. In fact, 13 will be enough here. The translation of this algorithm would bring
us to consider a table like this:
A 

r 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
B 
= 

3 
5 
7 
9 
11 
13 
15 
17 
19 
21 
23 
25 
Initialization 

2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 















*10 

20 
20 
20 
20 
20 
20 
20 
20 
20 
20 
20 
20 
20 
carried over 

10 
12 
12 
12 
10 
12 
7 
8 
9 
0 
0 
0 
0 
sum 
3 
30 
32 
32 
32 
30 
32 
27 
28 
29 
20 
20 
20 
20 
remainder 

0 
2 
2 
4 
3 
10 
1 
13 
12 
1 
20 
20 
20 















*10 

0 
20 
20 
40 
30 
100 
10 
130 
120 
10 
200 
200 
200 
carried over 

13 
20 
33 
40 
65 
48 
98 
88 
72 
150 
132 
96 
0 
sum 
1 
13 
40 
53 
80 
95 
148 
108 
218 
192 
160 
332 
296 
200 
remainder 

3 
1 
3 
3 
5 
5 
4 
8 
5 
8 
17 
20 
0 















*10 

30 
10 
30 
30 
50 
50 
40 
80 
50 
80 
170 
200 
0 
carried over 

11 
24 
30 
40 
40 
42 
63 
64 
90 
120 
88 
0 
0 
sum 
4 
41 
34 
60 
70 
90 
92 
103 
144 
140 
200 
258 
200 
0 
remainder 

1 
1 
0 
0 
0 
4 
12 
9 
4 
10 
6 
16 
0 















*10 

10 
10 
0 
0 
0 
40 
120 
90 
40 
100 
60 
160 
0 
carried over 

4 
2 
9 
24 
55 
84 
63 
48 
72 
60 
66 
0 
0 
sum 
1 
14 
12 
9 
24 
55 
124 
183 
138 
112 
160 
126 
160 
0 
remainder 

4 
0 
4 
3 
1 
3 
1 
3 
10 
8 
0 
22 
0 
And now let's explain how and why it works !
We recognise in the two first lines the numerator and the denominator of the steps
of the changing step base. In the third line, we find Pi 's expression in that base. And we fill the last column
of the remainder line with 0 's.
The algorithm to convert from right to left in the table is as follows :
Formally, the general principle consists in multiplying the digit by 10 and in computing
the remainder in an equivalent of a Euclidian division of this number by the step
of the changing step base. A number carried over will then appear at each calculation
and comes from the same phenomenon which appears when one multiplies 53 by 8; firstly
one multiplies 3 by 8, then one carries over 2 and adds it to the multiplication of 8 by 5 which is the multiplication
of the next digit.
Filling in the red column for example :
1) First we fill in the *10 line by multiplying the previous line by 10. No problem so far
!
2) We place in sum the sum of the *10 line with the carried over line, which is:
20+9=29
3) We do the Euclidian division of the sum by the
number in line B of the same column :
29=17*1+12
4) We place the 12 remainder in the remainder line. Then we multiply the 1 quotient
with line A and put the result in the next column's carried over line :
1*8=8
The 9 remainder of the red column comes from the previous column's
computation. And so here we pass on the 8 remainder to the next sum's computation.
By repeating the process with the other columns, we fill in the first table and we
get 30 as the last sum. But as we multiplied everything by 10, we take 3 as the first digit of Pi. That's it!
One little remark however : We can consider that a number in the last column is right
if it is not followed by a 9.
When the remainder in the r column is higher than 100, we can obtain (but it is very rare...) 10 after the last decimal
place that we take for Pi. We then take this digit plus 1 as digit for Pi, It's as simple as that!
Other formulas ?
Well I think any rational series should be all
right as long as in it's Horner form, we only find integrers as expression of Pi in the changing
step base (for the previous series, we only had 2's). Particularly, Ramanujan rational series should be able
to be used and give a better algorithm.
Gosper's page's formula is also valid :
As the fraction of two terms of the series is always smaller than 1/13, we will need Log_{13}(10^{n})=0.9 boxes to get a digit. So if we consider 6 boxes, we can hope
to get 6*1/0.9=6.6 decimal places, so 6 or even 7 decimal places at best. We can see that this is respected
perfectly in the table because we even obtain 7 décimal places :
A 

r 
1 
6 
15 
28 
45 
B 
= 

60 
168 
330 
546 
816 
Initialization 

3 
8 
13 
18 
23 
28 








*10 

30 
80 
130 
180 
230 
280 
carried over 

1 
0 
0 
0 
0 
0 
sum 
3 
31 
80 
130 
180 
230 
280 
remainder 

1 
20 
130 
180 
230 
280 








*10 

10 
200 
1300 
1800 
2300 
2800 
carried over 

4 
48 
75 
112 
135 
0 
sum 
1 
14 
248 
1375 
1912 
2435 
2800 
remainder 

4 
8 
31 
262 
251 
352 








*10 

40 
80 
310 
2620 
2510 
3520 
carried over 

41 
12 
120 
112 
180 
0 
sum 
4 
41 
92 
430 
2732 
2690 
3520 
remainder 

1 
32 
94 
92 
506 
256 








*10 

10 
320 
940 
920 
5060 
2560 
carried over 

5 
30 
45 
252 
135 
0 
sum 
1 
15 
350 
985 
1172 
5195 
2560 
remainder 

5 
50 
145 
182 
281 
112 








*10 

50 
500 
1450 
1820 
2810 
1120 
carried over 

9 
54 
75 
140 
45 
0 
sum 
5 
59 
554 
1525 
1960 
2855 
1120 
remainder 

9 
14 
13 
310 
125 
304 








*10 

90 
140 
130 
3100 
1250 
3040 
carried over 

2 
6 
135 
56 
135 
0 
sum 
9 
92 
146 
265 
3156 
1385 
3040 
remainder 

2 
26 
97 
186 
293 
592 








*10 

20 
260 
970 
1860 
2930 
5920 
carried over 

4 
36 
90 
140 
315 
0 
sum 
2 
24 
296 
1060 
2000 
3245 
5920 
remainder 

4 
56 
52 
20 
515 
208 
back to home page 