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

IMPLICIT NONE

REAL::eps1, eps2, x

WRITE (*, 20)
20 FORMAT (' ', 'Program to investigate the machine epsilon')

CALL myepsilon(eps1)
WRITE (*, 150) eps1
150 FORMAT (' ', 'We found epsilon to be about ', E17.9)

x = 1.0
eps2 = EPSILON(x)

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

END PROGRAM


! Given:   Nothing.
! Task:    To compute the approximate value of the machine epsilon and to
!          print out intermediate results showing how we got this value.
! Return:  eps   This approximate value of the machine epsilon.
SUBROUTINE myepsilon(eps)

IMPLICIT NONE
REAL, INTENT(OUT)::eps

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

xold = 1.0
delta = 1.0
n = 0

WRITE (*, 50)
50 FORMAT (' ', T5, 'n', T15, 'delta', T35, 'xnew')

DO 
   n = n + 1
   xnew = xold + delta
   WRITE (*, 100) n, delta, xnew
   IF ((xnew == xold) .OR. (n == 200)) THEN
      EXIT
   END IF
   delta = delta * 0.75
END DO

eps = delta

100 FORMAT (' ', T5, I3, T15, E17.9, T35, E17.9)

END SUBROUTINE