IEEE Inf and NaN support for numbers.c

[John W. Eaton]

The following patch adds basic IEEE Inf and NaN support for guile.  This
introduces two new predicates (inf? and nan?) and two new functions
(inf and nan, which return IEEE Infinity and NaN values,
respectively).

Example:

  guile> (/ 0)
  Inf
  guile> (/ 10 0)
  Inf
  guile> (/ 0 0) 
  NaN
  guile> (nan? (/ 0 0))
  #t
  guile> (nan? (/ 1 0))
  #f
  guile> (inf? (/ 1 0))
  #t
  guile> (inf? (/ 0 0))
  #f
  guile> (inf? (/ 1 2))
  #f
  guile> (inf)
  Inf
  guile> (nan)
  NaN
  guile> (+ 1 (nan))
  NaN
  guile> (/ 1 (inf))  
  0.0

etc.

The code for computing Infinity and NaN is borrowed from Octave.

jwe

-- 
www.octave.org        | Unfortunately we were hoplessly optimistic in 1954
www.che.wisc.edu/~jwe | about the problems of debugging FORTRAN programs.
                      |                                       -- J. Backus

ChangeLog:

2002-03-06  John W. Eaton  <address@hidden>

        * configure.in (AC_CHECK_HEADERS): Check for floatingpoint.h
        ieeefp.h, and nan.h.
        (AC_CHECK_FUNCS): Check for finite, isinf, and isnan.

libguile/ChangeLog:

2002-03-06  John W. Eaton  <address@hidden>

        * numbers.h: Conditionally include floatingpoint.h, ieeefp.h, and
        nan.h.  Provide declarations for scm_inf_p, scm_nan_p, scn_inf,
        and scm_nan.
        * numbers.c: [SCO && ! HAVE_ISNAN] (isnan): New function.
        [SCO && ! HAVE_ISINF] (isinf): New function.
        (xisinf, xisnan): New functions.
        (IS_INF): Delete.
        (isfinite): Define in terms of xisinf.
        (scm_inf_p, scm_nan_p): New functions.
        (guile_Inf, guile_NaN): New file-scope vars.
        (guile_ieee_init): New function.
        (scm_inf, scm_nan): New functions.
        (idbl2str): Handle Inf and NaN. Remove funny label and
        corresponding gotos.
        (ALLOW_DIVIDE_BY_ZERO): New macro.
        (scm_divide): Allow division by zero to occur if
        ALLOW_DIVIDE_BY_ZERO is defined.
        Handle bignums and ints as special cases.


Index: configure.in
===================================================================
RCS file: /cvs/guile/guile-core/configure.in,v
retrieving revision 1.183
diff -u -r1.183 configure.in
--- configure.in        4 Mar 2002 22:37:37 -0000       1.183
+++ configure.in        6 Mar 2002 21:38:28 -0000
@@ -468,6 +468,10 @@
 
 AC_REPLACE_FUNCS(inet_aton putenv strerror memmove mkstemp)
 
+AC_CHECK_HEADERS(floatingpoint.h ieeefp.h nan.h)
+
+AC_CHECK_FUNCS(finite isinf isnan)
+
 # When testing for the presence of alloca, we need to add alloca.o
 # explicitly to LIBOBJS to make sure that it is translated to
 # `alloca.lo' for libtool later on.  This can and should be done more cleanly.
Index: libguile/numbers.h
===================================================================
RCS file: /cvs/guile/guile-core/libguile/numbers.h,v
retrieving revision 1.63
diff -u -r1.63 numbers.h
--- libguile/numbers.h  25 Nov 2001 15:12:39 -0000      1.63
+++ libguile/numbers.h  6 Mar 2002 21:38:42 -0000
@@ -49,6 +49,24 @@
 #include "libguile/__scm.h"
 #include "libguile/print.h"
 
