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 *, "           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 *, "           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?

• It is long ago that I programmed in Fortran, but my guess would be that x**n is calculated in integer arithmetic and overflows (4^19 is too large for a 32-bit integer). REAL(x)**n might solve that. – Martin R Nov 6 '16 at 15:26
• Thanks! Fortran should include factorial as one of the intrinsic functions. – user3036719 Nov 6 '16 at 16:49
• Have you tried specifying a larger integer type? – Kyle Kanos Nov 7 '16 at 12:29
• No, eventually I used REAL(x) to solve the problem. – user3036719 Nov 12 '16 at 6:59

Well the exponential function is map from reals to reals (usually denoted f:RR), 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 4**16=4294967296. However, you certainly can specify a larger integer type (see the link) and easily fit 4**19=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.

• Thanks for the illustration and suggestions. Just have ocd to type keywords in capitals. – user3036719 Nov 20 '16 at 17:24
• I suppose I can understand wanting to differentiate keywords, but glad I can help! – Kyle Kanos Nov 20 '16 at 20:02