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 Log2(10n)=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 Log13(10n)=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 :