+#if defined (HAVE_FLOATINGPOINT_H)
+#include <floatingpoint.h>
+#endif
+
+#if defined (HAVE_IEEEFP_H)
+#include <ieeefp.h>
+#endif
+
+#if defined (HAVE_NAN_H)
+#if defined (SCO)
+#define _IEEE 1
+#endif
+#include <nan.h>
+#if defined (SCO)
+#undef _IEEE
+#endif
+#endif
+
 
 
 /* Immediate Numbers, also known as fixnums
@@ -202,6 +220,10 @@
 SCM_API SCM scm_exact_p (SCM x);
 SCM_API SCM scm_odd_p (SCM n);
 SCM_API SCM scm_even_p (SCM n);
+SCM_API SCM scm_inf_p (SCM n);
+SCM_API SCM scm_nan_p (SCM n);
+SCM_API SCM scm_inf (void);
+SCM_API SCM scm_nan (void);
 SCM_API SCM scm_abs (SCM x);
 SCM_API SCM scm_quotient (SCM x, SCM y);
 SCM_API SCM scm_remainder (SCM x, SCM y);
--- numbers.c   2002/03/06 19:03:02     1.1
+++ numbers.c   2002/03/06 22:30:00
@@ -68,20 +68,24 @@
  */
 #define FLOBUFLEN (10+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
 
-
-/* IS_INF tests its floating point number for infiniteness
-   Dirk:FIXME:: This test does not work if x == 0
- */
-#ifndef IS_INF
-#define IS_INF(x) ((x) == (x) / 2)
+#if defined (SCO)
+#if ! defined (HAVE_ISNAN)
+#define HAVE_ISNAN
+static int
+isnan (double x)
+{
+  return (IsNANorINF (x) && NaN (x) && ! IsINF (x)) ? 1 : 0;
+}
 #endif
+#if ! defined (HAVE_ISINF)
+#define HAVE_ISINF
+static int
+isinf (double x)
+{
+  return (IsNANorINF (x) && IsINF (x)) ? 1 : 0;
+}
 
-
-/* Return true if X is not infinite and is not a NaN
-   Dirk:FIXME:: Since IS_INF is broken, this test does not work if x == 0
- */
-#ifndef isfinite
-#define isfinite(x) (!IS_INF (x) && (x) == (x))
+#endif
 #endif
 
 
@@ -141,6 +145,140 @@
 }
 #undef FUNC_NAME
 
