Mathematically, why does this Python program give such an inaccurate result for the Taylor series of $x \mapsto \exp\left(x\right)$ at $x=-40$?

Question

from math import pow, factorial
if __name__ == '__main__':
    # Parameters
    x = -40
    n_terms = 124
    sum = 0
    for i in range(124):
        sum += pow(x, i) / factorial(i)
    print(sum)

Output:

5.881161314775606

We know $e^{-40}$ is basically $0$, but it somehow evaluated to ~6. I'm not sure where the source of this error is coming from. Sure, floats aren't 100% precise, but they're still pretty damn accurate at ~16 digits (base 10) of precision and 64 bits. So the errors from computing each term in the taylor series should be insignificantly small. I find it very hard to believe the errors would accumulate to ~6 after a sum of only 124 terms.

Could you produce a mathematical argument for why an error of ~6 is plausible?

Note: According to the Lagrange Error Bound, the actual sum of the first 124 terms should be extremely close to 0

[abs(pow(-40, i)/factorial(i) - (-40)**i/factorial(i)) for i in range(120)] has maximum entry 2 attained when i is 35 and 39. The integer exponentiation is exact. The error caused by computing the power as float is significant. Even worse, your terms are large and alternate in sign so you get catastrophic cancellation — Matthew Towers, Jan 11 '25 at 13:42
@Beatnik Dopa, out of curiosity, how do you use the value you compute with this code? For example, is it a part of some other expression? — Konstantin Sidorov, Jan 11 '25 at 18:47
Instead of using huge powers and immense factorials, which can't be accurately represented, you don't need power or factorial functions at all, but a running term that is derived from the previous term, and added to the sum. — Weather Vane, Jan 12 '25 at 11:43

Konstantin Sidorov · Answer 1 · 2025-01-11T18:23:52.903

Sure, floats aren't 100% precise, but they're still pretty damn accurate at ~16 digits (base 10) of precision and 64 bits.

Unfortunately, your notion of "pretty damn accurate" is far from enough for your computation; to see why, observe that your 40th and 41st terms are going to be $(-40)^{40} / 40!≈1.5 × 10^{16}$ and $(-40)^{41} / 41!≈-1.4 × 10^{16}$, so the error in 16th digit is going to introduce an absolute error of $≈0.5 × (1.5 × 10^{16}) × 10^{-16}=0.75$ in a single term of your sum. That means the the absolute error of your largest terms is going to be comparable with the value your observe; add this for the other few terms, and the absolute error of 6 becomes not so surprising.

Edit. As Matthew Towers points out in his comment, what you have observed is a common phenomenon known as catastrophic cancellation, with many more examples to be found in the referenced article. You have a sign-alternating series with large terms, meaning that your code has to subtract large numbers; alas, just because you have good approximations to the subtracted terms does not mean you will get a good approximation to their difference.

Can you explain how you got that error from? I don't understand the calculation — Beatnik Dopa, Jan 11 '25 at 13:44
It is the same as your comment that floats are good to a relative error of $16$ decimal digits. If the number is of order $10^{16}$, the step is about $10^{-16}$ of that, or $1$. He uses a more accurate value, but the point is the same. — Ross Millikan, Jan 11 '25 at 14:46
Or in other words, the floating point error of $(a_1+...+a_n)$ is in first order $(|a_1|+...+|a_n|)\mu$, with $\mu=2^{-53}\approx 10^{-16}$ the machine constant. With that the error estimate for the explored calculation is $\exp(40)\mu\approx 26.1$. This is an upper bound, the error contributions have a random signs and will thus partially cancel. — Lutz Lehmann, Jan 11 '25 at 18:39

score 2 · Answer 2 · answered Jan 12 '25 at 14:01

The error in the Taylor series for exp(-40) is about the same as for exp(+40). The only difference is that the same error is huge compared to exp(-40) and small compared to exp(40).

Just print out all the terms and see how big they become to make it really obvious. This is basic numerical mathematics. And it has nothing whatsoever to do with Python.

Mathematically, why does this Python program give such an inaccurate result for the Taylor series of $x \mapsto \exp\left(x\right)$ at $x=-40$?

2 Answers2