module utils
   use ISO_C_BINDING
   implicit none
   private
   public utils_init, ldmxcsr, stmxcsr, fstsw, fstcw, fldcw, &
      fclex, rdtsc, CPUID_type, cpuid
   integer(C_INT8_T) BAD_STUFF(224)
   data BAD_STUFF / &
      Z'31', Z'C0', Z'40', Z'75', Z'0A', Z'89', Z'4C', Z'24', &
      Z'08', Z'0F', Z'AE', Z'54', Z'24', Z'08', Z'C3', Z'0F', &
      Z'AE', Z'54', Z'24', Z'04', Z'C3', Z'90', Z'90', Z'90', &
      Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'40', Z'75', Z'0A', Z'0F', Z'AE', Z'5C', &
      Z'24', Z'08', Z'8B', Z'44', Z'24', Z'08', Z'C3', Z'50', &
      Z'0F', Z'AE', Z'1C', Z'24', Z'58', Z'C3', Z'90', Z'90', &
      Z'31', Z'C0', Z'DF', Z'E0', Z'C3', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'40', Z'75', Z'0A', Z'D9', Z'7C', Z'24', &
      Z'08', Z'66', Z'8B', Z'44', Z'24', Z'08', Z'C3', Z'50', &
      Z'D9', Z'3C', Z'24', Z'58', Z'C3', Z'90', Z'90', Z'90', &
      Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'40', Z'75', Z'0A', Z'66', Z'89', Z'4C', &
      Z'24', Z'08', Z'D9', Z'6C', Z'24', Z'08', Z'C3', Z'D9', &
      Z'6C', Z'24', Z'04', Z'C3', Z'90', Z'90', Z'90', Z'90', &
      Z'DB', Z'E2', Z'C3', Z'90', Z'90', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'40', Z'75', Z'0A', Z'0F', Z'31', Z'48', &
      Z'C1', Z'E2', Z'20', Z'48', Z'09', Z'D0', Z'C3', Z'0F', &
      Z'31', Z'C3', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', &
      Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'40', Z'75', Z'18', Z'53', Z'89', Z'D0', &
      Z'49', Z'87', Z'C8', Z'0F', Z'A2', Z'49', Z'90', Z'44', &
      Z'89', Z'00', Z'89', Z'58', Z'04', Z'89', Z'50', Z'08', &
      Z'89', Z'48', Z'0C', Z'5B', Z'C3', Z'53', Z'56', Z'8B', &
      Z'44', Z'24', Z'10', Z'8B', Z'4C', Z'24', Z'14', Z'0F', &
      Z'A2', Z'8B', Z'74', Z'24', Z'0C', Z'89', Z'06', Z'89', &
      Z'5E', Z'04', Z'89', Z'56', Z'08', Z'89', Z'4E', Z'0C', &
      Z'5E', Z'5B', Z'C3', Z'90', Z'90', Z'90', Z'90', Z'90' /
   type, bind(C) :: CPUID_type
      integer(C_INT32_T) eax
      integer(C_INT32_T) ebx
      integer(C_INT32_T) edx
      integer(C_INT32_T) ecx
   end type CPUID_type
   abstract interface
      subroutine ldmxcsr_template(source) bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT32_T), value :: source
      end subroutine ldmxcsr_template
      function stmxcsr_template() bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT32_T) :: stmxcsr_template
      end function stmxcsr_template
      function fstsw_template() bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT16_T) :: fstsw_template
      end function fstsw_template
      function fstcw_template() bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT16_T) :: fstcw_template
      end function fstcw_template
      subroutine fldcw_template(source) bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT16_T), value :: source
      end subroutine fldcw_template
      subroutine fclex_template() bind(C)
         use ISO_C_BINDING
         implicit none
      end subroutine fclex_template
      function rdtsc_template() bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT64_T) :: rdtsc_template
      end function rdtsc_template
      subroutine cpuid_template(result,eax,ecx) bind(C)
         use ISO_C_BINDING
         import CPUID_type
         implicit none
         type(CPUID_type) :: result
         integer(C_INT32_T), value :: eax
         integer(C_INT32_T), value :: ecx
      end subroutine cpuid_template
   end interface
   procedure(ldmxcsr_template), pointer :: ldmxcsr => NULL()
   procedure(stmxcsr_template), pointer :: stmxcsr => NULL()
   procedure(fstsw_template), pointer :: fstsw => NULL()
   procedure(fstcw_template), pointer :: fstcw => NULL()
   procedure(fldcw_template), pointer :: fldcw => NULL()
   procedure(fclex_template), pointer :: fclex => NULL()
   procedure(rdtsc_template), pointer :: rdtsc => NULL()
   procedure(cpuid_template), pointer :: cpuid => NULL()
   interface
      function VirtualAlloc(lpAddress, dwSize, flAllocationType, &
            flProtect) bind(C, name = 'VirtualAlloc')
         use ISO_C_BINDING
         implicit none
!GCC$ ATTRIBUTES STDCALL :: VirtualAlloc
         type(C_PTR) VirtualAlloc
         type(C_PTR), value :: lpAddress
         integer(C_SIZE_T), value :: dwSize
         integer(C_LONG), value :: flAllocationType
         integer(C_LONG), value :: flProtect
      end function VirtualAlloc
      function GetLastError() bind(C,name='GetLastError')
         use ISO_C_BINDING
         implicit none
