IEEE_arithmetic module

I have been experimenting with the relatively new IEEE arithmetic module.

The following code runs correctly with the 32 bit compiler. A run time error occurs with the 64 bit compiler at the second print statement when using CheckMate

program test
use, intrinsic :: iso_fortran_env, only : real32
use, intrinsic :: ieee_arithmetic
implicit none
real(real32) :: x1
  print*, ieee_value(1.0_real32, ieee_quiet_nan)
  x1 = ieee_value(1.0_real32, ieee_quiet_nan)
  print*, x1
end program test

The same problem occurs for real64 x

Ken

Thank you for the feedback. I will take a look at this.

Paul,

I am not sure you will thank me for this feedback - I suspect it will be a lot of work to put this right.

program propagation
use, intrinsic :: iso_fortran_env, only : real32
use, intrinsic :: ieee_arithmetic
implicit none
real(real32) :: xd, yd
  xd = ieee_value(1.0_real32, ieee_quiet_nan)
  
  ! NaN **should** propagate through almost all standard Fortran 
  ! intrinsic functions.
  
  ! These functions show the expected behaviour
  
  yd = atan(xd)        ; print*, 'atan       ', yd
  yd = cos(xd)         ; print*, 'cos        ', yd
  yd = sin(xd)         ; print*, 'sin        ', yd
  yd = tan(xd)         ; print*, 'tan        ', yd
  
  ! These are examples of intrinsic functions that do not propagate NaN
  ! and generate runtime errors
  
!  yd = atanh(xd)       ; print*, 'atanh      ', yd
!  yd = acos(xd)        ; print*, 'acos       ', yd
!  yd = acosh(xd)       ; print*, 'acosh      ', yd
!  yd = asin(xd)        ; print*, 'asin       ', yd
!  yd = asinh(xd)       ; print*, 'asinh      ', yd
!  yd = cosh(xd)        ; print*, 'cosh       ', yd 
!  yd = erf(xd)         ; print*, 'erf        ', yd 
!  yd = exp(xd)         ; print*, 'exp        ', yd

end program propagation 

Ken

Thank you anyway. Taking a positive view, we still don't claim to be more that Fortran 95 compliant. Anything else is an extension.

No worries. I had a bit of fun this afternoon with this work around for sqrt(NaN).

module safe_sqrt_mod
  use, intrinsic :: iso_fortran_env, only: real32, real64
  use, intrinsic :: ieee_arithmetic
  implicit none
  private
  public :: safe_sqrt   ! Cannot name this SQRT i.e. over ride the INTRINSIC since
  interface safe_sqrt   ! the functions below end up being called resursively.
    module procedure safe_sqrtS, safe_sqrtD
    module procedure safe_sqrtC, safe_sqrtZ
  end interface safe_sqrt
contains
  !---------------- Real, single precision ----------------
  elemental pure function safe_sqrtS(x) result(y)
    real(real32), intent(in) :: x
    real(real32) :: y
    if (ieee_is_nan(x)) then
      y = x
    else
      y = sqrt(x)
    end if
  end function safe_sqrtS
  !---------------- Real, double precision ----------------
  elemental pure function safe_sqrtD(x) result(y)
    real(real64), intent(in) :: x
    real(real64) :: y
    if (ieee_is_nan(x)) then
      y = x
    else
      y = sqrt(x)
    end if
  end function safe_sqrtD
  !---------------- Complex, single precision ----------------
  elemental pure function safe_sqrtC(z) result(w)
    complex(real32), intent(in) :: z
    complex(real32) :: w
    real(real32) :: zr, zi
    zr  = real(z)
    zi  = aimag(z)
    if (ieee_is_nan(zr) .or. ieee_is_nan(zi)) then
      w = z
    else
      w = sqrt(z)
    end if
  end function safe_sqrtC
  !---------------- Complex, double precision ----------------
  elemental pure function safe_sqrtZ(z) result(w)
    complex(real64), intent(in) :: z
    complex(real64) :: w
    real(real64) :: zr, zi
    zr   = real(z)
    zi   = aimag(z)
    if (ieee_is_nan(zr) .or. ieee_is_nan(zi)) then
      w = z
    else
      w = sqrt(z)
    end if
  end function safe_sqrtZ
