Computing exponential function by Taylor Series without overflow

The original question is to write a Fortran program to compute the sum of the first 20 terms in the exponential equation (for x=1,2,3,4,5):
$$\sum_{n=0}^\infty \frac{x^n}{n!} = 1 + \frac{x^1}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \frac{x^5}{5!} + \ldots$$

My original approach was:

PROGRAM Exponential_Function

IMPLICIT NONE

INTEGER :: x,n,i
REAL :: f=0.0,fact

! Print Header
PRINT *, "           x     Sum"

Do x = 1,5,1
    Do n = 0,19,1

        ! calculate factorial
        fact =1.0
        Do i=1,n,1
            fact = fact*i
        End Do

        f=f+(x**n)/(fact)
    END Do

    PRINT *, x, f
    f=0.0
End Do

END PROGRAM Exponential_Function

The code could run, but only displayed results up to x=3:

Then, I tried another approach:

IMPLICIT NONE

INTEGER :: x,n
REAL :: f=1.0,term=1.0

! Print Header
PRINT *, "           x     Sum"

Do x = 1,5,1
    Do n = 1,19,1
        term = term * (REAL(x)/REAL(n))
        f = f + term

    END Do

    PRINT *, x, f
    term = 1.0
    f = 1.0 

End Do

and it gave the right results:

x     Sum
           1     2.71828
           2     7.38906
           3     20.0855
           4     54.5982
           5     148.413

I wonder if the second solution is the only correct solution to this question?
Is there a way to rectify the first approach to avoid the arithmetic overflow problen in Fortran?

Well the exponential function is map from reals to reals (usually denoted f:R→R), so I would have expected the use of reals for all variables (i.e., x and n). However, I'll base my answer using integers, as that is what you've used.

Using larger integers

Fortran's basic integer precision has a largest integer value of 2147483647, which is exceeded for 416=4294967296. However, you certainly can specify a larger integer type (see the link) and easily fit 419=274877906944. This would be as easy as specifying the kind:

integer :: x,n

would become

integer, parameter :: i_8 = selected_int_kind(15)
integer(i_8) :: x, n

And you easily get the values you desire.

Alternative implementation

Instead of specifying larger precision, you could express the Taylor-expansion of the exponential as
$$
e^x = 1 + x\cdot \left(1 + \frac{x}{2}\cdot\left(1 + \frac{x}{3}\cdot \left(\cdots\right) \right) \right)
$$
and write your function as,

f=1.0
do i=n-1,1,-1
   f = 1.0 + x * f / i
end do

(here n is your number of terms you want). Fortran will implicitly cast i and x as reals (since they are multiplying/dividing a real in f), so you don't need to worry about larger precision in your integers. This version could also fit easily into a function where you can have the user specify the n and x values:

real function my_exp(n, x) result(f)
   integer, intent(in) :: n, x
   f=1.0
   do i=n-1,1,-1
      f = 1.0 + x * f / i
   end do
end function

which is small enough that a simple contains will be fine for your program (i.e., using a module isn't needed).

Other odds-and-ends

Modern Fortran (well, since Fortran 90) does not require the use of capitals for keywords, so you don't need to have your program yell at you that it is printing something or that you're going to specify the types of all variables (i.e., you can use implicit none, real and integer). Most IDEs I know have support for Fortran syntax, so you'll know when you're using keywords.

While initializing f=1.0 in the main will do what you expect, you may want to break the habit of doing so now because you may (accidentally) carry it over to functions where it will have different behaviour.

StackExchange Code Review Q#146302, answer score: 4