!GCC$ ATTRIBUTES STDCALL :: GetLastError
         integer(C_LONG) GetLastError
      end function GetLastError
   end interface
   contains
      subroutine utils_init()
         type(C_PTR) address
         integer(C_INTPTR_T) temp
         integer(C_LONG) error
         integer(C_INT8_T), pointer :: temp_ptr(:)
         type(C_FUNPTR) fun_address
         logical :: first = .TRUE.

         if(.NOT.first) return
         first = .FALSE.
         address = VirtualAlloc(C_NULL_PTR, &
            size(BAD_STUFF,kind=C_SIZE_T),int(Z'1000',C_LONG), &
            int(Z'40',C_LONG))
         if(.NOT.C_ASSOCIATED(address)) then
            error = GetLastError()
            write(*,'(a,z0,a)') "Error Z'",error,"' allocating memory"
            stop
         end if
         call C_F_POINTER(address, temp_ptr, shape(BAD_STUFF))
         temp_ptr = BAD_STUFF
         temp = transfer(address, temp)
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, ldmxcsr)
         temp = temp+32
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, stmxcsr)
         temp = temp+24
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, fstsw)
         temp = temp+8
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, fstcw)
         temp = temp+32
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, fldcw)
         temp = temp+24
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, fclex)
         temp = temp+8
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, rdtsc)
         temp = temp+32
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, cpuid)
      end subroutine utils_init
end module utils

module nudge_FSW
   implicit none
   integer, parameter :: dp = kind(1.0d0)
   integer, parameter :: ep_preferred = selected_real_kind(18,4931)
   integer, parameter :: ep = &
      (1+sign(1,ep_preferred))/2*ep_preferred+ &
      (1-sign(1,ep_preferred))/2*dp
   private
   public ep, set_inexact
   contains
      subroutine set_inexact(x,y,z)
         real(ep) x, y, z
         x = y*z
      end subroutine set_inexact
end module nudge_FSW

program mxcsr_test
   use ISO_C_BINDING
   use utils
   use nudge_FSW
   implicit none
   integer(C_INT32_T) reg
   integer(C_INT64_T) ts1, ts2
   integer(C_INT16_T) sw, cw
   integer(C_INT32_T) eax, ecx
   type(CPUID_type) result
!   character(C_SIZEOF(result)) string
!   character(C_SIZEOF(CPUID_type(0,0,0,0))) string
   character(16) string
   real(ep) alpha
   real(ep) :: beta = 4*atan(1.0_ep)
   real(ep) :: gamma = sqrt(2.0_ep)

   write(*,'(a,i0,a)') 'This is a ', &
      bit_size(1_C_INTPTR_T),'-bit processor'
   call utils_init
   ts1 = rdtsc()
   reg = stmxcsr()
   ts2 = rdtsc()
   write(*,'(a,i0)') 'Time read from TSC: ',ts2-ts1
   write(*,'(a,z0.8)') 'Initial value of MXCSR: ',reg
   reg = ibclr(reg,9)
   write(*,'(a,z0.8)') 'Value to be stored in MXCSR: ',reg
   call ldmxcsr(reg)
   reg = stmxcsr()
   write(*,'(a,z0.8)') 'Value read from MXCSR: ',reg
   cw = fstcw()
   write(*,'(a,z0.4)') 'Initial value of FCW: ',cw
   cw = ibclr(cw, 2)
   write(*,'(a,z0.4)') 'Value to be stored in FCW: ',cw
   call fldcw(cw)
   cw = fstcw()
   write(*,'(a,z0.4)') 'Value read from FCW: ',cw
   sw = fstsw()
   write(*,'(a,z0.4)') 'Initial value of FSW: ',sw
   call set_inexact(alpha, beta, gamma)
   sw = fstsw()
   write(*,'(a,z0.4)') 'Value of FSW after inexact: ',sw
   call fclex
   sw = fstsw()
   write(*,'(a,z0.4)') 'Value of FSW after FCLEX: ',sw
   eax = 0
   ecx = 0
   call cpuid(result, eax, ecx)
   write(*,'(a,i0)') 'Max input value for basic CPUID: ', result%eax
   string = transfer(result,string)
   write(*,'(a,a)') 'Brand name: ',string(5:16)
!write(*,*) len(string)
write(*,'(4(z8.8:1x))') result
! Get ep parameters
   write(*,'(a,g14.6)') 'epsilon = ', epsilon(alpha)
!   write(*,'(a,g14.6)') 'huge = ', huge(alpha)
   write(*,*) 'huge = ', huge(alpha)
   write(*,'(a,i0)') 'precision = ', precision(alpha)
   write(*,'(a,i0)') 'range = ', range(alpha)
   write(*,'(a,i0)') 'sizeof = ', sizeof(alpha)
!   write(*,'(a,e15.6)') 'tiny = ', tiny(alpha)
   write(*,*) 'tiny = ', tiny(alpha)
end program mxcsr_test