PROGRAM epsilon1
!
! Author:  Br. David Carlson
!
! Date:  December 31, 1999
!
! Revised:  May 2, 2004
!
! This program prints the value of the machine epsilon.

IMPLICIT NONE

REAL::epsvalue, x, xold, xnew, delta
INTEGER::n

xold = 1.0
delta = 1.0
xnew = xold + delta
n = 1

WRITE(*, 20)
20 FORMAT (' ', 'Program to investigate the machine epsilon')
WRITE (*, 50)
50 FORMAT (' ', T5, 'n', T15, 'delta', T35, 'xnew')

DO WHILE ((xold /= xnew) .AND. (n < 200))
   n = n + 1
   delta = delta * 0.7
   xnew = xold + delta
   WRITE (*, 100) n, delta, xnew
END DO

100 FORMAT (' ', T5, I3, T15, E17.9, T35, E17.9)
WRITE (*, 150) n, delta
150 FORMAT (' ', 'After ', I3, ' iterations we found epsilon to be about ',   &
   E17.9)

x = 1.0
epsvalue = EPSILON(x)

WRITE (*, 200) epsvalue
200 FORMAT (' ', 'Machine epsilon reported by EPSILON function is ', E17.9)

END PROGRAM