program machineEpsilon

!   estimate machine epsilon
!   For double precision, expect eps = 2^-53
!   For single precision, expect eps = 2^-24

integer :: i, n, nd

!single precision


real :: eps = 1.0, res, res2

!double precision
!All these work

double precision :: epsd = 1.d0 ! literal 1.d0 is double precision

real*8 :: resd

integer, parameter :: DP = kind(1.d0)
!equivalent
!integer, parameter :: DP = selected_real_kind(14) !ensure 14 decimal digits of prec.
!integer, parameter :: DP = selected_real_kind(15,307) !15 figs, range of exponent of 307

real(kind=DP) :: res2d 


real(kind=DP), parameter :: MAX64 = huge(0.0_DP)
real(kind=DP), parameter :: MIN64 = tiny(0.0_DP)
real, parameter :: MAX32 = huge(0.0)
real, parameter :: MIN32 = tiny(0.0)

!For subtratcion example
real ::  f1,f2
real(kind=DP) ::  d1,d2,diff

n=30
nd=60

print *,"Single precision"
print *, "Largest 32-bit floating point number is", MAX32
print *, "Smallest 32-bit floating point number is", MIN32
print *, ""

do i=1,n
  res = 1.0 + eps
  res2 = res - 1.0
  print *, i-1,  eps, res, res2
  eps=eps*0.5
enddo 

print *,"Double precision"
print *, "Largest 64-bit floating point number is", MAX64
print *, "Smallest 64-bit floating point number is", MIN64
print *, ""

do i=1,nd
  resd = 1.d0 + epsd
! or 1.0_DP
  res2d = resd - 1.d0
  print *, i-1,  epsd, resd, res2d
  epsd=epsd*0.5d0
enddo



d2     =1.12345678902000000000000d0
d1     =1.12345678901234567890123d0
!d2-d1 =0.00000000000765432109877d0
diff   =7.65432109877d-12

f2=d2
f1=d1

res = f2-f1
resd = d2-d1

print *, "Subtraction example "
print *, "original numbers in double", d2, d1
print *, "original numbers in float", f2, f1
 
print *,"float and double differences", res, resd
print *,"relative error using doubles",(resd-diff)/diff

end program