+static int
+xisinf (double x)
+{
+#if defined (HAVE_ISINF)
+  return isinf (x);
+#elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
+  return (! (finite (x) || isnan (x)));
+#else
+  return 0;
+#endif
+}
+
+static int
+xisnan (double x)
+{
+#if defined (HAVE_ISNAN)
+  return isnan (x);
+#else
+  return 0;
+#endif
+}
+
+#define isfinite(x) (! xisinf (x))
+
+SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0, 
+            (SCM n),
+           "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
+           "otherwise.")
+#define FUNC_NAME s_scm_inf_p
+{
+  if (SCM_REALP (n)) {
+    return SCM_BOOL (xisinf (SCM_REAL_VALUE (n)));
+  } else if (SCM_COMPLEXP (n)) {
+    return SCM_BOOL (xisinf (SCM_COMPLEX_REAL (n))
+                    || xisinf (SCM_COMPLEX_IMAG (n)));
+  } else {
+    return SCM_BOOL_F;
+  }
+}
+#undef FUNC_NAME
+
+SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0, 
+            (SCM n),
+           "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
+           "otherwise.")
+#define FUNC_NAME s_scm_nan_p
+{
+  if (SCM_REALP (n)) {
+    return SCM_BOOL (xisnan (SCM_REAL_VALUE (n)));
+  } else if (SCM_COMPLEXP (n)) {
+    return SCM_BOOL (xisnan (SCM_COMPLEX_REAL (n))
+                    || xisnan (SCM_COMPLEX_IMAG (n)));
+  } else {
+    return SCM_BOOL_F;
+  }
+}
+#undef FUNC_NAME
+
+/* Guile's idea of infinity.  */
+static double guile_Inf;
+
+/* Guile's idea of not a number.  */
+static double guile_NaN;
+
+static void
+guile_ieee_init (void)
+{
+#if defined (HAVE_ISINF) || defined (HAVE_FINITE)
+
+/* Some version of gcc on some old version of Linux used to crash when
+   trying to make Inf and NaN.  */
+
+#if defined (SCO)
+  double tmp = 1.0;
+  guile_Inf = 1.0 / (tmp - tmp);
+#elif defined (__alpha__) && ! defined (linux)
+  extern unsigned int DINFINITY[2];
+  guile_Inf = (*(X_CAST(double *, DINFINITY)));
+#else
+  double tmp = 1e+10;
+  guile_Inf = tmp;
+  for (;;)
+    {
+      guile_Inf *= 1e+10;
+      if (guile_Inf == tmp)
+       break;
+      tmp = guile_Inf;
+    }
+#endif
+
+#endif
+
+#if defined (HAVE_ISNAN)
+
+#if defined (__alpha__) && ! defined (linux)
+  extern unsigned int DQNAN[2];
+  guile_NaN =  (*(X_CAST(double *, DQNAN)));
+#else
+  guile_NaN = guile_Inf / guile_Inf;
+#endif
+
+#endif
+}
+
+SCM_DEFINE (scm_inf, "inf", 0, 0, 0, 
+            (void),
+           "Return Inf.")
+#define FUNC_NAME s_scm_inf
+{
+  static int initialized = 0;
+  if (! initialized)
+    {
+      guile_ieee_init ();
+      initialized = 1;
+    }
+  return scm_make_real (guile_Inf);
+}
+#undef FUNC_NAME
+
+SCM_DEFINE (scm_nan, "nan", 0, 0, 0, 
+            (void),
+           "Return NaN.")
+#define FUNC_NAME s_scm_nan
+{
+  static int initialized = 0;
+  if (! initialized)
+    {
+      guile_ieee_init ();
+      initialized = 1;
+    }
+  return scm_make_real (guile_NaN);
+}
+#undef FUNC_NAME
+
 
 SCM_GPROC (s_abs, "abs", 1, 0, 0, scm_abs, g_abs);
 /* "Return the absolute value of @var{x}."
@@ -1935,32 +2073,45 @@
       f = -f;
       a[ch++] = '-';
     }
-  else if (f > 0.0);
-  else
-    goto funny;
-  if (IS_INF (f))
+
+  if (xisinf (f))
     {
-      if (ch == 0)
-       a[ch++] = '+';
-    funny:
-      a[ch++] = '#';
-      a[ch++] = '.';
-      a[ch++] = '#';
+      a[ch++] = 'I';
+      a[ch++] = 'n';
+      a[ch++] = 'f';
       return ch;
     }
+  else if (xisnan (f))
+    {
+      a[ch++] = 'N';
+      a[ch++] = 'a';
+      a[ch++] = 'N';
+      return ch;
+    }
+
 #ifdef DBL_MIN_10_EXP  /* Prevent unnormalized values, as from 
                          make-uniform-vector, from causing infinite loops. */
   while (f < 1.0)
     {
       f *= 10.0;
       if (exp-- < DBL_MIN_10_EXP)
-       goto funny;
+       {
+         a[ch++] = '#';
+         a[ch++] = '.';
+         a[ch++] = '#';
+         return ch;
+       }
     }
   while (f > 10.0)
     {
       f *= 0.10;
       if (exp++ > DBL_MAX_10_EXP)
-       goto funny;
+       {
+         a[ch++] = '#';
+         a[ch++] = '.';
+         a[ch++] = '#';
+         return ch;
+       }
     }
 #else
   while (f < 1.0)
@@ -3690,6 +3841,10 @@
 }
 #undef FUNC_NAME
 
