MODULE numerics
!-----------------------------------------------------------------------
! NAME
!     numerics
!
! DESCRIPTION
!     Management of numerical precision and special values (NaN and Inf)
!     throughout the code.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!     The intrinsic function storage_size() returns the storage size of
!     argument in bits.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use, intrinsic :: iso_fortran_env, only : int8, int16, int32, int64, real32, real64, real128
!~ use, intrinsic :: ieee_arithmetic

! DECLARATION
! -----------
implicit none

! PARAMETERS
! ----------
! Integers
integer, parameter :: ti = int8                   ! Tiny integer      = 8 bits
integer, parameter :: si = int16                  ! Short integer     = 16 bits
integer, parameter :: di = int32                  ! Standard integer  = 32 bits (default)
integer, parameter :: li = int64                  ! Long integer      = 64 bits
integer, parameter :: lli = selected_int_kind(30) ! Long long integer = 128 bits
integer, parameter :: wi = di                     ! Working integer

! Reals
integer, parameter :: sp = real32  ! Simple precision    = 32  bits = 1-bit sign + 8-bits  exponent + 23-bits  significand
integer, parameter :: dp = real64  ! Double precision    = 64  bits = 1-bit sign + 11-bits exponent + 52-bits  significand
integer, parameter :: qp = real128 ! Quadruple precision = 128 bits = 1-bit sign + 15-bits exponent + 112-bits significand
integer, parameter :: wp = dp      ! Working precision

! Logicals
! The precision is given by the compiler to be the same as the standard integer precision (32 bits) for memory reasons.
! It can be interesting to have smaller logical to reduce memory load.
! But be careful: definition is compiler-dependent and processor-dependent! So the following kinds are not guaranteed for logicals!
integer, parameter :: k1 = 1 ! kind = 1 should be 8  bits
integer, parameter :: k2 = 2 ! kind = 2 should be 16 bits
integer, parameter :: k4 = 4 ! kind = 4 should be 32 bits (default)

! Smallest strictly positive representable real number
real(sp), parameter :: smallest_sp = tiny(1._sp)
real(dp), parameter :: smallest_dp = tiny(1._dp)
real(qp), parameter :: smallest_qp = tiny(1._qp)
real(wp), parameter :: smallest_nb = tiny(1._wp)

! Largest representable number
integer(ti),  parameter :: largest_ti = huge(1_ti)
integer(si),  parameter :: largest_si = huge(1_si)
integer(di),  parameter :: largest_di = huge(1_di)
integer(li),  parameter :: largest_li = huge(1_li)
integer(lli), parameter :: largest_lli = huge(1_lli)
integer(wi),  parameter :: largest_wi = huge(1_wi)
real(sp),     parameter :: largest_sp = huge(1._sp)
real(dp),     parameter :: largest_dp = huge(1._dp)
real(qp),     parameter :: largest_qp = huge(1._qp)
real(wp),     parameter :: largest_nb = huge(1._wp)

! Machine epsilon: smallest positive real number such that 1 + epsilon > 1
real(sp), parameter :: minieps_sp = epsilon(1._sp)
real(dp), parameter :: minieps_dp = epsilon(1._dp)
real(qp), parameter :: minieps_qp = epsilon(1._qp)
real(wp), parameter :: minieps = epsilon(1._wp)
real(wp), parameter :: tol = 100._wp*minieps

! Minimum number of significant decimal digits
integer, parameter :: prec_sp = precision(1._sp)
integer, parameter :: prec_dp = precision(1._dp)
integer, parameter :: prec_qp = precision(1._qp)
integer, parameter :: prec = precision(1._wp)

!~ contains
!~ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!~ !=======================================================================
!~ FUNCTION is_finite(var) RESULT(flag)
!~ !-----------------------------------------------------------------------
!~ ! NAME
!~ !     is_finite
!~ !
!~ ! DESCRIPTION
!~ !     Returns whether array values are finite.
!~ !
!~ ! AUTHORS & DATE
!~ !     JB Clement, 01/2026
!~ !
!~ ! NOTES
!~ !     Finite values are neither NaN nor Infinite.
!~ !-----------------------------------------------------------------------

!~ ! DECLARATION
!~ ! -----------
!~ implicit none

!~ ! ARGUMENTS
!~ ! ---------
!~ real(dp), dimension(..), intent(in)  :: var

!~ ! LOCAL VARIABLES
!~ ! ---------------
!~ logical(k4) :: flag

!~ ! CODE
!~ ! ----
!~ flag = .false.
!~ if (all(ieee_is_finite(var))) flag = .true.

!~ END FUNCTION is_finite
!~ !=======================================================================

!~ !=======================================================================
!~ FUNCTION is_nan(var) RESULT(flag)
!~ !-----------------------------------------------------------------------
!~ ! NAME
!~ !     is_nan
!~ !
!~ ! DESCRIPTION
!~ !     Returns whether array values has a Not-a-Number (NaN).
!~ !
!~ ! AUTHORS & DATE
!~ !     JB Clement, 01/2026
!~ !
!~ ! NOTES
!~ !
!~ !-----------------------------------------------------------------------

!~ ! DECLARATION
!~ ! -----------
!~ implicit none

!~ ! ARGUMENTS
!~ ! ---------
!~ real(dp), dimension(..), intent(in)  :: var

!~ ! LOCAL VARIABLES
!~ ! ---------------
!~ logical(k4) :: flag

!~ ! CODE
!~ ! ----
!~ flag = .false.
!~ if (any(ieee_is_nan(var))) flag = .true.

!~ END FUNCTION is_nan
!~ !=======================================================================

!~ !=======================================================================
!~ FUNCTION is_normal(var) RESULT(flag)
!~ !-----------------------------------------------------------------------
!~ ! NAME
!~ !     is_normal
!~ !
!~ ! DESCRIPTION
!~ !     Returns whether array values is normal.
!~ !
!~ ! AUTHORS & DATE
!~ !     JB Clement, 01/2026
!~ !
!~ ! NOTES
!~ !     A normal value is finite, non-zero and not subnormal/denormal (very
!~ !     small finite number with reduced precision).
!~ !-----------------------------------------------------------------------

!~ ! DECLARATION
!~ ! -----------
!~ implicit none

!~ ! ARGUMENTS
!~ ! ---------
!~ real(dp), dimension(..), intent(in)  :: var

!~ ! LOCAL VARIABLES
!~ ! ---------------
!~ logical(k4) :: flag

!~ ! CODE
!~ ! ----
!~ flag = .false.
!~ if (all(ieee_is_normal(var))) flag = .true.

!~ END FUNCTION is_normal
!~ !=======================================================================

END MODULE numerics