View on GitHub

Quorten Blog 1

First blog for all Quorten's blog-like writings

Number e… a very important number that most of us have memorized to about 10 digits, i.e. 2.71828182846, but how do you compute it? Let’s review the definition of number e as it is typically presented, in the problem of compound interest.

Non-compound interest is defined as follows:

A = P * r * t

“Annual equals principal times rate times time.” This tells us how much interest money is generated over the given time period. If we want to know the total money, we rephrase the equation as follows, using I instead of r.

P_t = P_0 * (1 + I * t)

Non-compound quarterly interest that is assigned a yearly interest rate is defined as follows.

P_t = P_0 * (1 + I/4 * (4*t))

This can be trivially seen to be the same as the annual interest equation.

Interest compounded annually is defined primarily as a recurrence relation.

P_(t+1) = P_t * (1 + I)

By some simple analysis, you see that you can rewrite as an equation this by raising (1 + I) to the power of t.

P_t = P_0 * (1 + I)^t

Now, what if we want to compound quarterly rather than annually? Simply put, if t is equal to one year, we want to effectively have compounded 4 times rather than just one, so we put a factor of 4 in the exponent. But since we want the overall yearly interest to be about the same, we also divide I by 4.

P_t = P_0 * (1 + I/4)^(4*t)

We can compound at smaller time intervals by defining a variable n for the compounding frequency.

P_t = P_0 * (1 + I/n)^(n*t)

Now, what happens if we simplify this to the case where I = 100%, P_0 = 1, and t = 1? We get the following simplified equation.

P_1 = (1 + 1/n)^n

If we evaluate this with different values of n, we can see that it converges to constant as n approaches infinity. That constant is number e.

lim_{n -> oo} (1 + 1/n)^n = e

However, if we evaluate this directly with a calculator using a conventional number of floating point bits of precision, what kinds of values do we get in practice?

n = 10, e = 2.594
n = 100, e = 2.705
n = 1000000, e = 2.7182805
n = 1000000000, e = 2.71828170574
n = 1000000000000, e = 1.0
n = 100000000000, e = 2.71828129627

These results are rather abysmal since 1/n gets so tiny as n grows large, and then we have to try to raise what looks like a roundoff error to a very high power. Surely, there’s a better way to compute number e, isn’t there? Indeed, there is. Wikipedia gives an unexplained magic solution that is much more numerically stable to evaluate and converges rapidly.

e = 1/0! + 1/1! + 1/2! + 1/3! + 1/4! + ...
= sum_{n = 0}^{oo} 1/n!

20200602/http://en.wikipedia.org/wiki/Number_e

But how and why does that work? Here, I present the explanation that is the missing link between these two equations.

The trick to reformulating the equation is the binomial theorem. The binomial theorem presents an easy way to expand any binomial raised to the nth power into a polynomial series of terms.

P(n, 2) = n*(n - 1)
P(n, 3) = n*(n - 1)*(n - 2)
P(n, r) = n!/(n - r)!
C(n, r) = P(n, r)/r! = n!/((n - r)! * r!)

(a + b)^n = C(n, 0) * a^n + C(n, 1) * a^(n - 1)*b +
  C(n, 2) * a^(n - 2)*b^2 + C(n, 3) * a^(n - 3)*b^3 + ... +
  C(n, n) * a^(n - n)*b^n
= sum_{x = 0}^n C(n, x) * a^(n - x)*b^x

For small binomial coefficients, they can also be easily computed using Pascal’s triangle.

         1
        1 1
       1 2 1
      1 3 3 1
     1 4 6 4 1
   1 5 10 10 5 1
  1 6 15 20 15 6 1
 1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1

Fortunately, the number e equation is a binomial raised to the nth power, so we can use the binomial theorem expansion to rephrase the equation. Let’s do that.

n = 5 expansion:

a^5 + 5*a^4*b + 10*a^3*b^2 + 10*a^2*b^3 + 5*a*b^4 + b^5
= 1^5 + 5*1^4*(1/n) + 10*1^3*(1/n)^2 + 10*1^2*(1/n)^3 + 5*1*(1/n)^4 +
  (1/n)^5
= 1 + 5/n + 10/n^2 + 10/n^3 + 5/n^4 + 1/n^5
= 1 + 5/5 + 10/5^2 + 10/5^3 + 5/5^4 + 1/5^5
= 1 + 1 + 2/5 + 2/5^2 + 1/5^3 + 1/5^5

= 1 + 1 + 0.4 + 0.08

Factorial series:

1/0! + 1/1! + 1/2! + 1/3!
= 1 + 1 + 0.5 + 0.167

n = 8 expansion:

1 + 8/8 + 28/8^2 + 56/8^3
= 1 + 1 + 0.438 + 0.109

What do we see? Looks like we’re seeing pretty much the start of the factorial series in the first few terms. What about the last few terms? Think about it, if n approaches infinity, you will be raising 1/n to a very high power, so the term will be insignificantly small and not really factor into the final result. So, let’s see what we can find out if we take the limit as n approaches infinity on the first few terms. Fortunately, since the first term in the binomial is one, that pretty much drops out of the simplified term equation.

