Digits of e

Larry Bartholdi & his gang

July 24, 2024

We give here a simple algorithm to compute digits of E. Its value is naturally available as a real number (type |(expt 1.)|), but it is significantly more difficult to determine more than, say, the 16 digits of 64-bit reals, for we must now develop the algorithm to compute the exponential.

It is a known fact (almost a definition !) that

ex = $\displaystyle \lim_{{n\rightarrow\infty}}^{}$$\displaystyle \left(\vphantom{1+\frac{x}{n}}\right.$1 + $\displaystyle {\frac{{x}}{{n}}}$$\displaystyle \left.\vphantom{1+\frac{x}{n}}\right)^{n}_{}$.

If we take x = 1, we get a formula which we can expand in a Newton series:

e = $\displaystyle \lim_{{n\rightarrow\infty}}^{}$$\displaystyle \left(\vphantom{1+\frac{n}{n\cdot 1!}+
\frac{n\cdot(n-1)}{n^2\cdot 2!}+\ldots}\right.$1 + $\displaystyle {\frac{{n}}{{n\cdot 1!}}}$ + $\displaystyle {\frac{{n\cdot(n-1)}}{{n^2\cdot 2!}}}$ + …$\displaystyle \left.\vphantom{1+\frac{n}{n\cdot 1!}+
\frac{n\cdot(n-1)}{n^2\cdot 2!}+\ldots}\right)$;

taking the limit yields (n - i)/n→1, thus

e = $\displaystyle \sum_{{i=0}}^{\infty}$$\displaystyle {\frac{{1}}{{i!}}}$.

Of course the numbers involved get extremely large or small, so we need some kind of computational trick to obtain a meaningful result. The idea is to represent e as a fixed-point multi-digit number in a variable base, and then convert it to decimal. Let us first refine these concepts.

A number x in b-base is a tuple (…x1x0.x-1…) such that

x = $\displaystyle \sum_{{i=-\infty}}^{{\infty}}$bixi,

with 0≤xi < b. This can be understood as a special case where b itself is a tuple whose values are all identical. We may then state:

Definition -base   A number x in $\bf b$-base is a tuple (…x1x0.x-1…) such that

x = $\displaystyle \sum_{{i=0}}^{{\infty}}$$\displaystyle \left(\vphantom{\prod_{j=1}^i b_j}\right.$$\displaystyle \prod_{{j=1}}^{i}$bj$\displaystyle \left.\vphantom{\prod_{j=1}^i b_j}\right)$xi + x0 + $\displaystyle \sum_{{i=-\infty}}^{0}$$\displaystyle \left(\vphantom{\prod_{j=i}^{-1}b_j^{-1}}\right.$$\displaystyle \prod_{{j=i}}^{{-1}}$bj-1$\displaystyle \left.\vphantom{\prod_{j=i}^{-1}b_j^{-1}}\right)$xi

with 0≤xi < bi, and of course all these values in N.

If we take the base (bi) with bi = - i, e will be represented as 1.111111… All we have to do is convert this tuple in a decimal representation. But that's enough of chatting! let's get to code.

The numbers are represented as infinite sequences extending one way using streams (quick, dash to your manual to check your knowledge). The digits will then be computed on demand, with no limitation except memory and time.

This routine will convert the e-base number |x| to decimal. The basic algorithm goes:

  1. leave the first digit unchanged.
  2. multiply the rest by 10, and shift it right one position. (this may be seen as a no-op in base 10. For instance, 1.11111111… will become 1.0(10)(10)(10)…). Using this initialization, the digits will always be a little late in the stream, so we have some time to make adjustments.
  3. normalize it; that is, if the next digit is too big, reduce it and increment the current one.
  4. recurse on the normalization.
  5. recurse back to (1).

Note the |(normalize 2 ...)| calls normalize with a source base of 2. Normalize increments this value for successive digits.

(define (convert x) (cons-stream (head x) (convert (normalize 2 (cons-stream 0 (mult (tail x)))))))

Multiply a number by 10. (define (mult x) (cons-stream (* 10 (head x)) (mult (tail x))))

Normalize the number. If its second digit is small enough, leave the current one unchanged; otherwise bump up the current digit and down the second one. Then call normalize recursively on the tail, incrementing the source base (|c|).

(define (normalize c x) (let ((d (head x)) (e (head (tail x)))) (if (< (+ e 9) c) (cons-stream d (normalize (+ c 1) (tail x))) (carry c (cons-stream d (normalize (+ c 1) (tail x)))))))

(define (carry c x) (cons-stream (+ (head x) (quotient (head (tail x)) c)) (cons-stream (remainder (head (tail x)) c) (tail (tail x)))))

This primitive will be useful in extracting a digit from the stream. |(nth 10)| would return the 10th digit of e.

(define (nth digit) ((named-lambda (n d s) (if (= d 0) (head s) (n (-1+ d) (tail s)))) digit e))

Here we cheat a little. E should be represented as 1.111…, but we prefer the representation 0.2111…, although this violates our definition. This is why the source base starts at 2 in |convert|.

(define ones (cons-stream 1 ones)) (define e (convert (cons-stream 2 ones)))

And finally something useful...

(define (digits n) (map nth ((named-lambda (loop from to) (if (> from to) '() (cons from (loop (1+ from) to)))) 0 n)))

This code could be described as giving an infinite number of digits in an infinite number of time. What we would like to know is how much time is needed to compute n digits.

Let us suppose k digits are used in e. Then |convert| gets called k times, calling itself |normalize| k times. The fact that a 0 gets inserted in the stream in |convert| shows us that k $\approx$ n/2. Then the running time is O(n2).