program double                                                      !test double precision numbers
  implicit none                                                     !enforce strong typing
  integer, parameter :: sp = selected_real_kind(p=precision(0.0))   !single precision kind
  integer, parameter :: dp = selected_real_kind(p=precision(0.0_sp)+1) !double precision kind
  integer, parameter :: k3 = selected_real_kind(p=precision(0.0_dp)+1)   !triple precision kind
!  integer, parameter :: k4 = selected_real_kind(p=precision(0.0_k3)+1)   !quadruple precision kind
!  integer, parameter :: k5 = selected_real_kind(p=precision(0.0_k4)+1)   !quintuple precision kind

  write(*,*)
  write(*,*) "Welcome to Double Precision demo"

  write(*,*)
  write(*,*) "KIND formerly was the number of (8-bit) bytes used to hold the number,"
  write(*,*) "  but now with larger word sizes (e.g., 64 bit), the KIND is "
  write(*,*) "  processor-dependent, and its value is arbitrary.   Caution."
  write(*,*) "  Here are the KIND values on my computer.  Yours might be different."
  write(*,*) "Single precision real kind is sp =", sp
  write(*,*) "Double precision real kind is dp =", dp
  write(*,*) "Triple precision real kind is tp =", k3
!  write(*,*) "Quad precision real kind is k4 =", k4

  write(*,*)
  write(*,*) "Precision is number of significant DECIMALS, range is +/- power of 10."
  write(*,*) "sp precision, range =", precision(0.0), range(0.0)
  write(*,*) "dp precision, range =", precision(0.0_dp), range(0.0_dp)
  write(*,*) "tp precision, range =", precision(0.0_k3), range(0.0_k3)
!  write(*,*) "qp precision, range =", precision(0.0_k4), range(0.0_k4)

  write(*,*)
  write(*,*) "Epsilon is the smallest resolvable difference between two decimal numbers."
  write(*,*) "sp epsilon =", epsilon(1.0)
  write(*,*) "dp epsilon =", epsilon(1.0_dp)
  write(*,*) "tp epsilon =", epsilon(1.0_k3)
!  write(*,*) "qp epsilon =", epsilon(1.0_k4)

  write(*,*)
  write(*,*) "Free format output shows one extra digit beyond the precision."
  write(*,*) "x = 1.0 in single precision", 1.0_sp
  write(*,*) "x = 1.0 in double precision", 1.0_dp
  write(*,*) "x = 1.0 in triple precision", 1.0_k3
!  write(*,*) "x = 1.0 in quad precision  ", 1.0_k4
  
  write(*,*)
endprogram double
