module mod_specialnumber
  implicit none

  real(4), parameter :: nan4 = transfer(Z'7FC00000', 0.e0)
  real(8), parameter :: nan8 = transfer(Z'7FF8000000000000', 0.d0)
  real(4), parameter :: pinf4 = transfer(Z'7F800000', 0.e0)
  real(4), parameter :: ninf4 = transfer(Z'FF800000', 0.e0)
  real(8), parameter :: pinf8 = transfer(Z'7FF0000000000000', 0.d0)
  real(8), parameter :: ninf8 = transfer(Z'FFF0000000000000', 0.d0)

  interface is_nan
    module procedure is_nan4
    module procedure is_nan8
  end interface

  interface is_inf
    module procedure is_inf4
    module procedure is_inf8
  end interface
contains

logical function is_nan4(f) result(is_nan)
  implicit none
  real(4), intent(in) :: f
  integer(4) :: n
  integer :: i

  n = transfer(f, 0_4)

  ! Exponent part
  is_nan = .true.
  do i = 23, 30
    if( .not. btest(n,i) )then
      is_nan = .false.
      return
    endif
  enddo

  ! Fraction part
  is_nan = .false.
  do i = 0, 22
    if( btest(n,i) )then
      is_nan = .true.
      return
    endif
  enddo
end function is_nan4

logical function is_nan8(f) result(is_nan)
  implicit none
  real(8), intent(in) :: f
  integer (8) :: n
  integer :: i

  n = transfer(f, 0_8)

  ! Exponent part
  is_nan = .true.
  do i = 52, 62
    if( .not. btest(n,i) )then
      is_nan = .false.
      return
    endif
  enddo

  ! Fraction part
  is_nan = .false.
  do i = 0, 51
    if( btest(n,i) )then
      is_nan = .true.
      return
    endif
  enddo
end function is_nan8

logical function is_inf4(f) result(is_inf)
  implicit none
  real(4), intent(in) :: f
  integer(4) :: n
  integer :: i

  n = transfer(f, 0_4)

  ! Exponent part
  is_inf = .true.
  do i = 23, 30
    if( .not. btest(n,i) )then
      is_inf = .false.
      return
    endif
  enddo

  ! Fraction part
  is_inf = .true.
  do i = 0, 22
    if( btest(n,i) )then
      is_inf = .false.
      return
    endif
  enddo
end function is_inf4

logical function is_inf8(f) result(is_inf)
  implicit none
  real(8), intent(in) :: f
  integer(8) :: n
  integer :: i

  n = transfer(f, 0_8)

  ! Exponent part
  is_inf = .true.
  do i = 52, 62
    if( .not. btest(n,i) )then
      is_inf = .false.
      return
    endif
  enddo

  ! Fraction part
  is_inf = .true.
  do i = 0, 51
    if( btest(n,i) )then
      is_inf = .false.
      return
    endif
  enddo
end function is_inf8

end module mod_nan


program main
  use mod_nan
  implicit none
  real(8), parameter :: nan8_ = transfer(-1_8, 0.d0)
  real(8), parameter :: pi = acos(-1.d0)

  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'pi   :', pi, &
         'is_nan:', is_nan(pi), 'is_inf:', is_inf(pi)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'nan8_:', nan8_, &
         'is_nan:', is_nan(nan8_), 'is_inf:', is_inf(nan8_)
  print*, ''

  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'nan4 :', nan4, 'is_nan:', is_nan(nan4), 'is_inf:', is_inf(nan4)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'nan8 :', nan8, 'is_nan:', is_nan(nan8), 'is_inf:', is_inf(nan8)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'pinf4:', pinf4, &
         'is_nan:', is_nan(pinf4), 'is_inf:', is_inf(pinf4)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'ninf4:', ninf4, &
         'is_nan:', is_nan(ninf4), 'is_inf:', is_inf(ninf4)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'pinf8:', pinf8, &
         'is_nan:', is_nan(pinf8), 'is_inf:', is_inf(pinf8)
  print"(1x,a,1x,f10.6,2(1x,a,1x,l1))",&
         'ninf8:', ninf8, &
         'is_nan:', is_nan(ninf8), 'is_inf:', is_inf(ninf8)
end program main