Approximating mathematical constant e in Python

Question

I'm attempting to approximate the constant e using 100,000 iterations and display the results where i is in multiples of 10,000 ranging from 10,000 to 100,000.

The code I have is this:

import math

e = 0
i = 0

while i < 100001:
    e += 1 / (math.factorial(i))
    i += 1
    if i % 10000 == 0:
        print("The approximation using", i, "iterations is: ", e)

There's two problems with this:

1.) The code takes forever to run, but I think that might be expected.

2.) I'm getting the same value for the approximation no matter how many iterations I use. Here's the beginning portion of my output:

The approximation using 10000 iterations is 2.7182818284590455
The approximation using 20000 iterations is 2.7182818284590455
The approximation using 30000 iterations is 2.7182818284590455

.. and so on. How can I fix this?

For context, this is an intro-level programming problem that only incorporates loops.

Edit: Noticed that I had my i and e sections backwards in my while loop. Fixed. Still encountering the issue of the number not changing, though.

Answer 1

Float numbers are encoded in 64 bits with means you can't keep iterating and getting new digits of e.

one alternative would be to use arbitrary precision arithmetic (GMP)

Install gmpy2 :

pip3 install gmpy2

Refractor your code :

import math
import gmpy2

e = gmpy2.mpq(1)
i = 0

while i < 100001:
    i += 1
    e += gmpy2.mpq(1) / gmpy2.mpq(math.factorial(i))
    print("The approximation using", i, "iterations is: ", gmpy2.mpfr(e,precision=i))

Answer 2

First thing, you really should start accumulating your value for e when i==0 - the snippet above starts from i==1 and that skips a full unit - that is why you get 1.718... instead of 2.718... - and the second thing is that, despite Python integer being able to have an arbitrary size (and thus, one can calculate fact(100_000) in a fraction of a second), for decimal numbers, Python will use 64bit floating point numbers by default - these can only hold that many decimal places, and you won't get any better precision past fact(25) (if that much).

You can resort to use Python's decimal.Decimal instead of ordinary floats - those have a configurable precision, and you could get to a few hundred decimal places using this method, at least, until everything gets too CPU or memory consuming.

from decimal import Decimal, getcontext

getcontext().prec = 200 # use 200 digits of precision:

def calc_e(iterations):
    to_e = 0
    for i in range(iterations + 1):
        to_e += Decimal(1) / math.factorial(i)

        if i % 20 == 0:
            print("The approximation using", i + 1, "iterations is: ", to_e)

Please note that the function displaying a lot of digits does not mean those displayed digits match the value for the correct digit in the e number already you have to either use a reference e value (other than the one shipped in Python's math.e which is also a float 64 number), or another formula, to check up to which digit your value is correct.

Answer 3

First of all, you should start with e = 1 to get correct values of e, your code right now results in value of e - 1.

You're working on wrong orders of magnitude - your code (after changing initial e value) converges to the most precise value floating point precision allows on 17th iteration, nowhere near the 10000 step you attempted:

The approximation using 15 iterations is:  2.71828182845823
The approximation using 16 iterations is:  2.718281828458995
The approximation using 17 iterations is:  2.718281828459043
The approximation using 18 iterations is:  2.7182818284590455
The approximation using 19 iterations is:  2.7182818284590455
The approximation using 20 iterations is:  2.7182818284590455
...

Any further iteration will end up on the same number, as the value of 1 / (math.factorial(i)) is so small that adding it to current value of e (2.7182818284590455) results in the same floating point number.

If you want arbitrary precision, you can use Decimal module:

import math
from decimal import *

getcontext().prec = 100

e = Decimal(1)
i = 0

while i < 100:
    i += 1
    e += Decimal(1) / Decimal(math.factorial(i))
    print("The approximation using", i, "iterations is: ", e)

Even then, with 100 digit precision you will achieve convergence around 70th iteration.

Answer 4

You could also just use scaled ints, e.g., change 1 / to 10**80 //. Then you end up with

271828182845904523536028747135266249775724709369995957496696762772407663035354729

which is accurate this far:

27182818284590452353602874713526624977572470936999595749669676277240766303535