end module safe_sqrt_mod

program p
use, intrinsic :: iso_fortran_env, only : real64, real32
use, intrinsic :: ieee_arithmetic
use safe_sqrt_mod
real(real32) :: r1, r2
real(real64) :: r3, r4
complex(real32) :: z1, z2
complex(real64) :: z3, z4
r1 = 1.0 ; r2 = safe_sqrt(r1) ; print*, r1,r2
r1 = ieee_value(1.0_real32, ieee_quiet_nan) ; r2 = safe_sqrt(r1) ; print*, r1,r2
r3 = 1.0 ; r4 = safe_sqrt(r3) ; print*, r3,r4
r3 = ieee_value(1.0_real64, ieee_quiet_nan) ; r4 = safe_sqrt(r3) ; print*, r3,r4
z1 = ieee_value(1.0_real32, ieee_quiet_nan) ; z2 = safe_sqrt(z1) ; print*, z1,z2
z3 = ieee_value(1.0_real64, ieee_quiet_nan) ; z4 = safe_sqrt(z3) ; print*, z3,z4
z4 = sqrt(z3) ;  print*, z3,z4  ! This is the problem we are trying to avoid
end program p

The penny dropped 5 minutes after my last post.

See: https://www.dropbox.com/scl/fi/rd3l48ainibsf3hl4njfj/mysqrt.f95?rlkey=g59uqv3g2awry7x44teqohxo0&st=ls134ged&dl=0

This implementation uses two separate modules to safely override sqrt while still retaining access to the intrinsic. The first module (intrinsic_func) defines kind-specific wrapper procedures (sqrt2S, sqrt2D, sqrt2C, sqrt2Z) that call the true intrinsic sqrt directly. These wrappers exist solely to provide an unambiguous path to the intrinsic and to avoid recursive calls. The second module (safe_sqrt_mod) publicly re-exports the generic name sqrt and implements NaN-safe logic. Each safe_sqrt* routine first checks the argument for IEEE NaNs and, if present, returns the NaN unchanged; otherwise it delegates the computation to the corresponding wrapper in intrinsic_func. This two-module structure is necessary because Fortran does not provide a standard mechanism to explicitly reference an intrinsic once it has been overridden by a generic interface, and attempting to call sqrt directly would result in infinite recursion. The separation therefore ensures correctness, purity, elemental behaviour, and portability, where NaN propagation in intrinsics is unreliable.

Of course, you would need to do this for every intrinsic procedure where you want to propagate NaN.

My unpure routines, which currently STOP program execution can remain that way for the time being !

Nevertheless this was a very instructive exercise for me.

Thanks for sharing your exercises. Small code examples are priceless to learn anything.

Users have not to hesitate and freely post their exercises, they for sure are helpful for someone

For example Salford/Silverfrost set of examples they created 30 years ago is essentially all what anyone needs to learn Clearwin in less than a day

I thought it would be instructive to extend the previous code to handle cases where the real argument, or one of both parts of a complex argument were plus or minus infinity.

The first hurdle to overcome was that the standard does not require the IEEE_arithmetic module to have a logical functions for testing if a real is plus infinity or negative infinity. So there is a new module infinity_tests with functions is_pos_inf and is_neg_inf to cover the real and double precision real cases.

It was fairly easy to deal with the real case. However the complex case required much more thought. See the decision tables in the comments in the code.

Here is the link to the updated code:https://www.dropbox.com/scl/fi/5lur9mtovz8tvnspa6x7m/mysqrt.2f95.f95?rlkey=w5ofokp9o7go5o08fyvnzskw6&st=sxwtwii8&dl=0

With the 64 bit compiler:

1. in release mode the code returns the expected results
2. in debug mode the code returns the expected results
3. with checkmate the core fails at run time for the original issue identified at the top of this post.

With the 32 bit compler:

1. in release mode the code returns the expected results
2. In debug mode the code fails to compile - error 163 - Internal compiler error
3. With checkmate the code fails to compile - error 163 - Internal compiler error

Hopefully, some part of this post will be of use to somebody at sometime in the future!