Should 0.0/0.0 be NaN?

Help
15

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.

(:warning: there may be errors remaining, as I don’t have R running on the computer where I typed this answer.)