Sections

Doing it on the computer
More decimals, please
Why large integers are good for small decimals

This post is part of a series, such as:

One of the nicest things that I learned in my freshman year was how to calculate the Euler’s number (that is, \(e\)) with as many digits as desired. Well, they didn’t teach me that, just enough math to do it on my own. We can calculate it using this formula

\[e= \frac{1}{1} + \frac{1}{1} + \frac{1}{1⋅ 2} + \frac{1}{1\cdot 2⋅ 3} + \frac{1}{1⋅ 2⋅ 3⋅ 4} + ⋯\] and this is easy to calculate. Moreover, the series converges very fast.

Sometimes is easier to understand this formula using the notation \(n! = 1⋅ 2⋅ 3⋯ n\). This is called \(n\) factorial, and it is just the product of all the integers between 1 and \(n\) When \(n=0\) we define \(0!=1\) and everything is coherent. In that case the formula for \(e\) can be beautifully written as \[e= \frac{1}{0!} + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + ⋯ + \frac{1}{k!}+⋯\]

To calculate \(e\) we will proceed step by step. In the step \(k\) we calculate \(e_k,\) the expression with the first \(k\) sums, that is \[e_k= \frac{1}{0!} + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + ⋯ + \frac{1}{k!}.\] It is easy to see that \(e_k\) will converge to \(e\) when \(k\) grows. It is also easy to see that \[e_k=e_{k-1} + \frac{1}{k!}\] so we can calculate a new value as a simple modification of the old value. Now we need to calculate that simple modification. If we give the name \(f_k\) to one over the factorial of \(k\), we have \[f_k = \frac{1}{k!} = \frac{1}{k}f_{k-1}.\] Therefore we have this simple procedure: we start with

\[\begin{aligned} e_0 & =1\\ f_0 & =1, \end{aligned}\] and then we calculate, for every \(k\), the values

\[\begin{aligned} f_k & =f_{k-1}/k\\ e_k & = e_{k-1} + f_k. \end{aligned}\] We keep iterating until \(f_k\) is small. That is, we stop when \(f_k<ε\) for some precision \(ε\) chosen a priori We can see that we need only two operations: dividing a real number by an integer and summing two real numbers.

Doing it on the computer

We can directly implement this algorithm in a program, for instance in Python. Since we only need the last version of \(e_k\) and \(f_k\), we will just call them e and f. They will be updated on every step, and we do not need to store the old values. We just need to print them.

Python is a nice language for this problem, as we will see below. There is one nice detail to be aware of if you are not familiar with the language. The expression k += 1 is a short version of k = k+1. It is also faster since the computer knows that it will update a fixed part of the memory and not create a new one. The same idea is valid for f /= k, representing f = f/k.

The code is this:

e = 1
f = 1
k = 1
while f>1e-15:
    f /= k
    e += f
    k += 1
    print(e)

2.0
2.5
2.6666666666666665
2.708333333333333
2.7166666666666663
2.7180555555555554
2.7182539682539684
2.71827876984127
2.7182815255731922
2.7182818011463845
2.718281826198493
2.7182818282861687
2.7182818284467594
2.71828182845823
2.718281828458995
2.718281828459043
2.7182818284590455
2.7182818284590455

We can get as many decimals as the precision of the floating-point representation of numbers. How can we go further?

More decimals, please

When I first learned this idea, I used a character array to represent e and f, and I wrote my own add and divide functions in the Turbo Pascal language. I remember that it took me 40 minutes to write the program and 20 minutes to run it and get ten thousand digits.

Today it is much easier to do it in Python, and it takes only a few seconds. This language can handle large integers, as large as we want, as long as we do not mix them with real numbers. For example, we can calculate \(2^{521} - 1\) and get a large prime number

2**521 - 1

6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151

Notice that the symbol ** means “to the power of”. If we need to divide an integer, we use integer division, written as //. It does not take any decimals. For example, compare

7 // 2

and

7 / 2

3.5

Why large integers are good for small decimals

Ok, so we can have large integers. How do they help us with our problem of small decimals?

