Digits Of E

June 19, 2012

We begin with the algorithm by Rabinowitz and Wagon:

(define-generator (e-spigot n)
  (define (times10 x) (* 10 x))
  (yield 2)
  (let loop1 ((ts (make-list n 10)) (k n))
    (if (< 1 k)
      (let loop2 ((ts ts) (rs (list)) (i (+ n 1)) (carry 0))
        (if (= i 1)
            (begin (yield carry) (loop1 (map times10 (reverse rs)) (- k 1)))
            (let* ((x (+ (car ts) carry)) (q (quotient x i)) (r (remainder x i)))
              (loop2 (cdr ts) (cons r rs) (- i 1) q))))))
  (let loop3 () (yield #f) (loop3)))

The outer loop processes the rows of the tableau one by one, producing another digit each time, and the inner loop processes each row. The third loop returns an indication that processing is complete each time it is called.

The Gibbons algorithm is nearly as pretty in Scheme as in Haskell; like Niemeijer, I made a few changes in translation:

(define-generator (make-spigot f lo hi)
  (define (split v)
    (values (vector-ref v 0) (vector-ref v 1) (vector-ref v 2)))
  (define (approx abc n)
    (let-values (((a b c) (split abc)))
      (quotient (+ (* a n) b) c)))
  (define (mul abc def)
    (let-values (((a b c) (split abc)) ((d e f) (split def)))
      (vector (* a d) (+ (* a e) (* b f)) (* c f))))
  (define (g k) (let-values (((n d a) (f k))) (vector n (* a d) d)))
  (let loop ((z (vector 1 0 1)) (k 1))
    (let ((lbound (approx z lo)))
      (cond ((= lbound (approx z hi))
              (yield lbound)
              (loop (mul (vector 10 (* -10 lbound) 1) z) k))
      (else (loop (mul z (g k)) (+ k 1)))))))

(define pi-spigot (make-spigot (lambda (k) (values k (+ k k 1) 2)) 3 4))
(define e-spigot (make-spigot (lambda (k) (values 1 k 1)) 1 2))

We use a three-slot vector for the tuple, and values to split tuples into pieces; Scheme’s lack of pattern-matching hurts us here.

We used generators from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/my56jfCf.

Pages: 1 2 3

Posted by programmingpraxis
Filed in Exercises
7 Comments »

7 Responses to “Digits Of E”

  1. Programming Praxis – Digits Of E « Bonsai Code said
    June 19, 2012 at 3:39 PM

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

  2. Remco Niemeijer said
    June 19, 2012 at 3:40 PM

    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. Programming Praxis – Billboard Challenge, Part 1 « Bonsai Code said
    June 22, 2012 at 10:19 AM

    […] First we reuse the unbounded spigot algorithm for calculating e from the last exercise; […]

  4. Miracle said
    June 22, 2012 at 5:02 PM

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

  5. Mike said
    June 23, 2012 at 11:40 PM

    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
    July 11, 2012 at 12:14 AM

    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
  7. Lucio said
    February 15, 2018 at 3:27 AM

    Hi Mike, when I tried your code it gave me error can’t iterate on a function, fixed the problem by changing the last two lines to

    print ”.join(str(d) for d in islice(pi_digits(), 10))
    print ”.join(imap(str, islice(e_digits(), 14)))

    in other words replacing the function pi_digits with pi_digits() which invokes the function and returns an array, similarly replacin e_digits with e_digits()

Leave a comment