One option would be to write the routines directly in assembly. This way the compiler is removed out of the picture.
@zaikunzhang out of curiosity I tried to play with your is_finite_sp implementation and got the following:
elemental pure function is_finite_sp(x) result(y)
use iso_fortran_env, only : wp=>real32
real(2*wp), parameter :: huge_2wp = 3.402823471D+38 ! increased a bit from huge(0.wp) in order to make the tests for huge pass
real(wp), intent(in) :: x
logical :: y
y = (huge_2wp - abs(x)>epsilon(0._wp))
end function is_finite_sp
With this change I get the following result:
True IS_FINITE: F F F F T T T T T T T T T T T T T T T
IS_FINITE: F F T T T T T T T T T T T T T T T T T
IEEE_IS_FINITE: T T T T T T T T T T T T T T T T T T T
The remaining issue here are the quite and signaling NaNs … here a small test Compiler Explorer what I see here is that with gfortran 13 the intrinsic isnan() returns False for both with -Ofast, but True with -O3 … No idea how to circumvent that one 
1 Like
Just a general comment, there are almost always some combination of compiler options that will cause a legitimate code to break. That is true for almost any compiler. So if the goal is to write fortran code that is bulletproof against all possible compiler option combinations, then that is doomed to failure.
In particular, there are options like -Ofast mentioned above, that tell the compiler to violate parentheses and to rearrange expressions in ways that will break almost any nontrivial fortran program. A programmer cannot expect his code to work in that situation. This is like trying to write portable C code with C compilers in the 1980s, it just wasn’t possible, and you would go crazy trying. C compilers of that era also did not respect parentheses, integer division i/j was allowed to either round up or down when an operand was negative, and the compiler could silently promote floats to doubles whenever it wanted. A programmer can’t do portable numerical computing under those conditions.
I think options like -Ofast are intended to be used only on small files, perhaps with a single subroutine, and perhaps with only a single line of code that does floating point arithmetic. If it works, then great, it speeds up that subroutine. But it was not really intended to compile a whole program, with millions of lines of floating point code, with those kinds of aggressive options.
I think the IEEE functions mostly work by simply comparing bit patterns, not by doing arithmetic on the operands and checking for hardware and/or software exceptions. The test for NaN is simply comparing to the NaN bit pattern, the test for Inf is comparing to the Inf bit pattern, and so on. The above code is trying to take the latter approach, which is not only inefficient, but nonportable since those exceptions depend on hardware characteristics, along with whatever compiler options are in effect at the time the code was compiled, along with whatever rounding and exception modes have been set at run time. And when different parts of the code are compiled with different sets of options, then the complications increase factorially.
It is a worthy goal to write bulletproof fortran code, but in practice it is somewhere between difficult to impossible.
3 Likes
Inspired by the PARANOIA program and issues similar to this we created a set of subroutines that started with some of the tests in paranoia that was extended over time to call routines that used numerically sensitive algorithms or exceptional values that we require to be called at the beginning of the program that generates a report including the compiler and compiler options, that also involves a key identifying that unique compilation for a program and the libraries used.
It is not available for release and is complicated so I am not suggesting that be done, but it might inspire a simpler similiar approach …
I think a simple subroutine making a unit test of these routines that you recommend in your documentation be called by programs calling the routines be added to your library that would identify when the tests do not work.
It is easy to even add a test of a saved value in a module to ensure the tests have been called and passed; although that can be unacceptable for performance reasons for procedures that are called heavily.
https://people.math.sc.edu/Burkardt/f_src/paranoia/paranoia.html
the approach is not for everyone, but if you have code that may be compiled in uncontrolled environments or ones where being paranoid is justified built-in and/or external compiler tests can help prevent surprises.
My biggest surprise is that you did not see similar issues with the other compilers.
Some of the examples of “incorrect” results above could easily be converted into a test procedure. If it is an acceptable performance hit the routines themselves can call the test on the first call and stop the program.
2 Likes
Here is a possible solution:
However, it contains non-standard Fortran code (sizeof(x) and integer(sizeof(x))).
This question inspired me to resurrect an old (incomplete) test program used to test behavior of special values and modernize it, as I THINK with a f2023 compiler with IEEE support most of the old results that were undefined now have defined standard-conforming behavior; but I could not
find if the MAX and MIN functions now specifically consider NAN and INF invalid values or if not,
what things like MAX(1.0,NAN) or MAX(INF,NAN) should produce. The first two compilers I tried produced different results with default optimization (did not even get to -Ofast). Is there a specified standard behavior or is that undefined. I could not find what I wanted in the standard to clarify that.
Just started modernizing the test so it is a WIP but ifort and gfortran already produced a few different results even without non-default optimization; which to me personally is more disturbing than when options that specifically allow for non-standard behavior are used.
test_special.f90
module M_nan
implicit none
private
public :: say_hello
public :: is_nan
public :: DIY_is_nan
contains
subroutine say_hello()
use, intrinsic :: iso_fortran_env, only : compiler_version
use, intrinsic :: iso_fortran_env, only : compiler_options
use, intrinsic :: IEEE_ARITHMETIC, only: IEEE_SUPPORT_DENORMAL, IEEE_SUPPORT_INF, IEEE_SUPPORT_NAN
character(len=:),allocatable :: string
integer :: ifound
write(*,'(a,/,1x,a)') 'This file was compiled by ', compiler_version()
write(*,'(a)') 'using the options'
string=compiler_options()
do
ifound=index(string,' -')
if(ifound.eq.0)then
write(*,'(1x,a)')string
exit
else
write(*,'(1x,a)')string(:ifound-1)
string=string(ifound+1:)
endif
enddo
write(*,*)'IEEE_SUPPORT_DENORMAL',IEEE_SUPPORT_DENORMAL(0.0)
write(*,*)'IEEE_SUPPORT_INF ',IEEE_SUPPORT_INF(0.0)
write(*,*)'IEEE_SUPPORT_NAN ',IEEE_SUPPORT_NAN (0.0)
end subroutine say_hello
elemental pure function is_nan(x)
!!use IEEE_EXCEPTIONS, only : ieee_support_nan ! is IEEE NaNs supported?
use, intrinsic :: IEEE_ARITHMETIC, only: IEEE_IS_NAN ! Determine if value is IEEE Not-a-Number.
use, intrinsic :: iso_fortran_env, only: int8, int16, int32, int64, real32, real64, real128
!$@ (#) M_units::is_nan(3f): determine if value is IEEE Not-a-Number
class(*), intent(in) :: x
logical :: is_nan
select type (x)
type is (real(kind=real32)); is_nan = ieee_is_nan(x)
type is (real(kind=real64)); is_nan = ieee_is_nan(x)
type is (real(kind=real128)); is_nan = ieee_is_nan(x)
type is (complex); is_nan = ieee_is_nan(real(x)) .and. ieee_is_nan(aimag(x))
!!type is (complex); is_nan=ieee_is_nan(x%re).and.ieee_is_nan(x%im)
end select
end function is_nan
elemental pure function DIY_is_nan(x)
!!use IEEE_EXCEPTIONS, only : ieee_support_nan ! is IEEE NaNs supported?
use, intrinsic :: iso_fortran_env, only: int8, int16, int32, int64, real32, real64, real128
!$@ (#) M_units::DIY_is_nan(3f): determine if value is IEEE Not-a-Number
class(*), intent(in) :: x
logical :: DIY_is_nan
select type (x)
type is (real(kind=real32)); DIY_is_nan = x.ne.x
type is (real(kind=real64)); DIY_is_nan = x.ne.x
type is (real(kind=real128)); DIY_is_nan = x.ne.x
type is (complex); DIY_is_nan = x.ne.x
end select
end function DIY_is_nan
end module M_nan
program main
use, intrinsic :: ieee_arithmetic
use M_nan, only: say_hello
use M_nan, only: is_nan
!use M_nan, only: is_nan=>DIY_is_nan
implicit none
real :: r1, r2, A(4), x, y, nan, inf
character(len=*),parameter :: g = '(*(g0,1x))'
logical,parameter :: T = .true., F = .false.
logical :: v(size(A))
! the cross platform way to get NaN on modern compilers
nan = ieee_value(0.0, ieee_quiet_nan)
inf = ieee_value(0.0, ieee_positive_inf)
call say_hello()
write(*,'(g0,g0,1x,z0)')'INFINITY=',nan,nan
write(*,'(g0,g0,1x,z0)')'INFINITY=',inf,inf
A = [-1.0, 0.0, 1.0, 2.0]
r1 = 1.0; call printit(r1, 'IS_NAN(1.0) = ', F)
r1 = 1.0; r2 = 0.0; call printit(r1/r2, 'IS_NAN(1.0/0.0) = ', F)
r1 = 0.0; r2 = 0.0; call printit(r1/r2, 'IS_NAN(0.0/0.0) = ', T)
r1 = 0.0; call printit(r1**0, 'IS_NAN(0^0) = ', F)
r1 = 0.0; r2 = 0.0; call printit(r1**r2, 'IS_NAN(0^0) = ', F)
r1 = -2.0; call printit(acos(r1), 'IS_NAN(acos(-2)) = ', T)
r1 = 1000.0; call printit(exp(r1), 'IS_NAN(exp(1000)) = ', F)
r1 = 0.0; call printit(log(r1), 'IS_NAN(log(0)) = ', F)
r1 = -1.0; call printit(sqrt(r1), 'IS_NAN(sqrt(-1)) = ', T)
write (*, g) 'For X=', A
call printits((acos(A)), [F, F, F, T], 'ACOS(X)')
call printits((asin(A)), [F, F, F, T], 'ASIN(X)')
call printits((atan(A)), [F, F, F, F], 'ATAN(X)')
call printits((log(A)), [T, F, F, F], 'LOG(X)')
call printits((sqrt(A)), [T, F, F, F], 'SQRT(X)')
call printits((1.0/(A)), [F, F, F, F], '1/X')
! Fortran handling of NaN vis-a-vis max and min
if(.not.ieee_is_nan(nan)) then
write(*,g) 'NaN: ',nan, 'NaN not working correctly'
endif
A = [0., 1., 2., nan]
print *, 'maximum of ',A, 'is', maxval(A)
if(maxval(A) /= 2) write(*,*) 'maxval() not working with NaN'
x = nan
y = 1.0
print *, 'maximum of',x,'and',y,'is',max(x,y)
if (max(x,y) /= 1) write(*,*) 'max() not working with NaN'
print *,'max of inf() and NaN is', max(inf, nan)
write(*,*)'inf>nan',inf>nan
write(*,*)'nam>inf',nan>inf
write(*,*)'inf>inf',inf>inf
write(*,*)'nan>nan',nan>nan
write(*,*)'inf<nan',inf<nan
write(*,*)'nam<inf',nan<inf
write(*,*)'inf<inf',inf<inf
write(*,*)'nan<nan',nan<nan
write(*,*)'inf==nan',inf==nan
write(*,*)'nam==inf',nan==inf
write(*,*)'inf==inf',inf==inf
write(*,*)'nan==nan',nan==nan
contains
subroutine printits(arg, expected, label)
real, intent(in) :: arg(:)
logical, intent(in) :: expected(:)
character(len=*), intent(in) :: label
write(*,g)merge('PASSED','FAILED',all(is_nan(arg).eqv.expected)),label,'expected',expected,'got',is_nan(arg),'for arg',arg
end subroutine printits
subroutine printit(arg, label, expected, notes)
real, intent(in) :: arg
character(len=*), intent(in) :: label
logical, intent(in) :: expected
character(len=*), intent(in), optional :: notes
if (present(notes)) then
write(*,g)merge('PASSED','FAILED',is_nan(arg).eqv.expected),label,'expected',expected,'got',is_nan(arg),'for arg',arg,notes
else
write(*,g)merge('PASSED','FAILED',is_nan(arg).eqv.expected),label,'expected',expected,'got',is_nan(arg),'for arg',arg
endif
end subroutine printit
end program main
1 Like
Hi @urbanjost ,
Thank you for raising this point. I have also observed similar problems during the development of PRIMA. As far as I recall, the behavior is unspecified in Fortran standards. This is one of the reasons I would like to have reliable tests of Inf and NaN.
In addition, I also found that maxval([a, b]) and max(a, b) may produce different results if a or b is NaN even without aggressive optimization. But I do not have a minimal example for the moment.
True. If my memory is correct, the behavior is unspecified in the standard. However, Inf and NaN are quite common in real-life scientific computing except for perfectly conditioned problems.
Maybe the non-standard sizeof can be converted to F2008 standard code using c_sizeof or storage_size? Note the latter gives the size in bits, the former in bytes.
Maybe. However, it is still nonstandard to relate the storage size with the kind, i.e., integer(storage_size(x)) :: y.
@zaikunzhang I managed to get all your tests to pass with gfortran 13 -Ofast, you can check and run the modifications here: Compiler Explorer
Full code
module consts_mod
implicit none
private
public :: SP, DP
integer, parameter :: SP = kind(0.0)
integer, parameter :: DP = kind(0.0D0)
end module consts_mod
module inf_mod
implicit none
private
public :: is_finite, is_inf, is_posinf, is_neginf
interface is_finite
module procedure is_finite_sp, is_finite_dp
end interface is_finite
interface is_posinf
module procedure is_posinf_sp, is_posinf_dp
end interface is_posinf
interface is_neginf
module procedure is_neginf_sp, is_neginf_dp
end interface is_neginf
interface is_inf
module procedure is_inf_sp, is_inf_dp
end interface is_inf
contains
elemental pure function is_finite_sp(x) result(y)
use consts_mod, only : wp=>SP
real(2*wp), parameter :: huge_2wp = 3.40282347D+38
real(wp), intent(in) :: x
logical :: y
y = huge_2wp>=abs(x) .and. -abs(x)>=-huge_2wp
end function is_finite_sp
elemental pure function is_finite_dp(x) result(y)
use consts_mod, only : wp=>DP
real(2*wp), parameter :: huge_2wp = 1.7976931348623158d+308
real(2*wp), parameter :: z = 0.0000000000000002d+308
real(wp), intent(in) :: x
logical :: y
y = huge_2wp + z>=abs(x) .and. -abs(x)>=-huge_2wp
end function is_finite_dp
elemental pure function is_inf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(2*wp), parameter :: huge_2wp = 3.402823471D+38
real(wp), intent(in) :: x
logical :: y
y = ( tiny(0._wp) >= huge_2wp-abs(x))
end function is_inf_sp
elemental pure function is_inf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(2*wp), parameter :: huge_2wp = 1.7976931348623158d+308
real(2*wp), parameter :: z = 0.0000000000000002d+308
real(wp), intent(in) :: x
logical :: y
y = (abs(x) > huge_2wp + z )
end function is_inf_dp
elemental pure function is_posinf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(2*wp), parameter :: huge_2wp = 3.402823471D+38
real(wp), intent(in) :: x
logical :: y
y = (abs(x) > huge_2wp) .and. (x > 0._wp)
end function is_posinf_sp
elemental pure function is_posinf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(2*wp), parameter :: huge_2wp = 1.7976931348623158d+308
real(2*wp), parameter :: z = 0.0000000000000002d+308
real(wp), intent(in) :: x
logical :: y
y = (abs(x) > huge_2wp + z) .and. (x > 0._wp)
end function is_posinf_dp
elemental pure function is_neginf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(2*wp), parameter :: huge_2wp = 3.402823471D+38
real(wp), intent(in) :: x
logical :: y
y = (abs(x) > huge_2wp) .and. (x < 0._wp)
end function is_neginf_sp
elemental pure function is_neginf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(2*wp), parameter :: huge_2wp = 1.7976931348623158d+308
real(2*wp), parameter :: z = 0.0000000000000002d+308
real(wp), intent(in) :: x
logical :: y
y = (abs(x) > huge_2wp + z) .and. (x < 0._wp)
end function is_neginf_dp
end module inf_mod
module ieee_infnan_mod
use, non_intrinsic :: consts_mod, only : SP, DP
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_is_finite
implicit none
private
public :: ieee_is_inf, ieee_is_posinf, ieee_is_neginf
interface ieee_is_inf
module procedure ieee_is_inf_sp, ieee_is_inf_dp
end interface ieee_is_inf
interface ieee_is_posinf
module procedure ieee_is_posinf_sp, ieee_is_posinf_dp
end interface ieee_is_posinf
interface ieee_is_neginf
module procedure ieee_is_neginf_sp, ieee_is_neginf_dp
end interface ieee_is_neginf
contains
elemental pure function ieee_is_inf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x))
end function ieee_is_inf_sp
elemental pure function ieee_is_inf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x))
end function ieee_is_inf_dp
elemental pure function ieee_is_posinf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x > 0)
end function ieee_is_posinf_sp
elemental pure function ieee_is_posinf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x > 0)
end function ieee_is_posinf_dp
elemental pure function ieee_is_neginf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x < 0)
end function ieee_is_neginf_sp
elemental pure function ieee_is_neginf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x < 0)
end function ieee_is_neginf_dp
end module ieee_infnan_mod
module infnan_mod
use inf_mod, only : is_finite, is_inf, is_posinf, is_neginf
implicit none
private
public :: is_finite, is_inf, is_posinf, is_neginf, is_nan
interface is_nan
module procedure is_nan_sp, is_nan_dp
end interface is_nan
contains
elemental pure function is_nan_sp(x) result(y)
use consts_mod, only : SP
implicit none
real(2*SP), parameter :: huge_2wp = 3.402823471D+38
real(SP), intent(in) :: x
logical :: y
y = .not.(abs(x)>huge_2wp) .and. .not.(abs(x)<huge_2wp)
end function is_nan_sp
elemental pure function is_nan_dp(x) result(y)
use consts_mod, only : DP
implicit none
real(2*dp), parameter :: huge_2wp = 1.7976931348623156d+308
real(DP), intent(in) :: x
logical :: y
y = .not.(abs(x)>huge_2wp) .and. .not.(abs(x)<huge_2wp)
end function is_nan_dp
end module infnan_mod
module test_infnan_mod
implicit none
private
public :: test_infnan
contains
function test_infnan() result(success)
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_is_finite, ieee_value, &
& ieee_quiet_nan, ieee_signaling_nan, ieee_positive_inf, ieee_negative_inf, &
& ieee_positive_normal, ieee_negative_normal, ieee_positive_denormal, ieee_negative_denormal, &
& ieee_positive_zero, ieee_negative_zero
use, non_intrinsic :: consts_mod, only : SP, DP
use, non_intrinsic :: ieee_infnan_mod, only : ieee_is_inf, ieee_is_posinf, ieee_is_neginf
use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_inf, is_posinf, is_neginf
use iso_c_binding
implicit none
logical :: success
real(SP) :: PINF_SP, NINF_SP, QNAN_SP, SNAN_SP, PNORMAL_SP, NNORMAL_SP, PDENORMAL_SP, NDENORMAL_SP, &
& PZERO_SP, NZERO_SP
real(DP) :: PINF_DP, NINF_DP, QNAN_DP, SNAN_DP, PNORMAL_DP, NNORMAL_DP, PDENORMAL_DP, NDENORMAL_DP, &
& PZERO_DP, NZERO_DP
real(SP) :: x(19)
real(DP) :: y(size(x))
logical :: check_nan(size(x))
logical :: check_finite(size(x))
logical :: check_inf(size(x))
logical :: check_posinf(size(x))
logical :: check_neginf(size(x))
logical :: success_SP
logical :: success_DP
PINF_SP = ieee_value(1.0, ieee_positive_inf)
NINF_SP = ieee_value(1.0, ieee_negative_inf)
QNAN_SP = ieee_value(1.0, ieee_quiet_nan)
SNAN_SP = ieee_value(1.0, ieee_signaling_nan)
PNORMAL_SP = ieee_value(1.0, ieee_positive_normal)
NNORMAL_SP = ieee_value(1.0, ieee_negative_normal)
PDENORMAL_SP = ieee_value(1.0, ieee_positive_denormal)
NDENORMAL_SP = ieee_value(1.0, ieee_negative_denormal)
PZERO_SP = ieee_value(1.0, ieee_positive_zero)
NZERO_SP = ieee_value(1.0, ieee_negative_zero)
PINF_DP = ieee_value(1.0D0, ieee_positive_inf)
NINF_DP = ieee_value(1.0D0, ieee_negative_inf)
QNAN_DP = ieee_value(1.0D0, ieee_quiet_nan)
SNAN_DP = ieee_value(1.0D0, ieee_signaling_nan)
PNORMAL_DP = ieee_value(1.0D0, ieee_positive_normal)
NNORMAL_DP = ieee_value(1.0D0, ieee_negative_normal)
PDENORMAL_DP = ieee_value(1.0D0, ieee_positive_denormal)
NDENORMAL_DP = ieee_value(1.0D0, ieee_negative_denormal)
PZERO_DP = ieee_value(1.0D0, ieee_positive_zero)
NZERO_DP = ieee_value(1.0D0, ieee_negative_zero)
x = [PINF_SP, NINF_SP, &
& QNAN_SP, SNAN_SP, &
& PNORMAL_SP, NNORMAL_SP, &
& PDENORMAL_SP, NDENORMAL_SP, &
& PZERO_SP, NZERO_SP, &
& huge(0.0), -huge(0.0), &
& epsilon(0.0), -epsilon(0.0), &
& tiny(0.0), -tiny(0.0), &
& 0.0, 1.0, -1.0]
y = [PINF_DP, NINF_DP, &
& QNAN_DP, SNAN_DP, &
& PNORMAL_DP, NNORMAL_DP, &
& PDENORMAL_DP, NDENORMAL_DP, &
& PZERO_DP, NZERO_DP, &
& huge(0.0D0), -huge(0.0D0), &
& epsilon(0.0D0), -epsilon(0.0D0), &
& tiny(0.0D0), -tiny(0.0D0), &
& 0.0D0, 1.0D0, -1.0D0]
check_finite = [.false., .false., &
& .false., .false., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., .true.]
check_inf = [.true., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_posinf = [.true., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_neginf = [.false., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_nan = [.false., .false., &
& .true., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
write (*, '(/(1A))') 'Single-precision tests:'
if (.not. all(is_finite(x) .eqv. check_finite)) then
write (*, *) 'True IS_FINITE: ', check_finite
write (*, *) 'IS_FINITE: ', is_finite(x)
write (*, *) 'IEEE_IS_FINITE: ', ieee_is_finite(x)
end if
if (.not. all(is_inf(x) .eqv. check_inf)) then
write (*, *) 'True IS_INF: ', check_inf
write (*, *) 'IS_INF: ', is_inf(x)
write (*, *) 'IEEE_IS_INF: ', ieee_is_inf(x)
end if
if (.not. all(is_posinf(x) .eqv. check_posinf)) then
write (*, *) 'True IS_POSINF: ', check_posinf
write (*, *) 'IS_POSINF: ', is_posinf(x)
write (*, *) 'IEEE_IS_POSINF: ', ieee_is_posinf(x)
end if
if (.not. all(is_neginf(x) .eqv. check_neginf)) then
write (*, *) 'True IS_NEGINF: ', check_neginf
write (*, *) 'IS_NEGINF: ', is_neginf(x)
write (*, *) 'IEEE_IS_NEGINF: ', ieee_is_neginf(x)
end if
if (.not. all(is_nan(x) .eqv. check_nan)) then
write (*, *) 'True IS_NAN: ', check_nan
write (*, *) 'IS_NAN: ', is_nan(x)
write (*, *) 'IEEE_IS_NAN: ', ieee_is_nan(x)
end if
success_SP = all(is_finite(x) .eqv. check_finite) .and. all(is_inf(x) .eqv. check_inf) .and. &
& all(is_posinf(x) .eqv. check_posinf) .and. all(is_neginf(x) .eqv. check_neginf) .and. &
& all(is_nan(x) .eqv. check_nan)
if (.not. success_SP) then
write (*, '(1A)') 'Some tests failed. The data is'
write (*, '(1P, 20E8.0)') x
else
write (*, '(1A)') 'All tests succeed.'
end if
write (*, '(/(1A))') 'Double-precision tests:'
if (.not. all(is_finite(y) .eqv. check_finite)) then
write (*, *) 'True IS_FINITE: ', check_finite
write (*, *) 'IS_FINITE: ', is_finite(y)
write (*, *) 'IEEE_IS_FINITE: ', ieee_is_finite(y)
end if
if (.not. all(is_inf(y) .eqv. check_inf)) then
write (*, *) 'True IS_INF: ', check_inf
write (*, *) 'IS_INF: ', is_inf(y)
write (*, *) 'IEEE_IS_INF: ', ieee_is_inf(y)
end if
if (.not. all(is_posinf(y) .eqv. check_posinf)) then
write (*, *) 'Truth IS_POSINF: ', check_posinf
write (*, *) 'IS_POSINF: ', is_posinf(y)
write (*, *) 'IEEE_IS_POSINF: ', ieee_is_posinf(y)
end if
if (.not. all(is_neginf(y) .eqv. check_neginf)) then
write (*, *) 'True IS_NEGINF: ', check_neginf
write (*, *) 'IS_NEGINF: ', is_neginf(y)
write (*, *) 'IEEE_IS_NEGINF: ', ieee_is_neginf(y)
end if
if (.not. all(is_nan(y) .eqv. check_nan)) then
write (*, *) 'True IS_NAN: ', check_nan
write (*, *) 'IS_NAN: ', is_nan(y)
write (*, *) 'IEEE_IS_NAN: ', ieee_is_nan(y)
end if
success_DP = all(is_finite(y) .eqv. check_finite) .and. all(is_inf(y) .eqv. check_inf) .and. &
& all(is_posinf(y) .eqv. check_posinf) .and. all(is_neginf(y) .eqv. check_neginf) .and. &
& all(is_nan(y) .eqv. check_nan)
if (.not. success_DP) then
write (*, '(1A)') 'Some tests failed. The data is'
write (*, '(1P, 20D8.0)') y
else
write (*, '(1A)') 'All tests succeed.'
end if
write (*, *) ''
success = success_SP .and. success_DP
end function test_infnan
end module test_infnan_mod
program test
use, non_intrinsic :: test_infnan_mod, only : test_infnan
if (.not. test_infnan()) then
stop 1
end if
end program test
Something that took me a while to get was testing for NaNs. Found out that it was better to test for the following simultaneous condition .not.(abs(x)>huge_2wp) .and. .not.(abs(x)<huge_2wp) it turns out that for both types of NaNs, both conditions will simultanously return False.
EDIT 2023-09-13
Full Code
module consts_mod
implicit none
private
public :: SP, DP
integer, parameter :: SP = kind(0.0)
integer, parameter :: DP = kind(0.0D0)
end module consts_mod
module inf_mod
implicit none
private
public :: is_finite, is_inf, is_posinf, is_neginf
interface is_finite
module procedure is_finite_sp, is_finite_dp
end interface is_finite
interface is_posinf
module procedure is_posinf_sp, is_posinf_dp
end interface is_posinf
interface is_neginf
module procedure is_neginf_sp, is_neginf_dp
end interface is_neginf
interface is_inf
module procedure is_inf_sp, is_inf_dp
end interface is_inf
contains
elemental pure function is_finite_sp(x) result(y)
use consts_mod, only : wp=>SP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = huge_2wp>=abs(x) .and. -abs(x)>=-huge_2wp
end function is_finite_sp
elemental pure function is_finite_dp(x) result(y)
use consts_mod, only : wp=>DP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = huge_2wp >=abs(x) .and. -abs(x)>=-huge_2wp
end function is_finite_dp
elemental pure function is_inf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = ( tiny(0._wp) >= huge_2wp-abs(x))
end function is_inf_sp
elemental pure function is_inf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = (abs(x) > huge_2wp )
end function is_inf_dp
elemental pure function is_posinf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = (abs(x) > huge_2wp) .and. (x > 0._wp)
end function is_posinf_sp
elemental pure function is_posinf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = (abs(x) > huge_2wp) .and. (x > 0._wp)
end function is_posinf_dp
elemental pure function is_neginf_sp(x) result(y)
use consts_mod, only : wp=>SP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = (abs(x) > huge_2wp) .and. (x < 0._wp)
end function is_neginf_sp
elemental pure function is_neginf_dp(x) result(y)
use consts_mod, only : wp=>DP
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = (abs(x) > huge_2wp) .and. (x < 0._wp)
end function is_neginf_dp
end module inf_mod
module ieee_infnan_mod
use, non_intrinsic :: consts_mod, only : SP, DP
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_is_finite
implicit none
private
public :: ieee_is_inf, ieee_is_posinf, ieee_is_neginf
interface ieee_is_inf
module procedure ieee_is_inf_sp, ieee_is_inf_dp
end interface ieee_is_inf
interface ieee_is_posinf
module procedure ieee_is_posinf_sp, ieee_is_posinf_dp
end interface ieee_is_posinf
interface ieee_is_neginf
module procedure ieee_is_neginf_sp, ieee_is_neginf_dp
end interface ieee_is_neginf
contains
elemental pure function ieee_is_inf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x))
end function ieee_is_inf_sp
elemental pure function ieee_is_inf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x))
end function ieee_is_inf_dp
elemental pure function ieee_is_posinf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x > 0)
end function ieee_is_posinf_sp
elemental pure function ieee_is_posinf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x > 0)
end function ieee_is_posinf_dp
elemental pure function ieee_is_neginf_sp(x) result(y)
real(SP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x < 0)
end function ieee_is_neginf_sp
elemental pure function ieee_is_neginf_dp(x) result(y)
real(DP), intent(in) :: x
logical :: y
y = (.not. ieee_is_nan(x)) .and. (.not. ieee_is_finite(x)) .and. (x < 0)
end function ieee_is_neginf_dp
end module ieee_infnan_mod
module infnan_mod
use inf_mod, only : is_finite, is_inf, is_posinf, is_neginf
implicit none
private
public :: is_finite, is_inf, is_posinf, is_neginf, is_nan
interface is_nan
module procedure is_nan_sp, is_nan_dp
end interface is_nan
contains
elemental pure function is_nan_sp(x) result(y)
use consts_mod, only : wp=>SP
implicit none
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = .not.(abs(x)>huge_2wp) .and. .not.(abs(x)<huge_2wp)
end function is_nan_sp
elemental pure function is_nan_dp(x) result(y)
use consts_mod, only : wp=>DP
implicit none
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
logical :: y
y = .not.(abs(x)>huge_2wp) .and. .not.(abs(x)<huge_2wp)
end function is_nan_dp
end module infnan_mod
module test_infnan_mod
implicit none
private
public :: test_infnan
contains
function test_infnan() result(success)
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_is_finite, ieee_value, &
& ieee_quiet_nan, ieee_signaling_nan, ieee_positive_inf, ieee_negative_inf, &
& ieee_positive_normal, ieee_negative_normal, ieee_positive_denormal, ieee_negative_denormal, &
& ieee_positive_zero, ieee_negative_zero
use, non_intrinsic :: consts_mod, only : SP, DP
use, non_intrinsic :: ieee_infnan_mod, only : ieee_is_inf, ieee_is_posinf, ieee_is_neginf
use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_inf, is_posinf, is_neginf
use iso_c_binding
implicit none
logical :: success
real(SP) :: PINF_SP, NINF_SP, QNAN_SP, SNAN_SP, PNORMAL_SP, NNORMAL_SP, PDENORMAL_SP, NDENORMAL_SP, &
& PZERO_SP, NZERO_SP
real(DP) :: PINF_DP, NINF_DP, QNAN_DP, SNAN_DP, PNORMAL_DP, NNORMAL_DP, PDENORMAL_DP, NDENORMAL_DP, &
& PZERO_DP, NZERO_DP
real(SP) :: x(19)
real(DP) :: y(size(x))
logical :: check_nan(size(x))
logical :: check_finite(size(x))
logical :: check_inf(size(x))
logical :: check_posinf(size(x))
logical :: check_neginf(size(x))
logical :: success_SP
logical :: success_DP
PINF_SP = ieee_value(1.0, ieee_positive_inf)
NINF_SP = ieee_value(1.0, ieee_negative_inf)
QNAN_SP = ieee_value(1.0, ieee_quiet_nan)
SNAN_SP = ieee_value(1.0, ieee_signaling_nan)
PNORMAL_SP = ieee_value(1.0, ieee_positive_normal)
NNORMAL_SP = ieee_value(1.0, ieee_negative_normal)
PDENORMAL_SP = ieee_value(1.0, ieee_positive_denormal)
NDENORMAL_SP = ieee_value(1.0, ieee_negative_denormal)
PZERO_SP = ieee_value(1.0, ieee_positive_zero)
NZERO_SP = ieee_value(1.0, ieee_negative_zero)
PINF_DP = ieee_value(1.0D0, ieee_positive_inf)
NINF_DP = ieee_value(1.0D0, ieee_negative_inf)
QNAN_DP = ieee_value(1.0D0, ieee_quiet_nan)
SNAN_DP = ieee_value(1.0D0, ieee_signaling_nan)
PNORMAL_DP = ieee_value(1.0D0, ieee_positive_normal)
NNORMAL_DP = ieee_value(1.0D0, ieee_negative_normal)
PDENORMAL_DP = ieee_value(1.0D0, ieee_positive_denormal)
NDENORMAL_DP = ieee_value(1.0D0, ieee_negative_denormal)
PZERO_DP = ieee_value(1.0D0, ieee_positive_zero)
NZERO_DP = ieee_value(1.0D0, ieee_negative_zero)
x = [PINF_SP, NINF_SP, &
& QNAN_SP, SNAN_SP, &
& PNORMAL_SP, NNORMAL_SP, &
& PDENORMAL_SP, NDENORMAL_SP, &
& PZERO_SP, NZERO_SP, &
& huge(0.0), -huge(0.0), &
& epsilon(0.0), -epsilon(0.0), &
& tiny(0.0), -tiny(0.0), &
& 0.0, 1.0, -1.0]
y = [PINF_DP, NINF_DP, &
& QNAN_DP, SNAN_DP, &
& PNORMAL_DP, NNORMAL_DP, &
& PDENORMAL_DP, NDENORMAL_DP, &
& PZERO_DP, NZERO_DP, &
& huge(0.0D0), -huge(0.0D0), &
& epsilon(0.0D0), -epsilon(0.0D0), &
& tiny(0.0D0), -tiny(0.0D0), &
& 0.0D0, 1.0D0, -1.0D0]
check_finite = [.false., .false., &
& .false., .false., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., &
& .true., .true., .true.]
check_inf = [.true., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_posinf = [.true., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_neginf = [.false., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
check_nan = [.false., .false., &
& .true., .true., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., &
& .false., .false., .false.]
write (*, '(/(1A))') 'Single-precision tests:'
if (.not. all(is_finite(x) .eqv. check_finite)) then
write (*, *) 'True IS_FINITE: ', check_finite
write (*, *) 'IS_FINITE: ', is_finite(x)
write (*, *) 'IEEE_IS_FINITE: ', ieee_is_finite(x)
end if
if (.not. all(is_inf(x) .eqv. check_inf)) then
write (*, *) 'True IS_INF: ', check_inf
write (*, *) 'IS_INF: ', is_inf(x)
write (*, *) 'IEEE_IS_INF: ', ieee_is_inf(x)
end if
if (.not. all(is_posinf(x) .eqv. check_posinf)) then
write (*, *) 'True IS_POSINF: ', check_posinf
write (*, *) 'IS_POSINF: ', is_posinf(x)
write (*, *) 'IEEE_IS_POSINF: ', ieee_is_posinf(x)
end if
if (.not. all(is_neginf(x) .eqv. check_neginf)) then
write (*, *) 'True IS_NEGINF: ', check_neginf
write (*, *) 'IS_NEGINF: ', is_neginf(x)
write (*, *) 'IEEE_IS_NEGINF: ', ieee_is_neginf(x)
end if
if (.not. all(is_nan(x) .eqv. check_nan)) then
write (*, *) 'True IS_NAN: ', check_nan
write (*, *) 'IS_NAN: ', is_nan(x)
write (*, *) 'IEEE_IS_NAN: ', ieee_is_nan(x)
end if
success_SP = all(is_finite(x) .eqv. check_finite) .and. all(is_inf(x) .eqv. check_inf) .and. &
& all(is_posinf(x) .eqv. check_posinf) .and. all(is_neginf(x) .eqv. check_neginf) .and. &
& all(is_nan(x) .eqv. check_nan)
if (.not. success_SP) then
write (*, '(1A)') 'Some tests failed. The data is'
write (*, '(1P, 20E8.0)') x
else
write (*, '(1A)') 'All tests succeed.'
end if
write (*, '(/(1A))') 'Double-precision tests:'
if (.not. all(is_finite(y) .eqv. check_finite)) then
write (*, *) 'True IS_FINITE: ', check_finite
write (*, *) 'IS_FINITE: ', is_finite(y)
write (*, *) 'IEEE_IS_FINITE: ', ieee_is_finite(y)
end if
if (.not. all(is_inf(y) .eqv. check_inf)) then
write (*, *) 'True IS_INF: ', check_inf
write (*, *) 'IS_INF: ', is_inf(y)
write (*, *) 'IEEE_IS_INF: ', ieee_is_inf(y)
end if
if (.not. all(is_posinf(y) .eqv. check_posinf)) then
write (*, *) 'Truth IS_POSINF: ', check_posinf
write (*, *) 'IS_POSINF: ', is_posinf(y)
write (*, *) 'IEEE_IS_POSINF: ', ieee_is_posinf(y)
end if
if (.not. all(is_neginf(y) .eqv. check_neginf)) then
write (*, *) 'True IS_NEGINF: ', check_neginf
write (*, *) 'IS_NEGINF: ', is_neginf(y)
write (*, *) 'IEEE_IS_NEGINF: ', ieee_is_neginf(y)
end if
if (.not. all(is_nan(y) .eqv. check_nan)) then
write (*, *) 'True IS_NAN: ', check_nan
write (*, *) 'IS_NAN: ', is_nan(y)
write (*, *) 'IEEE_IS_NAN: ', ieee_is_nan(y)
end if
success_DP = all(is_finite(y) .eqv. check_finite) .and. all(is_inf(y) .eqv. check_inf) .and. &
& all(is_posinf(y) .eqv. check_posinf) .and. all(is_neginf(y) .eqv. check_neginf) .and. &
& all(is_nan(y) .eqv. check_nan)
if (.not. success_DP) then
write (*, '(1A)') 'Some tests failed. The data is'
write (*, '(1P, 20D8.0)') y
else
write (*, '(1A)') 'All tests succeed.'
end if
write (*, *) ''
success = success_SP .and. success_DP
end function test_infnan
end module test_infnan_mod
program test
use, non_intrinsic :: test_infnan_mod, only : test_infnan
if (.not. test_infnan()) then
stop 1
end if
end program test
2 Likes
Hi @hkvzjal , thank you for your effort! The result is very impressive. However, I am afraid that `real(2*wp)ˋ is not standard conforming.
1 Like
I think what you really want here is
use consts_mod, only: wp => sp
real(wp), intent(in) :: x
integer, parameter :: mp = selected_real_kind(precision(x)+1) ! a more precise kind
real(mp), parameter :: huge_2wp = real(huge(x), mp)*2
2 Likes
Thanks @everythingfunctional I have edited the post to include a modified version of the code taking into account your suggestion!! At first I was going for something like
real(mp), parameter :: huge_2wp = real(huge(x), mp) + some_delta
Simply to avoid tolerance problems with large numbers. But indeed huge_2wp = real(huge(x), mp)*2 does the work even better 
My intention was just to show that by using a constant value storing the huge of the working_precission in a variable with twice wp, it was possible to circumvent the tolerance issues. These huge values can be stored as constants in the module, for instance by initializing
module consts_mod
use iso_fortran_env ,only: real32, real64, real128
implicit none
private
public :: SP, DP, QP, huge_2sp, huge_2dp
integer, parameter :: SP = real32
integer, parameter :: DP = real64
integer, parameter :: QP = real128
real(DP), parameter :: huge_2sp = real(huge(0._SP), DP)*2
real(QP), parameter :: huge_2dp = real(huge(0._DP), QP)*2
end module consts_mod
And then use huge_2sp or huge_2dp accordingly.
1 Like
Thank eveyone for the input.
I have finally found an implementation that passes all my tests. See
The key is to wrap the intrinsic huge into a function huge_value, which must be implemented in a file separated from is_inf, is_finite etc. This somehow fools gfortran-13 so that the compiler does not smartly optimize abs(x) <= huge(x) into true even if x is not finite.
N.B.: Indeed, my implementation supports the quadruple precision as well — even though I know it is rarely used in practice. The abovementioned repo does not handle the quadruple case, but the version included in PRIMA does so.
Did you try to combine -Ofast with -flto? Inlining the function might get you back to square zero … I’m just guessing here. I saw that the intrinsic isnan from gnu is bugged in gfortran 13 with -Ofast so there might be an intrinsic issue that is “just” bypassed by reducing optimization.
1 Like
Hi @hkvzjal , you are right : The functions do not work with `gfortran-13 -Ofast -flto` · Issue #29 · equipez/infnan · GitHub . Thank you for pointing this out.
I saw that the intrinsic isnan from gnu is bugged in gfortran 13 with -Ofast
However, I am not sure whether this is a bug. I doubt that gfortran intends to make isnan work properly when it is invoked with -Ofast.
1 Like
Is this because NaNs are not generated with this option, or because the tests for them are optimized away even if they are generated? I think if a compiler knows that a NaN is not generated, then it is allowed subsequently to optimize away the test for it.
I guess gfortran assumes that everything is finite when it is invoked with -Ofast, which is expected. So the behavior does not seem to be a bug.
Optimize Options (Using the GNU Compiler Collection (GCC))
Optimize Options (Using the GNU Compiler Collection (GCC))
-Ofast
Disregard strict standards compliance. -Ofast enables all -O3 optimizations. It also enables optimizations that are not valid for all standard-compliant programs. It turns on -ffast-math, -fallow-store-data-races and the Fortran-specific -fstack-arrays, unless -fmax-stack-var-size is specified, and -fno-protect-parens. It turns off -fsemantic-interposition.
-ffast-math
Sets the options -fno-math-errno, -funsafe-math-optimizations, -ffinite-math-only, -fno-rounding-math, -fno-signaling-nans, -fcx-limited-range and -fexcess-precision=fast.
It is very impressive that even the following version fails with gfortran-13 -Ofast -flto :
As you can see, this version uses posinf(x) and neginf(x) instead of huge_value(x) and -huge_value(x) when implementing is_inf, is_finite, etc. I thought that the definition of posinf and neginf should be complicated enough to avoid them being optimized out, but gfortran-13 turned out much more capable in optimizing the code than I expected.
1 Like