Commit 40dc0c58 authored by Radford Neal's avatar Radford Neal

new scheme for ISNAN, etc.

parent 3e062406
New approach to ISNAN, ISNA, etc. Assumes IEEE floating point double
corresponds to 64-bit int in consistent endian way. Implements these
as static inline functions on that basis, using union. No longer any
need for ENABLE_ISNAN_TRICK, since this should work on all current
processors (any with consistent endianness).
......@@ -140,7 +140,6 @@ likely to be useful are as follows:
compilers are used with an Intel/AMD 64-bit processor, the
following flags are recommended:
CFLAGS='-g -O2 -mtune=native -DENABLE_ISNAN_TRICK'
CXXLAGS='-g -O2 -mtune=native'
FCLAGS='-g -O2 -mtune=native'
FFLAGS='-g -O2 -mtune=native'
......
......@@ -3658,13 +3658,6 @@ Tuning compilation to a specific @acronym{CPU} family (e.g.@:
@samp{-mtune=native} for @command{gcc}) can give worthwhile performance
gains, especially on older architectures such as @cputype{ix86}.
When configuring pqR, @code{-DENABLE_ISNAN_TRICK} can be included in
@code{CFLAGS} in order to enable use of a ``trick'' for testing
whether a real value is @code{NA} or @code{NaN} more quickly. This
trick should work for @cputype{ix86} systems, and perhaps for others
as well. A run-time check is made for whether it does appear to work,
and an error message is printed if not.
@node Making manuals, , Compilation flags, Configuration variables
@subsection Making manuals
......
......@@ -1662,20 +1662,25 @@ void SET_ATTRIB_TO_ANYTHING (SEXP, SEXP);
#undef NA_LOGICAL
#define NA_LOGICAL INT_MIN
/* Redefine ISNAN to make use of a trick, if ENABLE_ISNAN_TRICK is
defined (with a -D argument in CFLAGS). This is done here only, not
in Arith.h, because the macro implementing the trick evaluates its
argument twice, which is bad if it has side effects. Such macro
calls are avoided in the interpreter, but may occur in packages.
The trick is faster for many non-NaN and non-NA numbers. It relies
on the results of converting NaN, NA, -NaN, and -NA to int all being
the same, which is checked for in InitArithmetic. */
#ifdef ENABLE_ISNAN_TRICK
# undef ISNAN
# define ISNAN(x) ((int)(x) == R_NaN_cast_to_int && isnan(x) != 0)
#endif
/* Define R_INFINITE and ISNAN_NOT_NA here as inline functions, using a
union with uint64_t. ISNAN, ISNA, and R_FINITE are defined similarly
in Arith.h. */
static inline int R_INFINITE (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) == ((uint64_t)0x7ff << 53);
}
static inline int ISNAN_NOT_NA (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) > ((uint64_t)0x7ff << 53) /* a NaN, but... */
&& (un.u & (((uint64_t)1<<32)-1)) != 1954; /* not NA */
}
#endif /* DEFN_H_ */
/*
......
/*
* pqR : A pretty quick version of R
* Copyright (C) 2013, 2014 by Radford M. Neal
* Copyright (C) 2013, 2014, 2015 by Radford M. Neal
*
* Based on R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
......@@ -29,7 +29,7 @@
/* Only for use where config.h has not already been included */
#if defined(HAVE_GLIBC2) && !defined(_BSD_SOURCE)
/* ensure that finite and isnan are declared */
/* ensure that finite and isnan are declared - probably now obsolete */
# define _BSD_SOURCE 1
#endif
......@@ -37,7 +37,7 @@
#ifdef __cplusplus
extern "C" {
#elif !defined(NO_C_HEADERS)
/* needed for isnan and isfinite, neither of which are used under C++ */
/* may now be unnecessary.. */
# include <math.h>
#endif
......@@ -72,32 +72,44 @@ LibExtern int R_NaInt; /* NA_INTEGER:= INT_MIN currently */
/* NA_STRING is a SEXP, so defined in Rinternals.h */
int R_IsNA(double); /* True for R's NA only */
int R_IsNaN(double); /* True for special NaN, *not* for NA */
int R_IsNaN(double); /* True for any NaN that is *not* NA */
int R_finite(double); /* True if none of NA, NaN, +/-Inf */
#define ISNA(x) R_IsNA(x)
/* ISNAN(): True for *both* NA and NaN.
NOTE: some systems do not return 1 for TRUE.
Also note that C++ math headers specifically undefine
isnan if it is a macro (it is on OS X and in C99),
hence the workaround. This code also appears in Rmath.h
In pqR, this definition may be changed in Defn.h when ENABLE_ISNAN_TRICK
is defined, to give a faster version in the interpreter. */
/* The code below also appears in Rmath.h */
#ifdef __cplusplus
int R_isnancpp(double); /* in arithmetic.c */
# define ISNAN(x) (R_isnancpp(x))
#else
# define ISNAN(x) (isnan(x) != 0)
#endif
/* The following is only defined inside R */
#ifdef HAVE_WORKING_ISFINITE
/* isfinite is defined in <math.h> according to C99 */
# define R_FINITE(x) isfinite(x)
int R_isnancpp(double); /* in arithmetic.c */
#define ISNAN(x) (R_isnancpp(x))
#define ISNA(x) (R_IsNA(x))
#define R_FINITE(x) (R_finite(x))
#else
# define R_FINITE(x) R_finite(x)
#include <stdint.h>
static inline int ISNAN (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) > ((uint64_t)0x7ff << 53);
}
static inline int ISNA (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return ((un.u >> 52) & 0x7ff) == 0x7ff
&& (un.u & (((uint64_t)1<<32)-1)) == 1954;
}
static inline int R_FINITE (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) < ((uint64_t)0x7ff << 53);
}
#endif
#ifdef __cplusplus
......
......@@ -617,20 +617,44 @@ double logspace_sub(double logx, double logy);
#if defined(MATHLIB_STANDALONE) && !defined(MATHLIB_PRIVATE_H)
/* second is defined by nmath.h */
/* If isnan is a macro, as C99 specifies, the C++
math header will undefine it. This happens on OS X.
The optimization of ISNAN done in Arith.h is not done here,
at least for now. */
# ifdef __cplusplus
int R_isnancpp(double); /* in mlutils.c */
# define ISNAN(x) R_isnancpp(x)
# else
# define ISNAN(x) (isnan(x)!=0)
# endif
/* The following definitions of ISNAN, ISNA, and R_FINITE are the same as
in Arith.h */
#ifdef __cplusplus
int R_isnancpp(double), R_IsNA(double), R_finite(double); /* in arithmetic.c */
#define ISNAN(x) (R_isnancpp(x))
#define ISNA(x) (R_IsNA(x))
#define R_FINITE(x) (R_finite(x))
#else
#include <stdint.h>
static inline int ISNAN (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) > ((uint64_t)0x7ff << 53);
}
static inline int ISNA (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return ((un.u >> 52) & 0x7ff) == 0x7ff
&& (un.u & (((uint64_t)1<<32)-1)) == 1954;
}
static inline int R_FINITE (double x)
{
union { double d; uint64_t u; } un;
un.d = x;
return (un.u << 1) < ((uint64_t)0x7ff << 53);
}
#endif
# define R_FINITE(x) R_finite(x)
int R_finite(double);
# ifdef WIN32 /* not Win32 as no config information */
# ifdef RMATH_DLL
......
......@@ -109,28 +109,6 @@ int matherr(struct exception *exc)
/* Hack to avoid possibly-incorrect constant folding. */
volatile double R_Zero_Hack = 0.0;
/* gcc had problems with static const on AIX and Solaris
Solaris was for gcc 3.1 and 3.2 under -O2 32-bit on 64-bit kernel */
#ifdef _AIX
#define CONST
#elif defined(sparc) && defined (__GNUC__) && __GNUC__ == 3
#define CONST
#else
#define CONST const
#endif
#ifdef WORDS_BIGENDIAN
static CONST int hw = 0;
static CONST int lw = 1;
#else /* !WORDS_BIGENDIAN */
static CONST int hw = 1;
static CONST int lw = 0;
#endif /* WORDS_BIGENDIAN */
typedef union
{ double value;
unsigned int word[2];
} R_ieee_double;
static double R_ValueOfNA(void)
{
......@@ -139,42 +117,24 @@ static double R_ValueOfNA(void)
int R_IsNA(double x)
{
if (isnan(x)) {
R_ieee_double y;
y.value = x;
return (y.word[lw] == 1954);
}
return 0;
return ISNA(x);
}
int R_IsNaN(double x)
{
if (isnan(x)) {
R_ieee_double y;
y.value = x;
return (y.word[lw] != 1954);
}
return 0;
return ISNAN_NOT_NA(x);
}
/* ISNAN uses isnan, which is undefined by C++ headers
This workaround is called only when ISNAN() is used
in a user code in a file with __cplusplus defined */
int R_isnancpp(double x)
int R_finite(double x)
{
return (isnan(x)!=0);
return R_FINITE(x);
}
/* Used to define ISNAN for C++ in Arith.h */
/* Mainly for use in packages */
int R_finite(double x)
int R_isnancpp(double x)
{
#ifdef HAVE_WORKING_ISFINITE
return isfinite(x);
#else
return (!isnan(x) & (x != R_PosInf) & (x != R_NegInf));
#endif
return ISNAN(x);
}
......@@ -188,13 +148,6 @@ void attribute_hidden InitArithmetic()
R_PosInf = 1.0/R_Zero_Hack;
R_NegInf = -1.0/R_Zero_Hack;
R_NaN_cast_to_int = (int) R_NaN;
#ifdef ENABLE_ISNAN_TRICK
if (R_NaN_cast_to_int != (int) R_NaReal
|| R_NaN_cast_to_int != (int) (-R_NaReal)
|| R_NaN_cast_to_int != (int) (-R_NaN))
error("Integer casts of NaN, NA, -NaN, -NA differ, don't define ENABLE_ISNAN_TRICK");
#endif
}
/* Keep these two in step */
......
......@@ -1961,23 +1961,8 @@ static SEXP do_fast_isna (SEXP call, SEXP op, SEXP x, SEXP rho, int variant)
ret = TRUE; goto vret;
}
if (VARIANT_KIND(variant) == VARIANT_OR) {
int i = 0;
if (n&1) {
for (i = 0; i < n; i++)
if (ISNAN(REAL(x)[i])) { ret = TRUE; goto vret; }
i += 1;
}
if (n&2) {
if (ISNAN(REAL(x)[i]+REAL(x)[i+1])
&& (ISNAN(REAL(x)[i])
|| ISNAN(REAL(x)[i+1]))) { ret = TRUE; goto vret; }
i += 2;
}
for ( ; i < n; i += 4)
if (ISNAN(REAL(x)[i]+REAL(x)[i+1]+REAL(x)[i+2]+REAL(x)[i+3])
&& (ISNAN(REAL(x)[i])
|| ISNAN(REAL(x)[i+1])
|| ISNAN(REAL(x)[i+2])
|| ISNAN(REAL(x)[i+3]))) { ret = TRUE; goto vret; }
ret = FALSE; goto vret;
}
for (i = 0; i < n; i++)
......@@ -2134,12 +2119,12 @@ static SEXP do_fast_isnan (SEXP call, SEXP op, SEXP x, SEXP rho, int variant)
case REALSXP:
if (VARIANT_KIND(variant) == VARIANT_AND) {
for (i = 0; i < n; i++)
if (!R_IsNaN(REAL(x)[i])) { ret = FALSE; goto vret; }
if (!ISNAN_NOT_NA(REAL(x)[i])) { ret = FALSE; goto vret; }
ret = TRUE; goto vret;
}
if (VARIANT_KIND(variant) == VARIANT_OR) {
for (i = 0; i < n; i++)
if (R_IsNaN(REAL(x)[i])) { ret = TRUE; goto vret; }
if (ISNAN_NOT_NA(REAL(x)[i])) { ret = TRUE; goto vret; }
ret = FALSE; goto vret;
}
for (i = 0; i < n; i++)
......@@ -2148,19 +2133,19 @@ static SEXP do_fast_isnan (SEXP call, SEXP op, SEXP x, SEXP rho, int variant)
case CPLXSXP:
if (VARIANT_KIND(variant) == VARIANT_AND) {
for (i = 0; i < n; i++)
if (! (R_IsNaN(COMPLEX(x)[i].r)
|| R_IsNaN(COMPLEX(x)[i].i))) { ret = FALSE; goto vret; }
if (! (ISNAN_NOT_NA(COMPLEX(x)[i].r)
|| ISNAN_NOT_NA(COMPLEX(x)[i].i))) { ret = FALSE; goto vret; }
ret = TRUE; goto vret;
}
if (VARIANT_KIND(variant) == VARIANT_OR) {
for (i = 0; i < n; i++)
if (R_IsNaN(COMPLEX(x)[i].r)
|| R_IsNaN(COMPLEX(x)[i].i)) { ret = TRUE; goto vret; }
if (ISNAN_NOT_NA(COMPLEX(x)[i].r)
|| ISNAN_NOT_NA(COMPLEX(x)[i].i)) { ret = TRUE; goto vret; }
ret = FALSE; goto vret;
}
for (i = 0; i < n; i++)
LOGICAL(ans)[i] = (R_IsNaN(COMPLEX(x)[i].r) ||
R_IsNaN(COMPLEX(x)[i].i));
LOGICAL(ans)[i] = (ISNAN_NOT_NA(COMPLEX(x)[i].r) ||
ISNAN_NOT_NA(COMPLEX(x)[i].i));
break;
default:
errorcall(call,
......@@ -2353,26 +2338,26 @@ static SEXP do_fast_isinfinite (SEXP call, SEXP op, SEXP x, SEXP rho,
case REALSXP:
if (VARIANT_KIND(variant) == VARIANT_AND) {
for (i = 0; i < n; i++)
if (ISNAN(REAL(x)[i]) || R_FINITE(REAL(x)[i]))
if (!R_INFINITE(REAL(x)[i]))
{ ret = FALSE; goto vret; }
ret = TRUE; goto vret;
}
if (VARIANT_KIND(variant) == VARIANT_OR) {
for (i = 0; i < n; i++)
if (!ISNAN(REAL(x)[i]) && !R_FINITE(REAL(x)[i]))
if (R_INFINITE(REAL(x)[i]))
{ ret = TRUE; goto vret; }
ret = FALSE; goto vret;
}
for (i = 0; i < n; i++)
LOGICAL(ans)[i] = !ISNAN(REAL(x)[i]) && !R_FINITE(REAL(x)[i]);
LOGICAL(ans)[i] = R_INFINITE(REAL(x)[i]);
break;
case CPLXSXP:
if (VARIANT_KIND(variant) == VARIANT_AND) {
for (i = 0; i < n; i++) {
xr = COMPLEX(x)[i].r;
xi = COMPLEX(x)[i].i;
if ((ISNAN(xr) || R_FINITE(xr))
&& (ISNAN(xi) || R_FINITE(xi))) { ret = FALSE; goto vret; }
if (!R_INFINITE(xr)
&& !R_INFINITE(xi)) { ret = FALSE; goto vret; }
}
ret = TRUE; goto vret;
}
......@@ -2380,16 +2365,14 @@ static SEXP do_fast_isinfinite (SEXP call, SEXP op, SEXP x, SEXP rho,
for (i = 0; i < n; i++) {
xr = COMPLEX(x)[i].r;
xi = COMPLEX(x)[i].i;
if (!ISNAN(xr) && !R_FINITE(xr)
|| !ISNAN(xi) && !R_FINITE(xi)) { ret = TRUE; goto vret; }
if (R_INFINITE(xr) || R_INFINITE(xi)) { ret = TRUE; goto vret; }
}
ret = FALSE; goto vret;
}
for (i = 0; i < n; i++) {
xr = COMPLEX(x)[i].r;
xi = COMPLEX(x)[i].i;
LOGICAL(ans)[i] =
!ISNAN(xr) && !R_FINITE(xr) || !ISNAN(xi) && !R_FINITE(xi);
LOGICAL(ans)[i] = R_INFINITE(xr) || R_INFINITE(xi);
}
break;
default:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment