It has been almost a year since my last post here. I haven’t been idle code-wise, but I have been too distracted by Real Life to do any blogging here. My publication rate is way down on my main blog, too.
To get me posting here more, I think I need more of a focus on casually documenting my own little projects rather than Code Wise articles. Those take more work than I seem willing to put in these days. With that in mind, I have some trivial Python fun to share…
This little project began when I needed a way to illustrate where the special number e (2.71828…) comes from. The special number pi (3.14159…) is easy; it’s just the ratio between a circle’s diameter and circumference. But e is a bit more mysterious. [See this page on my main blog for details.]
One way to calculate the number is to implement the exp function, which is defined:
For as many terms as required by the desired precision. The code is simple and allows for calculating as many terms as desired:
002| from math import e, factorial
003|
004| def exp_factory (nbr_of_terms):
005| def exp_function (x):
006| terms = []
007| for n in range(nbr_of_terms):
008| term = Fraction(pow(x,n), factorial(n))
009| terms.append(term)
010| return float(sum(terms))
011| return exp_function
012|
013| fmt = ‘%3d: %.16f %.16e’
014|
015| for nterms in [1,2,3,5,10,15,20,25]:
016| # Get exp function using NT terms…
017| exp_func = exp_factory(nterms)
018| # Calculate exp(1)…
019| approx_e = exp_func(1)
020| print(fmt % (nterms, approx_e, e–approx_e))
021| print(‘ e: %.16f’ % e)
022| print()
Call exp_factory and pass it the number of terms desired. It returns a function that implements exp with that many terms. Note that, in Python, this could easily be a lambda returning a lambda, making it a one-liner, but shown here as a multi-line function for clarity.
Lines 15–22 illustrate using the factory function. It generates exp(1) for functions with 1, 2, 3, 5, 10, 15, 20, and 25 terms. When run it produces the following output:
1: 1.0000000000000000 1.7182818284590451e+00 2: 2.0000000000000000 7.1828182845904509e-01 3: 2.5000000000000000 2.1828182845904509e-01 5: 2.7083333333333335 9.9484951257116094e-03 10: 2.7182815255731922 3.0288585284310443e-07 15: 2.7182818284582297 8.1534778928471496e-13 20: 2.7182818284590451 0.0000000000000000e+00 25: 2.7182818284590451 0.0000000000000000e+00 e: 2.7182818284590451
The first column is the number of terms used; the second is the result obtained with that many terms; the third is the difference between the calculated value and the library-supplied value (shown on the bottom row). As the data shows, the precision of the result is indistinguishable from the library value with just 20 terms.
§
Which is all well and good, but I wanted to implement the formula that e comes from, the formula for exponential growth:
![\displaystyle{e}=\lim_{{n}\rightarrow\infty}\left[\left({1}+\left(\frac{1}{n}\right)\!\!\right)^{n}\right]](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7Be%7D%3D%5Clim_%7B%7Bn%7D%5Crightarrow%5Cinfty%7D%5Cleft%5B%5Cleft%28%7B1%7D%2B%5Cleft%28%5Cfrac%7B1%7D%7Bn%7D%5Cright%29%5C%21%5C%21%5Cright%29%5E%7Bn%7D%5Cright%5D+&bg=ffffff&fg=000080&s=2&c=20201002)
[See this page for the details.] In particular, I was interested in the terms themselves. I also wanted to create the terms for the complex version. (See them at the bottom of the linked page.)
The first pass at implementing the equation gave me the intermediate steps leading to the final result. Which wasn’t quite what I was going for, but fine for the first pass. Here’s the code:
002| tot = 1.0
003| data = [tot]
004| # Calculate…
005| for ix in range(n):
006| tot *= (1.0 + (1.0/n))
007| if capture_steps: data.append(tot)
008| return data if capture_steps else tot
009|
010| def calc_e_imag (n, capture_steps=0):
011| tot = complex(1,0)
012| data = [tot]
013| # Calculate…
014| for ix in range(n):
015| tot *= complex(1.0, 1.0/n)
016| if capture_steps: data.append(tot)
017| return data if capture_steps else tot
018|
019| terms = 20
020|
021| print(calc_e_real(terms, capture_steps=0))
022| print(calc_e_real(terms, capture_steps=1))
023| print()
024| print(calc_e_imag(terms, capture_steps=0))
025| print(calc_e_imag(terms, capture_steps=1))
026| print()
The first function calculates the real value of e, the second calculates its complex value. The code at the bottom, lines 21–26, demonstrate the two functions. The capture_steps flag determines whether these functions return the list of approximations or just the final result.
With terms set to only 20 (as shown in the code), the two approximations of e are not particularly close. The actual values (per the Python library) are:
- exp(1) = 2.7182818284590451
- exp(i) = (0.5403023058681398+0.8414709848078965i)
And the code above, set to 20 terms, generates:
- calc_e_real(1) = 2.6532977051444226
- calc_e_imag(i) = (0.5546805276912777+0.8622847647277038j)
Close, but no cigar. In fact, it takes a lot of iterations to come close using these (a million doesn’t even come that close). But deriving the value wasn’t the point. The point was generating the intermediate terms, and that the code does that just fine.
§
Which is all well and good, but I still needed the fractions for my page, not just the intermediate (floating-point) values. So, back to the drawing board. Complicating the issue is that deriving the terms requires algebraically expanding the series for as many terms as desired, and only plugging in the value of n once you’ve generated the series.
For example, if you wanted an expression with five terms, then set n=4:
And then do the math. But to get those terms, you have to do the algebra, and implementing all of algebra isn’t a trivial task. Even a decent sub-set would take some doing. But there is a pattern we can leverage that makes it easy. It starts by noticing the numerators of the fractions.
Which I’ll come back to, but we can also take a stab at implementing something algebra-like to generate the desired series. And we also want it to work with complex values. That opens a can of worms I’ll get to in a bit.
The thing to consider is what’s fundamental in the algebra involved. Consider the last iteration of (1+1/4)^4. At that point we have:
To carry out that final multiplication first multiply the left side by one (the first value on the right side) and then multiply the left side by 1/4 (the second value on the right side). This gives us:
Adding like terms gives us the five-term expression above. Effectively, the second set of terms is a copy of the first set shifted to the right one term (because the denominators change due to multiplying by 1/n — in this case, 1/4).
We can leverage this by assuming an initial seed, (1+1/n), and using a shift technique on the numerators. (The denominators always follow the same pattern: n⁰, n¹, n², n³…) It ends up looking like this:
| × 1 | 1 | 3 | 3 | 1 | |
|---|---|---|---|---|---|
| × ¼ | 1 | 3 | 3 | 1 | |
| numerators | 1 | 4 | 6 | 4 | 1 |
| denominators | 40 | 41 | 42 | 43 | 44 |
This time the code is a little more involved:
002|
003| def get_coefficients (depth=9):
004| def next_row (row):
005| left = row + [0]
006| right = [0] + row
007| return [a+b for a,b in zip(left,right)]
008|
009| table = [[1]]
010| for ix in range(depth):
011| row = next_row(table[–1])
012| table.append(row)
013| return table
014|
015| tab = get_coefficients(4)
016| row = tab[–1]
017| print(row)
018| print()
019| base = len(row) – 1
020|
021| terms = [(n,base**ix) for ix,n in enumerate(row)]
022| for t in terms:
023| print(repr(t))
024| print()
025|
026| fracs = [Fraction(a,b) for a,b in terms]
027| for f in fracs:
028| print(repr(f))
029| print()
030|
031| value = float(sum(fracs))
032| print(value)
033| print()
And this is only for the real calculation. The complex version, as implied above, is even more involved (because we have to deal with complex math).
The meat is in the get_coefficients function. It builds a list of rows of coefficients, starting with the one-term zeroth degree (which is just one). The next row has two terms, (1, 1), then comes (1, 2, 1), followed by (1, 3, 3, 1), and finally the five-term (1, 4, 6, 4, 1) row of interest. The list would continue if we specified a greater value for depth.
The for loop (lines 10–12) builds the list with values it gets from the next_row function, which takes the current row, adds it to a right-shifted copy, and returns a new row with one more element.
Lines 16–18 grab the last (and therefore longest) row and print it. Line 21 converts the list of coefficients to a list of tuples with numerators and denominators. (Lines 22–24 print the new list.) Line 26 converts the tuples to Python Fraction objects, which allows us to do math on them. (Lines 27–29 print the new list.)
Lastly, line 31 sums the Fraction terms into a floating-point value, which line 32 prints. The output looks like this:
[1, 4, 6, 4, 1] (1, 1) (4, 4) (6, 16) (4, 64) (1, 256) Fraction(1, 1) Fraction(1, 1) Fraction(3, 8) Fraction(1, 16) Fraction(1, 256) 2.44140625
The result value, 2.4414…, is pretty far from the expected value (2.71828…). But that’s to be expected with only five terms. And note how the Fractions collapse 4/4 into 1/1 and 4/64 into 1/16. Same value, different spelling.
§
Before we get to the complex version, let’s take a moment to look at the pattern of coefficients as the degree goes from zero to ten. They look like this:
[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] [1, 9, 36, 84, 126, 126, 84, 36, 9, 1] [1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1]
You might recognize this pattern as Pascal’s Triangle. It turns out that the coefficients for (1+1/n)^ⁿ follow the same pattern. So, for example, we can know that (1+1/10)¹⁰ gives us the same coefficients as found in the bottom row above.
And by the way, this similarity means we can build the Triangle using the same shift technique we just used above. In fact, all we’re really doing in the get_coefficients function is building Pascal’s Triangle and returning it. Other code turns some row of the list into the terms to approximate the number e.
This works in reverse. It’s easy to generate Pascal’s Triangle, and we can use the terms of any row as the coefficients for a series of terms. The denominators, of course, are just increasing powers of n.
§
Two related factors complicate the complex version of this. Firstly, (1+i/n)ⁿ for any n greater than zero results in alternating imaginary terms. Secondly, those imaginary terms multiple to create negative coefficients. It turns out those alternate, too, but in pairs. Here’s an example:
It’s the expansion of (1+i/8)⁸. It has enough terms to illustrate the alternating pattern of imaginary terms (every other) and negative terms (two every two). We can see the effect of these by revisiting the five-term table from above:
| × 1 | 1 | 3i | -3 | -1i | |
|---|---|---|---|---|---|
| × ¼ | 1i | 3i2 | -3i | -1i2 | |
| numerators | 1 | 4i | -6 | -4i | 1 |
| denominators | 40 | 41 | 42 | 43 | 44 |
Note that the starting coefficients are (1, 3i, -3, -1i) — the negative values arising from previous shift operations (you can test this yourself). And every other term is imaginary. We multiply each term of this by (1+i/n), so again we add a copy multiplied by i/n (which effectively shifts it to the right because the powers of the denominators all go up by one).
Multiplying by i makes the real coefficients imaginary and the imaginary ones real. But in the latter case, it also negates the coefficient value because we end up with i², which is -1. Two of the coefficients in the second row of the table above (shown with i²) are negated, and their actual values are -3 and +1, respectively.
Now, here’s the (first part of the) algorithm for generating complex terms:
002|
003| def get_coefficients (degree=9):
004| “””Terms are tuples: (num,i_flag).”””
005| def shift (term):
006| # Return a shifted version of term…
007| if term[1]:
008| # Imaginary term: i² = -1…
009| return (–term[0], 0)
010| # Real term: becomes imaginary…
011| return (term[0], 1)
012|
013| def adder (term_a, term_b):
014| # Add the numeric terms…
015| c0 = term_a[0] + term_b[0]
016| # Set i_flag if set in A or B…
017| c1 = 1 if (term_a[1] or term_b[1]) else 0
018| return (c0, c1)
019|
020| def next_row (row):
021| ts0 = row + [(0,0)]
022| ts1 = [(0,0)] + [shift(t) for t in row]
023| # Return ts0+ts1…
024| return [adder(a,b) for a,b in zip(ts0,ts1)]
025|
026| table = [[(1,0)]]
027| for ix in range(degree):
028| row = next_row(table[–1])
029| table.append(row)
030| return table
031|
032| tab = get_coefficients(4)
033| for row in tab:
034| print(row)
035| print()
The function get_coefficients again only generates the list of rows of coefficients (same as found in Pascal’s Triangle, remember). But it treats the coefficients as tuples consisting of a value and an i_flag that determines whether the term is imaginary (true if non-zero). So, for example: (-3,0)=-3 and (5,1)=5i.
The output looks like this:
[(1, 0)] [(1, 0), (1, 1)] [(1, 0), (2, 1), (-1, 0)] [(1, 0), (3, 1), (-3, 0), (-1, 1)] [(1, 0), (4, 1), (-6, 0), (-4, 1), (1, 0)]
And it’s easy to see the negative coefficients and those flagged as imaginary. By ignoring those and taking the absolute values of the coefficients, we end up with a list of coefficients for calculating the real value of e (or equally, the rows of Pascal’s Triangle). That allows us to use the get_coefficients function for both real and complex values.
Below is the second part of the algorithm. There are several functions for doing different things with the coefficients we get from get_coefficients.
002|
003| def coeffs_2_nmbrs (coeffs, prn=False):
004| nmbrs = [abs(c[0]) for c in coeffs]
005| if prn:
006| print(‘,’.join([str(n) for n in nmbrs]))
007| print()
008| return nmbrs
009|
010| def coeffs_2_tuples (coeffs, prn=False):
011| denom = len(coeffs) – 1
012| terms = []
013| for ix,c in enumerate(coeffs):
014| term = (c[0], pow(denom,ix))
015| terms.append(term)
016| if prn:
017| for t in terms: print(t)
018| print()
019| return terms
020|
021| def coeffs_2_terms_real (coeffs, prn=False):
022| denom = len(coeffs) – 1
023| terms = []
024| for ix,c in enumerate(coeffs):
025| term = Fraction(abs(c[0]), pow(denom,ix))
026| terms.append(term)
027| real_value = float(sum(terms))
028| if prn:
029| for t in terms: print(term2str(t))
030| print(‘\n%.12f\n’ % real_value)
031| return (real_value, terms)
032|
033| def coeffs_2_terms_imag (coeffs, prn=False):
034| denom = len(coeffs) – 1
035| real_terms = []
036| imag_terms = []
037| for ix,c in enumerate(coeffs):
038| term = Fraction(c[0], pow(denom,ix))
039| if c[1]:
040| imag_terms.append(term)
041| else:
042| real_terms.append(term)
043| real_value = float(sum(real_terms))
044| imag_value = float(sum(imag_terms))
045| z = complex(real_value, imag_value)
046|
047| xr,xi = list(real_terms), list(imag_terms)
048| terms = []
049| while 0 < (len(xr)+len(xi)):
050| if 0 < len(xr): terms.append(xr.pop(0))
051| if 0 < len(xi): terms.append(xi.pop(0))
052| if prn:
053| for t in terms: print(term2str(t))
054| print()
055| return (z, real_terms, imag_terms)
056|
057| coeffs_list = get_coefficients(4, prn=False)
058| coeffs_row = coeffs_list[–1]
059|
060| coeffs_2_nmbrs(coeffs_row, prn=True)
061| coeffs_2_tuples(coeffs_row, prn=True)
062|
063| coeffs_2_terms_real(coeffs_row, prn=True)
064| coeffs_2_terms_imag(coeffs_row, prn=True)
The first function, coeffs_2_nmbrs, returns a list of the (absolute) numeric values for any given row of coefficients. These would be the values found in Pascal’s Triangle (or the numerators of the real series). This function allows us to use the coefficients generator to generate as much of the Triangle as we desire.
The second function, coeffs_2_tuples, returns the series terms as (numerator, denominator) tuples. I wanted this to see the actual fractions, not the reduced versions the next two functions return. I needed the full terms for that page on my other blog.
The third function, coeffs_2_terms_real, uses a row of coefficients to calculate the real value of e¹. (Which is 2.718281828459… but remember that approximation won’t be close without a lot of terms.)
The last function, coeffs_2_terms_imag, uses the coefficients to calculate the complex value ei. (Which is (0.540302…+0.841470…i), but again the approximation won’t be very close without out a lot of terms.)
The last lines just exercise the above functions. When run the output is:
1,4,6,4,1 (1, 1) (4, 4) (-6, 16) (-4, 64) (1, 256) +1.000000000e+00 = 1 +1.000000000e+00 = 1 +3.750000000e-01 = 3/8 +6.250000000e-02 = 1/16 +3.906250000e-03 = 1/256 2.441406250000 +1.000000000e+00 = 1 +1.000000000e+00 = 1 -3.750000000e-01 = -3/8 -6.250000000e-02 = -1/16 +3.906250000e-03 = 1/256
And this was exactly what I needed!
§
This ran long, so I’ll end abruptly. I’m always happy to answer questions about anything that wasn’t clear.
Ø
The Hard-Core Coder 
I probably should have picked something much shorter for getting back into blogging here. This one took me three days to write!
ATTENTION: The WordPress Reader strips the style information from posts, which can destroy certain important formatting elements. If you’re reading this in the Reader, I highly recommend (and urge) you to [A] stop using the Reader and [B] always read blog posts on their website.
This post is: Calculating the Number e (in Python)