+#if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
+     || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
+#define ALLOW_DIVIDE_BY_ZERO
+#endif
 
 /* The code below for complex division is adapted from the GNU
    libstdc++, which is adapted it from f2c's libF77, and is subject to
@@ -3735,8 +3890,10 @@
       long xx = SCM_INUM (x);
       if (xx == 1 || xx == -1) {
        return x;
+#ifndef ALLOW_DIVIDE_BY_ZERO
       } else if (xx == 0) {
        scm_num_overflow (s_divide);
+#endif
       } else {
        return scm_make_real (1.0 / (double) xx);
       }
@@ -3744,9 +3901,11 @@
       return scm_make_real (1.0 / scm_i_big2dbl (x));
     } else if (SCM_REALP (x)) {
       double xx = SCM_REAL_VALUE (x);
-      if (xx == 0.0)
+#ifndef ALLOW_DIVIDE_BY_ZERO
+      if (xx == 0.0) {
        scm_num_overflow (s_divide);
       else
+#endif
        return scm_make_real (1.0 / xx);
     } else if (SCM_COMPLEXP (x)) {
       double r = SCM_COMPLEX_REAL (x);
@@ -3770,7 +3929,11 @@
     if (SCM_INUMP (y)) {
       long yy = SCM_INUM (y);
       if (yy == 0) {
+#ifndef ALLOW_DIVIDE_BY_ZERO
        scm_num_overflow (s_divide);
+#else
+       return scm_make_real ((double) xx / (double) yy);
+#endif
       } else if (xx % yy != 0) {
        return scm_make_real ((double) xx / (double) yy);
       } else {
@@ -3789,9 +3952,11 @@
       return scm_make_real ((double) xx / scm_i_big2dbl (y));
     } else if (SCM_REALP (y)) {
       double yy = SCM_REAL_VALUE (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
       if (yy == 0.0)
        scm_num_overflow (s_divide);
       else
+#endif
        return scm_make_real ((double) xx / yy);
     } else if (SCM_COMPLEXP (y)) {
       a = xx;
@@ -3816,7 +3981,14 @@
     if (SCM_INUMP (y)) {
       long int yy = SCM_INUM (y);
       if (yy == 0) {
+#ifndef ALLOW_DIVIDE_BY_ZERO
        scm_num_overflow (s_divide);
+#else
+       if (scm_bigcomp (x, scm_i_int2big (0)) == 0)
+         return scm_nan ();
+       else
+         return scm_inf ();
+#endif
       } else if (yy == 1) {
        return x;
       } else {
@@ -3855,9 +4027,11 @@
        : scm_make_real (scm_i_big2dbl (x) / scm_i_big2dbl (y));
     } else if (SCM_REALP (y)) {
       double yy = SCM_REAL_VALUE (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
       if (yy == 0.0)
        scm_num_overflow (s_divide);
       else
+#endif
        return scm_make_real (scm_i_big2dbl (x) / yy);
     } else if (SCM_COMPLEXP (y)) {
       a = scm_i_big2dbl (x);
@@ -3869,18 +4043,21 @@
     double rx = SCM_REAL_VALUE (x);
     if (SCM_INUMP (y)) {
       long int yy = SCM_INUM (y);
-      if (yy == 0) {
+#ifndef ALLOW_DIVIDE_BY_ZERO
+      if (yy == 0)
        scm_num_overflow (s_divide);
-      } else {
+      else
+#endif
        return scm_make_real (rx / (double) yy);
-      }
     } else if (SCM_BIGP (y)) {
       return scm_make_real (rx / scm_i_big2dbl (y));
     } else if (SCM_REALP (y)) {
       double yy = SCM_REAL_VALUE (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
       if (yy == 0.0)
        scm_num_overflow (s_divide);
       else
+#endif
        return scm_make_real (rx / yy);
     } else if (SCM_COMPLEXP (y)) {
       a = rx;
@@ -3893,9 +4070,12 @@
     double ix = SCM_COMPLEX_IMAG (x);
     if (SCM_INUMP (y)) {
       long int yy = SCM_INUM (y);
-      if (yy == 0) {
+#ifndef ALLOW_DIVIDE_BY_ZERO
+      if (yy == 0)
        scm_num_overflow (s_divide);
-      } else {
+      else
+#endif
+      {
        double d = yy;
        return scm_make_complex (rx / d, ix / d);
       }
@@ -3904,9 +4084,11 @@
       return scm_make_complex (rx / d, ix / d);
     } else if (SCM_REALP (y)) {
       double yy = SCM_REAL_VALUE (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
       if (yy == 0.0)
        scm_num_overflow (s_divide);
       else
+#endif
        return scm_make_complex (rx / yy, ix / yy);
     } else if (SCM_COMPLEXP (y)) {
       double ry = SCM_COMPLEX_REAL (y);