Approximating e with the magical Log Gamma function

Previously, I wrote a high-level overview on the log-sum-exp function, and how it can be used to reduce overflow and underflow errors. In this post, we’ll demonstrate how you can leverage the Log Gamma function to accurately approximate e without overflow error.

Computationally Approximating e

There are many representations of e. The infinite series representation is very elegant:

infinite series representation of e

So let’s try to approximate e as close as possible! That is, let’s keep adding terms in the infinite series until our computers explode.

Here’s what a naive implementation in Python might look like.

from math import factorial# input n is the number of terms in 
# your series. higher => more accurate
factorials = [factorial(i) for i in range(n)]
one_overs = [1.0/i for i in factorials]
return sum(one_overs)

And this isn’t that bad. It gets us pretty dang close to e (2.7182818284590452353…)

brute_force(150)Out[142]: 2.7182818284590455

But suppose that we want even more accuracy. Well, then let’s just increase the number of terms we sum from 150 to 1000.

brute_force(1000)OverflowError                             Traceback (most recent call last)
<ipython-input-146-8ec3a5075e84> in <module>()
1 # Brute force has overflow error at n=1000
----> 2 brute_force(1000)

Uh-oh. We got an OverflowError. What’s intuitively going on here is that, as you can imagine, the factorial terms are blowing up really fast and causing our computer to complain. Intuitively, we can fix our overflow problem by avoiding direct calculations of the factorial terms.

Log Gamma Function to the Rescue

Thankfully, we have mathematics. Without spoiling the good stuff, we will convert our elegant representation of e into a mathematically equivalent representation. What our new representation will lack in elegance, it will make up for in computability!

First, we shall derive an alternate expression for each element in the series.

The first two steps are just using simple properties of exponents and logs. The last step is just converting the factorial function into its equivalent log gamma function, which is just an extension of the gamma function, which is itself an extension of the factorial function. ¯\_(ツ)_/¯

Now let us plug in our new expression into the infinite series.

Once we plug in our new expression, we essentially have a sum of exponentiations of the log-gamma function. Now let’s code it up in Python.

from math import exp
from scipy.special import gammaln
def smart_way(n):
#note: we shift everything up one bc we are using logGamma
neg_log_gammas = [-gammaln(float(i)) for i in range(1,n+1)]
exponentiated = [exp(i) for i in neg_log_gammas]
return sum(exponentiated)
# we're still going at n = 1 million

Notice how now we are able to approximate e to 1 million terms, whereas the naive implementation couldn’t even calculate 1000. Note that I chose “1 million” arbitrarily — it can still probably go even higher!

So why does this work? Well, by using the log-gamma function, we have successfully avoided directly calculating the factorial terms, so we avoided that source of overflow error.

how log-gamma is calculated. no factorial terms!

In summary, math is fun and useful.

Data Scientist at Instagram.