The idea is to multiply \(e_k\) and \(f_k\) by a big constant \(Ω\) (let’s say \(10^{10000}\)), and represent them by large integers. We will get numbers with ten thousand digits that we can add, and we can divide by a small integer. So we will have \(e_k\) converging to \(e⋅ Ω\).

Let’s test this idea with sixty digits, and we can see each step

num_decimals = 60
e = 10**num_decimals
f = 10**num_decimals
k = 1
while f>0:
    f //= k
    e += f
    k += 1
    print(e, k)

2000000000000000000000000000000000000000000000000000000000000 2
2500000000000000000000000000000000000000000000000000000000000 3
2666666666666666666666666666666666666666666666666666666666666 4
2708333333333333333333333333333333333333333333333333333333332 5
2716666666666666666666666666666666666666666666666666666666665 6
2718055555555555555555555555555555555555555555555555555555553 7
2718253968253968253968253968253968253968253968253968253968251 8
2718278769841269841269841269841269841269841269841269841269838 9
2718281525573192239858906525573192239858906525573192239858903 10
2718281801146384479717813051146384479717813051146384479717809 11
2718281826198492865159531826198492865159531826198492865159527 12
2718281828286168563946341724119501897279675057452835230613003 13
2718281828446759002314557870113425668981224536780092335647885 14
2718281828458229747912287594827277366959906642446324986007519 15
2718281828458994464285469576474867480158485449490740496031494 16
2718281828459042259058793450327841862233396624931016465407992 17
2718281828459045070516047795848605061178979635251032698900727 18
2718281828459045226708117481710869683342623135824366934094767 19
2718281828459045234928752728335199400298604372696647683315505 20
2718281828459045235339784490666415886146403434540261720776541 21
2718281828459045235359357431729807147377251008913767151131828 22
2718281828459045235360247110869052204705925898658017397966159 23
2718281828459045235360285792570758511546303067777332626089390 24
2718281828459045235360287404308329607664652116490637427261191 25
2718281828459045235360287468777832451509386078439169619308063 26
2718281828459045235360287471257428714734183538514113165156019 27
2718281828459045235360287471349265613372138999998370333520758 28
2718281828459045235360287471352545502609208837908522375248070 29
2718281828459045235360287471352658602238073315077837962893839 30
2718281828459045235360287471352662372225702130983481815815364 31
2718281828459045235360287471352662493838206286335276778812832 32
2718281828459045235360287471352662497638597041190020371406502 33
2718281828459045235360287471352662497753760397397739874212370 34
2718281828459045235360287471352662497757147554933261036059601 35
2718281828459045235360287471352662497757244330862847354969521 36
2718281828459045235360287471352662497757247019083113641605907 37
2718281828459045235360287471352662497757247091737715433136620 38
2718281828459045235360287471352662497757247093649678638176901 39
2718281828459045235360287471352662497757247093698703335742036 40
2718281828459045235360287471352662497757247093699928953181164 41
2718281828459045235360287471352662497757247093699958846289435 42
2718281828459045235360287471352662497757247093699959558030108 43
2718281828459045235360287471352662497757247093699959574582216 44
2718281828459045235360287471352662497757247093699959574958400 45
2718281828459045235360287471352662497757247093699959574966759 46
2718281828459045235360287471352662497757247093699959574966940 47
2718281828459045235360287471352662497757247093699959574966943 48
2718281828459045235360287471352662497757247093699959574966943 49

We can see that we need only a few iterations since f gets small very quickly. Let’s do it now with ten thousand digits. We will print only the final result.

num_decimals = 10000
e = 10**num_decimals
f = 10**num_decimals
k = 1
while f>0:
    f //= k
    e += f
    k += 1
print("we did", k, "steps")

we did 3250 steps

If everything is right, we should have the result in e. It will be boring to print it here in this webpage. Instead, we write it into a file that you can download if you like.

with open("e-10000-digits.txt", mode="w") as out:
    print(e, file=out)

Now the result should be here. What other famous number can we calculate now?