It is not in GFortran.
It is NaN in gfortran if a variable x is set to 0.0 and you evaluate x/x, but there is a compile-time error if x is a named constant or if you evaluate 0.0/0.0.
implicit none
real :: x
x = 0.0
print*,x/x ! NaN
! print*,0.0/0.0
! Line above with gfortran gives 'Error: Division by zero at (1)'
end
How to set a named constant to NaN is discussed at floating point - Having parameter (constant) variable with NaN value in Fortran - Stack Overflow.
1 Like
The second case is detected by the compiler, which apparently is concerned only by the value of the denominator.
The first NaN is occurring inside the FPU, according to IEEE754 and independently of the compiler.
Maybe it would more coherent for GFortran to consider 0.0/0.0 as a NaN, and to accept to compile the code, with a warning.
Note that ifx is compiling the code print *,0.0/0.0 without even a warning.
@vmagnin said āNote that ifx is compiling the code print *,0.0/0.0 without even a warning.ā That depends on oneās choice of options. When compiling a new program I always use
-assume protect_parens -check all,nouninit -fmath-errno -g -ftrapuv -traceback -warn interface -standard-semantics
and that gave me a compile-time error and told me where to look in the program. Other ifx options will stop it objecting. When using a Fortran 95 compiler with some f2003 features I read the character string āNaNā into a real variable when I want a NaN. Nowadays one can still do that, or use ieee_arithmetic, which offers both a signaling NaN and a quiet one. Neither of them need be the same as the read-in NaN.
1 Like
This feels really inconsistent to me. Is this a bug?
In R and C++, 0.0/0.0 is NAN. In Python, it is an error. Not sure how this behavior should be.
There are many situations where the compile time evaluation of an expression is different from the run time evaluation of the equivalent expression. This is kind of an oddball corner case that behaves that way. After all, from the compilerās perspective, why would anyone intentionally put the constant expression 0.0/0.0 in their source code?
If I need to create a NAN but donāt want to use the ieee_arithmetic, I can use 0.0/0.0.
1 Like
In the wikipedia page, 0.0/0.0 is cited several times as an example of NaN:
Do we have GFortran developers here? I believe this should be considered a bug.
If you want to have a NaN parameter for a conditional evaluation without using ieee_arithmetic module you could define your NaN parameters as follows:
real(4), parameter :: NaN_sp = real(z'7fc00000',4)
real(8), parameter :: NaN_dp = real(z'7ff8000000000000',8)
real(16), parameter :: NaN_qp = real(z'7fff8000000000000000000000000000',16)
print *, NaN_sp
print *, NaN_dp
print *, NaN_qp
end
Yet, I do agree that it might be a bug or at least a missing feature from gfortran. ifort/ifx will print a remark if you do define the parameter real(4), parameter :: NaN_sp = 0.0/0.0 but it will effectively affect NaN to the parameter.
remark #7954: The floating invalid condition was detected while evaluating this operation; the result is a NaN.
EDIT
this code snippet is not meant to be portable across all compilers/architectures, it is a quick idea for the user asking for ideas regarding the specific compiler targeted at the beginning of the discussion.
@hkvzjalās NaN example is, I suspect, not portable. I donāt now have access to NAG but when I did I had to use different kind numbers from the 4, 8, 16 in the program. And two compilers I do have access to gave unexpected output. AOCC flang:
NaN
NaN
5.9722037830479493269174770714937169E-4947
G95 warned at compile time about integer overflow of real(ā¦ ,8) and real(ā¦ ,16) but not real(ā¦ ,4) and printed
2.1432893E+9
9.221120237041091E+18
1.701385873120399643178730384677195E+38
Sadly, G95 is no longer maintained but itās still useful for f95 and many f2003 features.
It is not portable for (at least) three reasons. First, is the fact that it has hard-coded integer literals as KIND values. It is trivial to change that part of the code to use parameters in a portable way. In fact, I donāt know why anyone posting code here would not post code with parameters; it is an often-discussed topic, and it is easy to do. The next reason involves the conversion of the integer strings to real values using REAL() (which I think is equivalent to the TRANSFER() convention in this context). That is not portable due to the underlying byte-ordering issues. The posted code assumes little-endian conventions with consecutive bytes constituting, low to high, the mantissa, the exponent, and the sign. As a counter example, we all know of one legacy floating point format, on a little-endian machine, that put the sign and exponent in the middle of the floating point word. The third reason is that it assumes IEEE floating point format conventions. Although that is the most common floating point format these days, that does not make it portable to all fortran compilers.
Having said all that, with the state of the fortran language right now, that code (using parameters for the KINDs, of course) is probably the most portable code that one could write. It is portable to little-endian machines that support IEEE floating point. The obvious way to change the language to make it more portable would be to allow the relevant functions in the IEEE intrinsic modules to be used in initializations. Of course, that would only make the code portable on the subset of fortran compilers that support those modules and that floating point format. Universal support does not seem to be practical, or even possible, beyond that.
I have written testnans.f90 below to be fairly portable, and to test what one gets from reading āNaNā, calculating 0.0/0.0, and the two IEEE NaNs, in each of the 3 or 4 real kinds available in these 4 compilers: g95, gfortran, ifort, ifx, in my Ubuntu system. I wonāt quote the output as itās fairly long, but I found it instructive. The program:
program testnans
use ieee_arithmetic, only: ieee_value,ieee_quiet_nan,ieee_signaling_nan
implicit none
integer,parameter:: sp = kind(1.0),dp = kind(1d0), &
eq = selected_real_kind(precision(1.0_dp)+1),ep = merge(eq,dp,eq>0), &
qq = selected_real_kind(precision(1.0_ep)+1),qp = merge(qq,ep,qq>0)
character:: cnan*3 = 'NaN', &
which(4)*9=[' read-in',' 0/0',' quiet','signaling']
character(*),parameter:: fmt = "(A,F4.0,I3,A,99(Z9.8:))"
real(sp):: spnan(4), sp0 = 0
real(dp):: dpnan(4), dp0 = 0
real(ep):: epnan(4), ep0 = 0
real(qp):: qpnan(4), qp0 = 0
integer :: iw,n(1) = 0, bits(4), bytes(4)
bits = [nint(log(real(maxexponent(sp0)))/log(2.0))+digits(sp0)+1 &
, nint(log(real(maxexponent(dp0)))/log(2.0))+digits(dp0)+1 &
, nint(log(real(maxexponent(ep0)))/log(2.0))+digits(ep0)+1 &
, nint(log(real(maxexponent(qp0)))/log(2.0))+digits(qp0)+1 ]
bytes = (bits(:)+1)/8
read(cnan,*) spnan(1)
spnan(2) = sp0/sp0
spnan(3)=ieee_value(sp0,ieee_quiet_nan)
spnan(4)=ieee_value(sp0,ieee_signaling_nan)
read(cnan,*) dpnan(1)
dpnan(2) = dp0/dp0
dpnan(3)=ieee_value(dp0,ieee_quiet_nan)
dpnan(4)=ieee_value(dp0,ieee_signaling_nan)
read(cnan,*) epnan(1)
epnan(2) = ep0/ep0
epnan(3)=ieee_value(ep0,ieee_quiet_nan)
epnan(4)=ieee_value(ep0,ieee_signaling_nan)
read(cnan,*) qpnan(1)
qpnan(2) = qp0/qp0
qpnan(3)=ieee_value(qp0,ieee_quiet_nan)
qpnan(4)=ieee_value(qp0,ieee_signaling_nan)
do iw=1,4
print fmt,which(iw),spnan(iw),bytes(1),' bytes',transfer(spnan(iw),n)
end do
print *
do iw=1,4
print fmt,which(iw),dpnan(iw),bytes(2),' bytes',transfer(dpnan(iw),n)
end do
print *
do iw=1,4
print fmt,which(iw),epnan(iw),bytes(3),' bytes',transfer(epnan(iw),n)
end do
if(qp==ep) stop ! as ifort, ifx have only 3 real kinds but gfortran, g95 have 4.
print *
do iw=1,4
print fmt,which(iw),qpnan(iw),bytes(4),' bytes',transfer(qpnan(iw),n)
end do
end program testnans
1 Like
I agree with this, especially if the programs are to be shared with others. But I think we also need to respect that āwe are all adults here,ā and everyone is free to make their own choices.
It is also hard to sell this point when even some vendors perpetuate the use of the legacy REAL* syntax (.e.g. Nvidia, and Intel (which at leasts uses green to denote this as an extension)), and itās a staple in various exemplary materials. (On the other hand, not documenting it would also be bad.)
As the name of the OP indicates, they use Fortran from R and primarily with gfortran if Iām not mistaken. For a long time, all R platforms were using gfortran, and when you search for tutorials about using Fortran from R, they generally recommend using double precision, or real(8), without any further consideration of why. Iāve just noticed now that the Writing R Extensions manual, has extended the section on portable Fortran code to discuss this issue.
In the past, Iāve suggested using a module like this one for new development:
module rkinds
use, intrinsic :: iso_c_binding
private
integer, parameter, public :: rp = c_double
end module
Ideally, this would be distributed together with the R platform, and documented with a table:
| Value | R Vector | Fortran |
|---|---|---|
| Real | numeric |
real(rp) |
| Complex | complex |
complex(rp) |
But I guess that the use of _rp on literals is perceived as unappetizing, especially if compared to the C++ Rcpp package, since in C++ floating-point literals are double by default (and instead single precision literals need to have the f postfix).
Other types in R can be problematic to inter-operate with Fortran. For example the boolean type is implemented like this:
typedef enum { FALSE = 0, TRUE /*, MAYBE */ } Rboolean;
So in principle it shares the representation of logical that gfortran uses, but not the one of Intel Fortran (unless certain compiler flags are used). Even if the representation is the same, you need an additional wrapper layer to get around the type-casting:
module rkinds
use, intrinsic :: iso_c_binding, sexp => c_ptr
private
public :: sexp
integer, parameter, public :: rp = c_double
enum, bind(c)
enumerator :: rfalse = 0
enumerator :: rtrue = 1
end enum
integer, parameter, public :: rboolean = kind(rfalse)
public :: as_rboolean, rtrue, rfalse
contains
integer(rboolean) function as_rboolean(x)
logical, value :: x
if (x) then
as_rboolean = rtrue
else
as_rboolean = rfalse
end if
end function
end module
Returning to the original issue, and assuming you want to use NaN from R, you could simply use the constant, defined in include/R_ext/Arith.h:
/* implementation of these : ../../main/arithmetic.c */
LibExtern double R_NaN; /* IEEE NaN */
LibExtern double R_PosInf; /* IEEE Inf */
LibExtern double R_NegInf; /* IEEE -Inf */
LibExtern double R_NaReal; /* NA_REAL: IEEE */
LibExtern int R_NaInt; /* NA_INTEGER:= INT_MIN currently */
which gets initialized by a function in src/main/arithmetic.c:
/* Arithmetic Initialization */
attribute_hidden void InitArithmetic(void)
{
R_NaInt = INT_MIN;
R_NaReal = R_ValueOfNA();
// we assume C99, so
#ifndef OLD
R_NaN = NAN;
R_PosInf = INFINITY;
R_NegInf = -INFINITY;
#else
R_NaN = 0.0/R_Zero_Hack;
R_PosInf = 1.0/R_Zero_Hack;
R_NegInf = -1.0/R_Zero_Hack;
#endif
}
I havenāt traced exactly where do the macros come from, but it should be enough to bind it like this:
module rparams
use, intrinsic :: iso_c_binding
private
implicit none
public :: r_nan, r_posinf, r_neginf
public :: r_isna, r_isnan, r_finite
real(c_double), bind(c,name="R_NaN") :: r_nan
real(c_double), bind(c,name="R_PosInf") :: r_posinf
real(c_double), bind(c,name="R_NegInf") :: r_neginf
protected :: r_nan, r_posinf, r_neginf
interface
function r_isna(x) bind(c,name='R_IsNA')
import c_double, c_int
real(c_double), value :: x
integer(c_int) :: r_isna
end function
function r_isnan(x) bind(c,name='R_IsNaN')
import c_double, c_int
real(c_double), value :: x
integer(c_int) :: r_isnan
end function
function r_finite(x) bind(c,name='R_finite')
import c_double, c_int
real(c_double), value :: x
integer(c_int) :: r_finite
end function
end interface
end module
Caveat: this will only work within an R session, where the hidden InitArithmetic function has been called (presumably when the R interpreter is launched). It is also tightly bound to the main R implementation.
(
there may be errors remaining, as I donāt have R running on the computer where I typed this answer.)
I was reading the AOCC documentation (never used it) and just found that they offer flang āas the default Fortran front-end compiler supporting F2008 and Real128 featuresāā¦ Not sure which flang do they refer to ( ? ) ā¦ I tried the flang in Compiler Explorer and I do get a proper NaN for the quad precision conversionā¦
For G95, Iāve never used it, I understand from here G95 - Wikipedia that gfortran was forked from it, as such I guess that it is safe to take it as the reference.
I agree, there is really no excuse for this. However, I do see a reason that compiler vendors might use hard coded literals in their other documentation (i.e. INTEGER(4), REAL(8), LOGICAL(2), COMPLEX(4), etc.). This allows them to kill two birds with one stone, they can document whatever feature they are describing, and at the same time document their KIND values for the various intrinsic types.
The downside is when new fortran users read this documentation, they might tend to follow the same conventions, which is not the best practice for programmers writing portable code.
1 Like
The previous program I posted here would not compile with AOCC flang because that compiler did not have ieee_support_nan for all its real kinds. This program is shorter because it only uses one real kind but it shows that different compilers give different NaNs. It compiled and ran with AOCC flang, gfortran, g95, ifort and ifx if suitable options were used. (For example avoid an option that would stop execution on attempting 0.0/0.0). The program:
program testnans2 ! F2003. Various single precision NaNs
use ieee_arithmetic, only: ieee_value,ieee_quiet_nan,ieee_signaling_nan,&
ieee_support_nan
implicit none
integer,parameter:: nans = 4
character:: cnan*3 = 'NaN', &
which(nans)*9 = [' read-in','zero/zero',' quiet','signaling']
character(*),parameter:: fmt = "(A,1X,F0.0,9(Z9.8:))"
real:: nan1(nans), zero = 0
integer :: iw, n(1) = 1
read(cnan,*) nan1(1)
nan1(2) = zero/zero
if(ieee_support_nan(zero)) then
nan1(3) = ieee_value(zero,ieee_quiet_nan)
nan1(4) = ieee_value(zero,ieee_signaling_nan)
else
nan1(3:4) = 0 ! some compilers want all array elements evaluated
end if
do iw=1,merge(nans,nans-2,ieee_support_nan(zero))
print fmt,which(iw),nan1(iw),transfer(nan1(iw),n)
end do
end program testnans2
Output from ifx or ifort: 3 different NaNs, read-in == zero/zero:
read-in NaN FFC00000
zero/zero NaN FFC00000
quiet NaN 7FC00000
signaling NaN 7FA00000
Output from gfortran or AOCC flang: 3 different NaNs, read-in == quiet
read-in NaN 7FC00000
zero/zero NaN FFC00000
quiet NaN 7FC00000
signaling NaN 7FA00000
Output from g95: 4 differrent NaNs:
read-in NaN 7F8EA46B
zero/zero NaN FFC00000
quiet NaN 7FFFFFFF
signaling NaN 7F800001