Can't calculate e^(-x) for "high" x

Question

This is my code:

def fact(y):
    if y == 0:
        fact=1
        return fact
    else:
        fact=1
        for k in range (1, y+1):
            fact = fact * k
        return fact

def e_negative_x(x):
    n=0
    numerical_precesion=1
    numerical_precesion_ideal= 10**(-8)
    while numerical_precesion > numerical_precesion_ideal:
        sum=0
        for i in range (0, n+1):
            ind=((-x)**i)/fact(i)
            if i> 0 & n-i == 0:
                ind_2=((-x)**(i-1))/fact(i-1)
                numerical_precesion_1 = ind_2 - ind
                numerical_precesion = abs(numerical_precesion_1)
            sum = sum + ind
        if   numerical_precesion > numerical_precesion_ideal:
            n += 1
        elif sum < 0:
            n += 1
    return sum

I tried to use this for x=0.1;1;10;15 as a exercise and I got "correct" for my precision but when I tried x=30 it got stuck at -8e-5 (incorrect answer). I tried to increase my terms but it was still stuck at -8e-5 with a different decimal places and after a big increase in terms didn't change at all.

Edit: It was -8e-5 which is wrong because at infinity this series tend to 0. I made a infinite loop (for infinite terms) and made it print out how many terms it had and the sum for that terms. I had 600 + terms before i shut down. At 89term i got stuck at -8e5 and after 117 i got stuck at -8.553016433669241e-05 till 600 + term.

Answer 1

How to do it better:

Well I would say summing up positive and negative terms of comparable size is prone to error since it leads to summing small and big numbers which is known to be problematic. But you can avoid it since exp(-30)=1/exp(30).

1/math.fsum(((30)**n/math.factorial(n) for n in range(0,100+1)))

is

9.357622968840174e-14

i.e. positve and small like you would expect. And

1/math.fsum(((30)**n/math.factorial(n) for n in range(0,100+1)))-np.exp(-30)

gives

-1.262177448353619e-29

So we get bascially the same as numpy does.

Where exactly is the error:

Here is a table of different ways to calculate exp(-30) approximately with the taylor series summing up N terms:

• Naive menas I just calculated the terms as floats and sumed them up.

• In the fractions column I did the same but used python's fraction instead of float division.

• Then I thought it would be nice to distinguish error in division/muliplication and in summing up. So I calculated each summand with fraction but turned them into floats before summing up.

• Finally I calculated the values of exp(30) nativly and then 1/exp(30) which I am calling the 1/e^30_trick.

            native     fractions  float(fraction)  1/e^30_trick
0    1.000000e+00  1.000000e+00     1.000000e+00  1.000000e+00
10   1.212548e+08  1.212548e+08     1.212548e+08  4.187085e-09
20   8.529171e+10  8.529171e+10     8.529171e+10  2.652040e-12
30   3.848426e+11  3.848426e+11     3.848426e+11  1.706501e-13
40   6.333654e+10  6.333654e+10     6.333654e+10  9.670058e-14
50   8.782292e+08  8.782292e+08     8.782292e+08  9.360412e-14
60   1.685584e+06  1.685584e+06     1.685584e+06  9.357627e-14
70   6.225246e+02  6.225247e+02     6.225246e+02  9.357623e-14
80   5.588481e-02  5.595324e-02     5.588481e-02  9.357623e-14
90  -6.697346e-05  1.459487e-06    -6.697346e-05  9.357623e-14
100 -6.843293e-05  1.276217e-11    -6.843293e-05  9.357623e-14
110 -6.843294e-05  9.361706e-14    -6.843294e-05  9.357623e-14
120 -6.843294e-05  9.357623e-14    -6.843294e-05  9.357623e-14

I was quite happy with the result. Since it shows both the naiv and even the float(fractions) way are terrible so the error must be in summing up the negative and positive terms. Also even though both the 1/exp(30) trick and the factions do converge to the "correct" value of 9.357623e-14. The former converges much faster than the latter.

Table Code:

series = pd.Series(np.arange(0,151,10))

@np.vectorize
def naiv(N):
    return math.fsum(((-30)**n/math.factorial(n) for n in range(0,N+1)))

@np.vectorize
def fractions(N):
    return float(sum(Fraction((-30)**n,math.factorial(n)) for n in range(0,N+1)))

@np.vectorize
def float_fractions(N):
    return math.fsum(float(Fraction((-30)**n,math.factorial(n))) for n in range(0,N+1))

@np.vectorize
def one_over_30_trick(N):
    return 1/math.fsum(((30)**n/math.factorial(n) for n in range(0,N+1)))

pd.DataFrame({"native":naiv(series),"fractions":fractions(series),"float(fraction)":float_fractions(series),"1/e^30_trick":one_over_30_trick(series)},index=series)