C(n, x) * (1/n)^x

x = 2
n*(n - 1) / (2! * n^2)
lim_{n -> oo} n*(n - 1) / (2! * n^2) = 1/2!

x = 3
n*(n - 1)*(n - 2) / (3! * n^3)
lim_{n -> oo} n*(n - 1)*(n - 2) / (3! * n^3) = 1/3!

Okay, the pattern should be clear now. Computing permutations at the limit? When you try to add or subtract a constant from infinity, it is negligibly small compared to your very large number, so it can basically be dropped out and ignored. Then you can just substitute n anywhere you see n - c. Happy, because of the structure of our binomial, we are also dividing by n raised to an equivalent power as the factors of the permutation, so those factors cancel out, leaving us with just the 1/x! factor from our combinatorial computation. Now that we know each term can be simplified to 1/x! when n approaches infinity, we can rewrite our limit-transformed binomial theorem expansion:

e = sum_{x = 0}^{oo} 1/x!

And that’s how we derive the factorial series that Wikipedia just magically presented without fully explaining. (Except, in hindsight, for giving a vague hint.)

With that knowledge in hand, we can not only happily compute e to whatever number of bits of precision we may need, but we also have an excellent numerically stable computational tool for computing any binomial raised to a power. When a binomial has rational terms, we can rephrase it using the binomial theorem to get a more efficient computational form.


Now, what about that vague hint that Wikipedia gives? Taking a second look at a Calculus book, I see what they are hinting at. They are hinting at using the Maclaurin series expansion of e^x, which results in a power series. Then, the power series is evaluated at one to compute e. However, the problem with this technique is that it requires you to know the derivate of e^x in advance, and how can you even determine the derivative rule for that accurately if you don’t even have an accurate computational result for the number?

Well, maybe it’s worth a try just for posterity. The function e^x can be expressed as follows:

e^x = lim_{n -> oo} (1 + 1/n)^(n*x)

Now we want to use Newton’s difference quotient on this.

lim_{h -> 0} (e^(x + h) - e^x) / h

Now this is looking ugly, but let’s try it anyways.

lim_{h -> 0} lim_{n -> oo} ((1 + 1/n)^(n*(x+h)) - (1 + 1/n)^(n*x)) / h
= lim_{h -> 0} lim_{n -> oo} ((1 + 1/n)^(n*x+n*h) - (1 + 1/n)^(n*x)) / h
= lim_{h -> 0} lim_{n -> oo} ((1 + 1/n)^(n*x) * (1 + 1/n)^(n*h) -
                              (1 + 1/n)^(n*x)) / h
= lim_{h -> 0} lim_{n -> oo} (1 + 1/n)^(n*x) * ((1 + 1/n)^(n*h) - 1) / h

Now if h is approaching zero, one of the exponents will approach zero. On the other hand, the other exponent is approaching infinity. But hold on, we can use a limit rule to simplify this right here.

= lim_{n -> oo} (1 + 1/n)^(n*x) * lim_{h -> 0} ((1 + 1/n)^(n*h) - 1) / h
= e^x * lim_{n -> oo} lim_{h -> 0} ((1 + 1/n)^(n*h) - 1) / h

Now we can clearly see there is a factor of e^x in this mix, independent of the difference quotient limit. It’s hard to do very much more here, other than to show that with some simple rephrasing, we’re almost right back to where we started.

= e^x * lim_{n -> oo} lim_{h -> 0} (1 + 1/n)^(n*h) / h - 1 / h
= e^x * lim_{h -> 0} (e^h - 1) / h

I’m pretty much stuck trying to attempt that.

By the way, my silly Calculus textbook (Single Variable Calculus: Early Transcendentals, Seventh Edition, by James Stewart, published by Brooks/Cole Copyright 2012, 2008, ISBN-10 0-538-49867-6) simply defined number e such that the following holds true:

lim_{h -> 0} (e^h - 1) / h = 1

Well that sure makes it easy to finish this limit computation!

Now, if we computed number e in advance, with sufficient accuracy (thus only possible using the rephrased formula), that is sure easy to believe why it would hold true.

2.71828182846^0.000001 - 1 = 0.000001
0.000001 / 0.000001 = 1.0

And indeed, my derivation using limits on terms of the binomial theorem does not rely on assuming you know the derivative of e^x. But, how you would prove the previous symbolically without prior knowledge of the computational result, I am not sure of such a method at the moment.

Well… now here’s the real trick. By circular reasoning about why we definitely know power series to hold true regardless of not knowing the value of number e, we can start my making that assumption and then prove that the result of evaluating the power series does in fact result in a consistent computation formula for number e. And that’s how we prove that the derivative of e^x is in fact e^x, and therefore the limit also holds true.

Long story short, it’s easier to be able to directly compute number e without needing to rely on any assumptions about its derivative when phrasing up your computational formula.