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 first n decimal digits of e:

1. Initialize: Let the first digit be 2 and initialize an array A of length n + 1 to (1, 1, 1, . . . , 1).

2. Repeat n − 1 times:

Multiply by 10: Multiply each entry of A by 10.

Take the fractional part: Starting from the right, reduce the ith entry of A modulo i + 1, carrying the quotient one place left.

Output the next digit: The final quotient is the next digit of e.

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 n2.

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 y used 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 π and e, so remove the argument and inline the function.

- The definition of π from the paper calls stream with some starting arguments. Since e will 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. For e, each term is 1 + (1/k)x (or rather, it defines e − 1 starting from k = 2, but if you start from k = 1 you add the term 1 + 1/1x, which gives you 1 + 1/1*(e-1), which of course is e). So both are of the form a + (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 produces n, d and a for 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.

About these ads

Pages: 1 2 3

6 Responses to “Digits Of E”

  1. […] today’s Programming Praxis exercise, our goal is to implement two algorithms to calculate the digits of e […]

  2. Here’s my Haskell solution for the first algorithm (since my solution for the second one has already been posted in the exercise). A version with comments can be found at http://bonsaicode.wordpress.com/2012/06/19/programming-praxis-digits-of-e/ .

    import Data.List
    
    spigot_e :: Int -> [Int]
    spigot_e n = 2 : take (n - 1) (f $ replicate (n + 1) 1) where
        f = (\(d,xs) -> d : f xs) .
            mapAccumR (\a (i,x) -> divMod (10*x+a) i) 0 . zip [2..]
    
  3. […] First we reuse the unbounded spigot algorithm for calculating e from the last exercise; […]

  4. Miracle said

    Does anybody share a Java or C# code for this exercises?

  5. Mike said

    Basically a direct translation of the haskell code into Python 2.7. I create the input stream and initialize the state vector (z) in ‘stream()’ and eliminated ‘streamDigits()’.

    
    from itertools import count, imap
    
    def stream(lo, hi, f):
        def approx((a,b,c), n): return (a*n + b)//c
        def mul((a,b,c),(d,e,f)): return a*d, a*e + b*f, c*f
        
        xs = ((n, a*d, d) for n,d,a in imap(f, count(1)))
        z = 1, 0, 1
        while True:
            lbound = approx(z, lo)
            if lbound == approx(z, hi):
                yield lbound
                z = mul((10, -10*lbound, 1), z)
    
            else:
                z = mul(z, next(xs))
    
    			
    def pi_digits():
        return stream(3, 4, lambda k: (k, 2*k + 1, 2))
    
    def e_digits():
        return stream(1, 2, lambda k: (1, k, 1))
    
    # test
    from itertools import islice
    
    print ''.join(str(d) for d in islice(pi_digits, 10))   # returns "3141592653"
    print ''.join(imap(str, islice(e_digits, 14)))         # returns "27182818284590"
    
    
  6. David said

    Here is FORTH code for the first algorithm (by Stan & Stanley) Though space is proportional to n, not sure why you mention n**2.

    100 constant #digits
    : int-array  create cells allot  does> swap cells + ;
    
    #digits 1+ int-array e-digits[]
    
    : init-e ( -- )
       [ #digits 1+ ] literal 0 DO
          1  i e-digits[]  !
       LOOP
       cr ." 2." ;
    
    : .e  ( -- )
       init-e
       [ #digits 1- ] literal 0 DO
          0  \ carry
          0 #digits DO
             i e-digits[] dup @  10 *  rot +  i 2 + /mod -rot  swap !
          -1 +LOOP
          0 .r
       LOOP ;
    
    .e
    2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427 ok
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 634 other followers

%d bloggers like this: