!     Last change:  DOS  24 Jul 2000    8:26 pm
!                                       *** copyright 2000 ***
! *** filename macheps.f95              *** John F. Monahan **
!                                       **********************
      program macheps
!  compute machine epsilon in a particularly simple manner
!
      implicit none
      integer i,ieps,ione,id,ionep
      real one,onep,eps,d
!
  20  format(4x,'one in both formats',4x,f15.10,4x,z8)
  21  format(/2x,'epsilon',i4,'*(2**(',i3,'))',4x,f15.10,4x,z8)
  22  format(2x,'1 + eps ',4x,f15.10,4x,z8)
  23  format(/3x,'epsilon = 2**(',i3,')+ 2**(',i3,')',4x,f15.10,4x,z8)
!
!                   save in same spot
!  *** some compilers only permit integers written in z format ***
!
      equivalence (eps,ieps),(one,ione),(d,id),(onep,ionep)
!
!                   write the results here
      open( unit=6, file='macheps.out' )
!
!                   assign value to one (note equivalence)
      one = 1.
!                   assign initial value of 2**(-25)
      eps = scale( 1., -25)
!                   write out 1. in floating and internal
!                    representations 
      write(6,20) one,ione
!                   start with eps small and increase
      do i = 0,10
      d = float(i)*eps
!                   write out the epsilon
      write(6,21) i,-25,d,id
      onep = one + d
!                   write the sum if it's different from one
      if( onep .gt. one ) write(6,22) onep,ionep
      end do  ! loop on i
!                   try again but with 2**(-30)
      eps = scale( 1., -30)
!                   start with eps small and increase
      do i = 60,70
      d = float(i)*eps
!                   write out the epsilon
      write(6,21) i,-30,d,id
      onep = one + d
!                   write the sum if it's different from one
      if( onep .gt. one ) write(6,22) onep,ionep
      end do  ! loop on i
!                   try two more times
!                   first 1 + 2**(-24) + 2**(-47)   (different)
      eps = scale( 1., -24) + scale( 1., -47)
      write(6,23) -24,-47,eps,ieps
      onep = one + eps
      write(6,22) onep,ionep
!                   second ! + 2**(-24) + 2**(-48)  (same)
      eps = scale( 1., -24) + scale( 1., -48)
      write(6,23) -24,-48,eps,ieps
      onep = one + eps
      write(6,22) onep,ionep
      stop
      end program macheps
! *** end of filename macheps.f95       *********************