# Challenge: Testing Inf and NaN with \`gfortran-13 -Ofast\`

**URL:** <https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481>\
**Category:** Uncategorized\
**Created:** [September 7, 2023, 6:16am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481 "2023-09-07T06:16:06Z")\
**Posts on this page:** 19\
**Page:** 2

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [September 9, 2023, 2:57pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/21 "2023-09-09T14:57:22Z")

</div>

One option would be to write the routines directly in assembly. This way the compiler is removed out of the picture.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [September 9, 2023, 3:36pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/22 "2023-09-09T15:36:10Z")

</div>

@zaikunzhang out of curiosity I tried to play with your `is_finite_sp` implementation and got the following:

```auto
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:

```auto
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](https://godbolt.org/z/fTde37hxT) 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 🤔

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [September 9, 2023, 6:12pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/23 "2023-09-09T18:12:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![urbanjost](https://avatars.discourse-cdn.com/v4/letter/u/0ea827/32.png) [@urbanjost](https://fortran-lang.discourse.group/u/urbanjost)\
**Post date:** [September 10, 2023, 1:28pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/24 "2023-09-10T13:28:13Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [September 11, 2023, 12:15am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/25 "2023-09-11T00:15:53Z")

</div>

Here is a possible solution:

> <https://github.com/equipez/infnan/pull/24>
>
> Replace ad-hoc comparisons with comparison to IEEE defined values.
> 
> Additional…ly, the Makefile had dos line endings for some reason...

However, [it contains non-standard Fortran code](https://github.com/equipez/infnan/pull/24#issuecomment-1711188808) (`sizeof(x)` and `integer(sizeof(x))`).

---

<div class="post-metadata">

**Author:** ![urbanjost](https://avatars.discourse-cdn.com/v4/letter/u/0ea827/32.png) [@urbanjost](https://fortran-lang.discourse.group/u/urbanjost)\
**Post date:** [September 11, 2023, 12:56am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/26 "2023-09-11T00:56:28Z")

</div>

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**
>
> ```fortran
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [September 11, 2023, 1:03am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/27 "2023-09-11T01:03:21Z")

</div>

Hi @urbanjost ,

> [@urbanjost](#):
>
> 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

Thank you for raising this point. I have also observed similar problems during the development of [PRIMA](https://github.com/libprima/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.

> [@urbanjost](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![gareth](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@gareth](https://fortran-lang.discourse.group/u/gareth)\
**Post date:** [September 11, 2023, 2:41am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/28 "2023-09-11T02:41:25Z")

</div>

Maybe the non-standard `sizeof` can be converted to F2008 standard code using [c\_sizeof](https://gcc.gnu.org/onlinedocs/gfortran/C_005fSIZEOF.html) or [storage\_size](https://gcc.gnu.org/onlinedocs/gfortran/STORAGE_005fSIZE.html)? Note the latter gives the size in bits, the former in bytes.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [September 11, 2023, 5:50am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/29 "2023-09-11T05:50:22Z")

</div>

Maybe. However, it is still nonstandard to relate the storage size with the kind, i.e., `integer(storage_size(x)) :: y`.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [September 12, 2023, 9:10pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/30 "2023-09-12T21:10:50Z")

</div>

@zaikunzhang I managed to get all your tests to pass with gfortran 13 -Ofast, you can check and run the modifications here: [Compiler Explorer](https://godbolt.org/z/d1YxPGaoK)

 ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/f/f88444ace610dbc3f5734c32d01da4be030e8826.png)

> **Full code**
>
> ```auto
> 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**
>
> ```auto
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [September 13, 2023, 12:18am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/31 "2023-09-13T00:18:30Z")

</div>

Hi @hkvzjal , thank you for your effort! The result is very impressive. However, I am afraid that `real(2\*wp)ˋ is not standard conforming.

---

<div class="post-metadata">

**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\
**Post date:** [September 13, 2023, 12:40am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/32 "2023-09-13T00:40:37Z")

</div>

> [@hkvzjal](#):
>
> ```auto
> use consts_mod, only : wp=>SP
> real(2*wp), parameter :: huge_2wp = 3.40282347D+38
> real(wp), intent(in) :: x
> 
> ```

I think what you really want here is

```auto
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

```

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [September 13, 2023, 5:23am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/33 "2023-09-13T05:23:08Z")

</div>

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

```auto
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 🙂

> [@zaikunzhang](#):
>
> I am afraid that `real(2\*wp)ˋ is not standard conforming.

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

```auto
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.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [October 2, 2023, 4:59pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/34 "2023-10-02T16:59:01Z")

</div>

Thank eveyone for the input.

I have finally found an implementation that passes all my tests. See

> **[GitHub - equipez/infnan: Fortran functions for checking Inf/NaN. They should...](https://github.com/equipez/infnan)**
>
> Fortran functions for checking Inf/NaN. They should work even if the compilers are invoked with aggressive optimization flags such as \`gfortran -Ofast\`. - GitHub - equipez/infnan: Fortran function...

The key is to wrap the intrinsic `huge` into a function [`huge_value`](https://github.com/equipez/infnan/blob/main/huge_value.f90), 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](https://github.com/libprima/prima/blob/main/fortran/common/infnan.F90) does so.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 2, 2023, 5:14pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/35 "2023-10-02T17:14:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [October 2, 2023, 5:34pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/36 "2023-10-02T17:34:53Z")

</div>

Hi @hkvzjal , you are right : [The functions do not work with `gfortran-13 -Ofast -flto` · Issue #29 · equipez/infnan · GitHub](https://github.com/equipez/infnan/issues/29) . 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`.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [October 2, 2023, 8:20pm UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/37 "2023-10-02T20:20:46Z")

</div>

> [@zaikunzhang](#):
>
> I doubt that `gfortran` intends to make `isnan` work properly when it is invoked with `-Ofast`.

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.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [October 3, 2023, 4:48am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/39 "2023-10-03T04:48:29Z")

</div>

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.

[gcc.gnu.org](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)

### [Optimize Options (Using the GNU Compiler Collection (GCC))](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)

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.

---

<div class="post-metadata">

**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\
**Post date:** [October 3, 2023, 4:52am UTC](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481/40 "2023-10-03T04:52:55Z")

</div>

It is very impressive that even the following version fails with `gfortran-13 -Ofast -flto` :

> **[GitHub - equipez/infnan at d9504b0effa658fba03e4bbf772354404b03db55](https://github.com/equipez/infnan/tree/d9504b0effa658fba03e4bbf772354404b03db55)**
>
> d9504b0effa658fba03e4bbf772354404b03db55

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.

[Previous page](https://fortran-lang.discourse.group/t/challenge-testing-inf-and-nan-with-gfortran-13-ofast/6481.md?page=1)