Calculate e in Python [duplicate]

Question 1

The obvious solution would be

import math

def e(n=10):
    return sum(1 / float(math.factorial(i)) for i in range(n))

but it loses precision around n=20 (the error as compared to math.e is around 10^-16)

40-digits precision might be a challenge, and might require arbitrary precision arithmetic

I do not really see the point to have such precise "e" value since you won't be able to perform any calculations with such precision (if you do anything to it, you will lose that precision, unless you do everything in some arbitrary precision arithmetic).

Question 2

If you need more precision, you could try to use Fraction:

from fractions import Fraction # use rational numbers, they are more precise than floats
e = Fraction(0)
f = Fraction(1)
n = Fraction(1)
while True:
    d = Fraction(1) / f # this ...
    if d < Fraction(1, 10**40): # don't continue if the advancement is too small
        break
    e += d # ... and this are the formula you wrote for "e"
    f *= n # calculate factorial incrementally, faster than calling "factorial()" all the time
    n += Fraction(1) # we will use this for calculating the next factorial
print(float(e))

or Decimal:

from decimal import Decimal, getcontext
getcontext().prec = 40 # set the precision to 40 places
e = Decimal(0)
f = Decimal(1)
n = Decimal(1)
while True:
    olde = e
    e += Decimal(1) / f
    if e == olde: # if there was no change in the 40 places, stop.
        break
    f *= n
    n += Decimal(1)
print(float(e))

So here is e in 1000 places:

2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035044

To see more clearly what it does, here is its simplified version:

e = f = 1.0
for i in range(2, 16):
    e += 1.0 / f
    f *= i
print(e)

Question 3

To produce a similar result to my answer from a question on approximating pi:

from functools import wraps

def memoize(f):
    """Store and retrive previous results of the decorated function f."""
    cache = {}
    @wraps(f)
    def func(*args):
        if args not in cache:
            cache[args] = f(*args)
        return cache[args]
    return func

@memoize
def fact(n):
    """Recursively calculates n!."""
    if n <= 1:
        return 1
    return n * fact(n - 1)

def inverse_fact_n(start_n=0):
    """Generator producing the infinite series 1/n!."""
    numerator = 1.0
    denominator = start_n
    while True:
        yield numerator / fact(denominator)
        denominator += 1

def approximate_e(steps=None, tolerance=None):
    """Calculate an approximation of e from summation of 1/n!."""
    if steps is None and tolerance is None:
        raise ValueError("Must supply one of steps or tolerance.")
    series = inverse_fact_n()
    if steps is not None: # stepwise method
        return sum(next(series) for _ in range(steps))
    output = 0 # tolerance method
    term = next(series)
    while abs(term) > tolerance: 
        output += term
        term = next(series)
    return output

if __name__ == "__main__":
    from math import e
    print("math.e:\t\t{0:.20f}.".format(e))
    stepwise = approximate_e(steps=100)
    print("Stepwise:\t{0:.20f}.".format(stepwise))
    tolerated = approximate_e(tolerance=0.0000000001)
    print("Tolerated:\t{0:.20f}.".format(tolerated))

The function approximate_e allows you to specify either:

A number of steps ("I want it to take this long"); or
A desired tolerance ("I want it to be this accurate").

There is some relatively advanced Python around this (e.g. the memoizing decorator function and the generator function to produce the series), but you can just focus on the main function, where next(series) gives you the next term of the summation.

This gives me the output:

math.e:     2.71828182845904509080.
Stepwise:   2.71828182845904553488.
Tolerated:  2.71828182844675936281.

Question 4

The most effective way is to use the properties of the exponential function.

exp(x)=(exp(x/N))^N

Thus you compute x=exp(2^(-n)) with 2n bits more precision than required in the final result, and compute e by squaring the result n times.

For small numbers x, the error of truncating the series for exp(x) at the term with power m-1 is smaller than two times the next term with power m.

To summarize, to compute e with a precision/accuracy of d bit, you select some medium large n and select m such that

2^(1-mn)/m! is smaller than 2^(-d-2n)

This determination of m can also be done dynamically (using Decimal as in the answer of user22698)

from decimal import Decimal, getcontext


def eulernumber(d):
    dd=d
    n=4
    while dd > 1:
        dd /= 8
        n += 1
    getcontext().prec = d+n
    x = Decimal(1)/Decimal(1 << n)
    eps = Decimal(1)/Decimal(1 << (1 + (10*d)/3 ))
    term = x
    expsum = Decimal(1) + x
    m = 2
    while term > eps:
        term *= x / Decimal(m)
        m += 1
        expsum += term
    for k in range(n):
        expsum *= expsum
    getcontext().prec = d
    expsum += Decimal(0)
    return expsum


if __name__ == "__main__":
    for k in range(1,6):
        print(k,eulernumber(4*k))

    for k in range(10,13):
        print(k,eulernumber(4*k))

with output

( 1, Decimal('2.718'))
( 2, Decimal('2.7182818'))
( 3, Decimal('2.71828182846'))
( 4, Decimal('2.718281828459045'))
( 5, Decimal('2.7182818284590452354'))
(10, Decimal('2.718281828459045235360287471352662497757'))
(11, Decimal('2.7182818284590452353602874713526624977572471'))
(12, Decimal('2.71828182845904523536028747135266249775724709370'))

See the (unix/posix) bc math library for a more professional implementation of this idea, also for the logarithm and trig functions. The code of the exponential function is even given as example in the man page.