Taylor approximation of exp(x) function for large x

Question

import numpy
from scipy.special import factorial

def Tn_exp(x):

    # Number of terms for Taylor Series
    MAX_N = 25 + numpy.abs(x)*3 
    
    # Coefficients of Taylor Series
    p = numpy.arange(MAX_N-1, -1, -1)

    # check the sign of x and do the calculation
    if x >= 0:
        Tn = 1/factorial(MAX_N)
        for coefficient in p:
            Tn = Tn * x + 1/factorial(coefficient)
    else:
        Tn = 1/factorial(MAX_N) * (-1)**MAX_N
        for coefficient in p:
            Tn = Tn * x + 1/factorial(coefficient) * (-1)**coefficient
        Tn = 1/Tn
    
    return Tn

I write a Python function to approximate the value of $\exp(x)$ by Taylor Series. It works well for small $x$ but it cannot deal with very large $x$ due to the problem of $x!$. It only shows inf for any $x > 170$ and will cause missing values for $\frac{1}{x!}$. Could someone give some ideas how to deal with this problem since I also consider the situation that $x$ is not an integer? Thank you very much.

Answer 1

Using: $$a_0=1, a_{n+1}=\frac{xa_n}{n+1},$$ we avoid dealing with huge factorials. Then we sum up enough $a_i.$

def Tn_exp(x):
    if x<0: 
        return 1/Tn_exp(-x)

    MaxN = int(3*x+25)
    a = 1.0
    sum = 1.0
    for i in  range(1,MaxN+1):
        a = a/i*x
        sum += a
    return sum

We could quit early if i>2x and a is very small compared to sum. That’s because the result of the terms will add up to less than $a.$

If you are not specifically trying to compute the Taylor series, just want $\exp$ implemented for larger $x,$ you won’t be able to do so using floating points, because the reason the built-in Python function can’t compute $\exp(710)$ is that the value overflows Python’s floating point.

If your only goal is computing $e^x$ for larger $x,$ it is easier to use $\log_{10}e$ and return pairs $(y,n)$ where $y$ is a real number, $1\leq y<10$ and $n$ is an integer, and the pair represents $y\cdot 10^n.$ Specifically:

import math

Loge10=math.log(10)
Log10e=1/Loge10

def better_exp(x):
    a = x*Log10e
    n = math.floor(a)
    b = a-n
    return (math.exp(b*Loge10),n)

This will be much faster.

The pairs $(y,n)$ are necessary because for larger $x,$ $e^x$ is larger than the maximum floating point value allowed for a Python floating point number.

Answer 2

Perhaps use the identity $${\displaystyle e^{x}=1+{\cfrac {x}{1-{\cfrac {x}{x+2-{\cfrac {2x}{x+3-{\cfrac {3x}{x+4-\ddots }}}}}}}}} $$ Due to Euler.