## Digits Of E

### June 19, 2012

We gave an algorithm for computing the digits of π in a previous exercise. Today, we look at two algorithms for computing the digits of *e*.

We begin with an algorithm due to Stanley Rabinowitz and Stan Wagon:

Algorithm

e-spigot: compute the firstndecimal digits ofe:1.

Initialize:Let the ﬁrst digit be 2 and initialize an arrayAof lengthn+ 1 to (1, 1, 1, . . . , 1).2. Repeat

n− 1 times:

Multiply by 10:Multiply each entry ofAby 10.

Take the fractional part: Starting from the right, reduce theith entry ofAmoduloi+ 1, carrying the quotient one place left.

Output the next digit:The ﬁnal quotient is the next digit ofe.

We give an example of the calculation for the first ten digits of *e* on the next page, where we compute the first ten digits of *e* as 2 7 1 8 2 8 1 8 2 6. As you can see, this algorithm suffers from the fact that the last digit may be wrong (it should be 8, not 6); in fact, as the paper suggests, there are circumstances where several of the last digits may be wrong. The algorithm also suffers from the fact that it is bounded, meaning that the number of digits must be specified in advance, and it needs space proportional to *n*^{2}.

A different algorithm comes from Jeremy Gibbons, and is both unbounded and requires only constant space; Gibbons gave the sequence for π, but Tom Moertel adapts it to *e*. When I was unable to work out the algorithm from Moertel’s description, I asked Remco Niemeijer, a regular contributor to Programming Praxis, for help, and he responded with this gorgeous hunk of Haskell code, which computes both π and *e*:

`stream :: Integral a => (a, a) -> (a, a, a) -> [(a, a, a)] -> [a]`

stream (lo, hi) z ~(x:xs) = if lbound == approx z hi

then lbound : stream (lo, hi) (mul (10, -10*lbound, 1) z) (x:xs)

else stream (lo, hi) (mul z x) xs

where lbound = approx z lo

approx (a,b,c) n = div (a*n + b) c

mul (a,b,c) (d,e,f) = (a*d, a*e + b*f, c*f)

```
```streamDigits :: (Num a, Integral a, Enum a) => (a -> (a, a, a)) -> (a, a) -> [a]

streamDigits f range = stream range (1,0,1) [(n, a*d, d) | (n,d,a) <- map f [1..]]

stream_pi, stream_e :: [Integer]

stream_pi = streamDigits (\k -> (k, 2*k + 1, 2)) (3, 4)

stream_e = streamDigits (\k -> (1, k , 1)) (1, 2)

`main :: IO ()`

main = do print $ take 30 stream_pi

print $ take 30 stream_e

Niemeijer explained that he had adapted the algorithm in Gibbons’ paper as follows:

–

`unit`

is a constant used in only one place. Remove the definition and just use the value directly.–

`comp`

is basic 2×2 matrix multiplication. Rename to`mul`

to make purpose clearer.–

`extr`

as defined in the paper didn’t compile for me: it should be`fromInteger (q * x)`

instead of`fromInteger q * x`

. But there’s more to be gained: the only place`extr`

is used is preceded by`floor`

, so essentially we have`floor (a / b)`

, which is equal to`div a b`

. Since it’s used to approximate the real value I named it`approx`

.– The code used 2×2 matrices expressed as a 4-tuple. However, the bottom-left element is always 0. Remove this element everywhere, leaving 3-tuples and making

`mul`

and`approx`

simpler.– For both π and

e, the`next`

and`safe`

functions passed to`stream`

are identical save for the bounds used, so pass in the bounds and integrate the functions in`stream`

.– The

yused in`stream`

is the lower bound. Rename accordingly.– The

`cons`

passed to stream is always`comp`

(or in my case`mul`

). Remove the argument and use`mul`

directly.–

`prod`

is also the same for π ande, so remove the argument and inline the function.– The definition of π from the paper calls

`stream`

with some starting arguments. Sinceewill be doing the same, make a function`streamDigits`

to abstract this out.– Since

`stream`

now takes far fewer parameters, there’s no longer a reason to name them all.– According to the hulver site, each term in the π function is 2 +

k/(2k+ 1)x. Fore, each term is 1 + (1/k)x(or rather, it definese− 1 starting fromk= 2, but if you start fromk= 1 you add the term 1 + 1/1x, which gives you 1 + 1/1*(e-1), which of course ise). So both are of the forma+ (n/d)x. The matrix used in`lfts`

is (n,a*d, 0,d). I didn’t like having to do the multiplication myself, so I made`streamDigits`

take a function that producesn,dandafor a given term to stay closer to the mathemetical definition. This function is used to produce the actual matrix.– Add

`stream_e`

and`stream_pi`

functions.

Your task is to write the two spigot functions for *e* described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.