Re-worked NaN and Inf support

1 files changed, 538 insertions, 0 deletions

diff --git a/trionan.c b/trionan.c
new file mode 100644
index 00000000..d8bd99b2
--- /dev/null
+++ b/trionan.c
@@ -0,0 +1,538 @@
+/*************************************************************************
+ *
+ * $Id$
+ *
+ * Copyright (C) 2001 Bjorn Reese <breese@users.sourceforge.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
+ * MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE AUTHORS AND
+ * CONTRIBUTORS ACCEPT NO RESPONSIBILITY IN ANY CONCEIVABLE MANNER.
+ *
+ ************************************************************************
+ *
+ * Functions to handle special quantities in floating-point numbers
+ * (that is, NaNs and infinity). They provide the capability to detect
+ * and fabricate special quantities.
+ *
+ * Although written to be as portable as possible, it can never be
+ * guaranteed to work on all platforms, as not all hardware supports
+ * special quantities.
+ *
+ * The approach used here (approximately) is to:
+ *
+ *   1. Use C99 functionality when available.
+ *   2. Use IEEE 754 bit-patterns if possible.
+ *   3. Use platform-specific techniques.
+ *
+ * This program has been tested on the following platforms (in
+ * alphabetic order)
+ *
+ *   OS              CPU          Compiler
+ * -------------------------------------------------
+ *   AIX 4.1.4       PowerPC      gcc
+ *   Darwin 1.3.7    PowerPC      gcc
+ *   FreeBSD 2.2     x86          gcc
+ *   FreeBSD 3.3     x86          gcc
+ *   FreeBSD 4.3     x86          gcc
+ *   FreeBSD 4.3     Alpha        gcc
+ *   HP-UX 10.20     PA-RISC      gcc
+ *   HP-UX 10.20     PA-RISC      HP C++
+ *   IRIX 6.5        MIPS         MIPSpro C
+ *   Linux 2.2       x86          gcc
+ *   Linux 2.2       Alpha        gcc
+ *   Linux 2.4       IA64         gcc
+ *   Linux 2.4       StrongARM    gcc
+ *   NetBSD 1.4      x86          gcc
+ *   NetBSD 1.4      StrongARM    gcc
+ *   NetBSD 1.5      Alpha        gcc
+ *   RISC OS 4       StrongARM    Norcroft C
+ *   Solaris 2.5.1   x86          gcc
+ *   Solaris 2.5.1   Sparc        gcc
+ *   Solaris 2.6     Sparc        WorkShop 4.2
+ *   Solaris 8       Sparc        Forte C 6
+ *   Tru64 4.0D      Alpha        gcc
+ *   Tru64 5.1       Alpha        gcc
+ *   WinNT           x86          MSVC 5.0 & 6.0
+ *
+ ************************************************************************/
+
+static const char rcsid[] = "@(#)$Id$";
+
+
+/*************************************************************************
+ * Include files
+ */
+#include "triodef.h"
+#include "trionan.h"
+
+#include <math.h>
+#include <string.h>
+#include <limits.h>
+#include <float.h>
+#if defined(TRIO_PLATFORM_UNIX)
+# include <signal.h>
+#endif
+#include <assert.h>
+
+#ifndef __STDC__
+# define volatile
+# define const
+#endif
+
+/*************************************************************************
+ * Definitions
+ */
+
+/* We must enable IEEE floating-point on Alpha */
+#if defined(__alpha) && !defined(_IEEE_FP)
+# if defined(TRIO_COMPILER_DECC)
+#  error "Must be compiled with option -ieee"
+# elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))
+#  error "Must be compiled with option -mieee"
+# endif
+#endif /* __alpha && ! _IEEE_FP */
+
+/*
+ * In ANSI/IEEE 754-1985 64-bits double format numbers have the
+ * following properties (amoungst others)
+ *
+ *   o FLT_RADIX == 2: binary encoding
+ *   o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used
+ *     to indicate special numbers (e.g. NaN and Infinity), so the
+ *     maximum exponent is 10 bits wide (2^10 == 1024).
+ *   o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because
+ *     numbers are normalized the initial binary 1 is represented
+ *     implictly (the so-called "hidden bit"), which leaves us with
+ *     the ability to represent 53 bits wide mantissa.
+ */
+#if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)
+# define USE_IEEE_754
+#endif
+
+
+/*************************************************************************
+ * Data
+ */
+
+#if defined(USE_IEEE_754)
+
+/*
+ * Endian-agnostic indexing macro.
+ *
+ * The value of internalEndianMagic, when converted into a 64-bit
+ * integer, becomes 0x0001020304050607 (we could have used a 64-bit
+ * integer value instead of a double, but not all platforms supports
+ * that type). The value is automatically encoded with the correct
+ * endianess by the compiler, which means that we can support any
+ * kind of endianess. The individual bytes are then used as an index
+ * for the IEEE 754 bit-patterns and masks.
+ */
+#define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[(x)])
+
+static const double internalEndianMagic = 1.4015997730788920e-309;
+
+/* Mask for the exponent */
+static const unsigned char ieee_754_exponent_mask[] = {
+  0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
+};
+
+/* Mask for the mantissa */
+static const unsigned char ieee_754_mantissa_mask[] = {
+  0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
+};
+
+/* Bit-pattern for infinity */
+static const unsigned char ieee_754_infinity_array[] = {
+  0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
+};
+
+/* Bit-pattern for quiet NaN */
+static const unsigned char ieee_754_qnan_array[] = {
+  0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
+};
+
+
+/*************************************************************************
+ * trio_make_double
+ */
+static double
+trio_make_double(const unsigned char *values)
+{
+  volatile double result;
+  int i;
+
+  for (i = 0; i < (int)sizeof(double); i++) {
+    ((volatile unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
+  }
+  return result;
+}
+
+/*************************************************************************
+ * trio_examine_double
+ */
+static int
+trio_is_special_quantity(double number,
+			 int *has_mantissa)
+{
+  unsigned int i;
+  unsigned char current;
+  int is_special_quantity = (1 == 1);
+
+  *has_mantissa = 0;
+
+  for (i = 0; i < (unsigned int)sizeof(double); i++) {
+    current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
+    is_special_quantity
+      &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
+    *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
+  }
+  return is_special_quantity;
+}
+
+#endif /* USE_IEEE_754 */
+
+
+/*************************************************************************
+ * trio_pinf
+ */
+double
+trio_pinf(void)
+{
+  /* Cache the result */
+  static double result = 0.0;
+
+  if (result == 0.0) {
+    
+#if defined(INFINITY) && defined(__STDC_IEC_559__)
+    result = (double)INFINITY;
+
+#elif defined(USE_IEEE_754)
+    result = trio_make_double(ieee_754_infinity_array);
+
+#else
+    /*
+     * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
+     * as infinity. Otherwise we have to resort to an overflow
+     * operation to generate infinity.
+     */
+# if defined(TRIO_PLATFORM_UNIX)
+    void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
+# endif
+
+    result = HUGE_VAL;
+    if (HUGE_VAL == DBL_MAX) {
+      /* Force overflow */
+      result += HUGE_VAL;
+    }
+    
+# if defined(TRIO_PLATFORM_UNIX)
+    signal(SIGFPE, signal_handler);
+# endif
+
+#endif
+  }
+  return result;
+}
+
+/*************************************************************************
+ * trio_ninf
+ */
+double
+trio_ninf(void)
+{
+  static double result = 0.0;
+
+  if (result == 0.0) {
+    /*
+     * Negative infinity is calculated by negating positive infinity,
+     * which can be done because it is legal to do calculations on
+     * infinity (for example,  1 / infinity == 0).
+     */
+    result = -trio_pinf();
+  }
+  return result;
+}
+
+/*************************************************************************
+ * trio_nan
+ */
+double
+trio_nan(void)
+{
+  /* Cache the result */
+  static double result = 0.0;
+
+  if (result == 0.0) {
+    
+#if defined(TRIO_COMPILER_SUPPORTS_C99)
+    result = nan(NULL);
+
+#elif defined(NAN) && defined(__STDC_IEC_559__)
+    result = (double)NAN;
+  
+#elif defined(USE_IEEE_754)
+    result = trio_make_double(ieee_754_qnan_array);
+
+#else
+    /*
+     * There are several ways to generate NaN. The one used here is
+     * to divide infinity by infinity. I would have preferred to add
+     * negative infinity to positive infinity, but that yields wrong
+     * result (infinity) on FreeBSD.
+     *
+     * This may fail if the hardware does not support NaN, or if
+     * the Invalid Operation floating-point exception is unmasked.
+     */
+# if defined(TRIO_PLATFORM_UNIX)
+    void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
+# endif
+    
+    result = trio_pinf() / trio_pinf();
+    
+# if defined(TRIO_PLATFORM_UNIX)
+    signal(SIGFPE, signal_handler);
+# endif
+    
+#endif
+  }
+  return result;
+}
+
+/*************************************************************************
+ * trio_isnan
+ */
+int
+trio_isnan(volatile double number)
+{
+#if defined(isnan) || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
+  /*
+   * C99 defines isnan() as a macro. UNIX95 defines isnan() as a
+   * function. This function was already present in XPG4, but this
+   * is a bit tricky to detect with compiler defines, so we choose
+   * the conservative approach and only use it for UNIX95.
+   */
+  return isnan(number);
+  
+#elif defined(TRIO_COMPILER_MSVC)
+  /*
+   * MSC has an _isnan() function
+   */
+  return _isnan(number);
+
+#elif defined(USE_IEEE_754)
+  /*
+   * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
+   * pattern, and a non-empty mantissa.
+   */
+  int has_mantissa;
+  int is_special_quantity;
+
+  is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
+  
+  return (is_special_quantity && has_mantissa);
+  
+#else
+  /*
+   * Fallback solution
+   */
+  int status;
+  double integral, fraction;
+  
+# if defined(TRIO_PLATFORM_UNIX)
+  void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
+# endif
+  
+  status = (/*
+	     * NaN is the only number which does not compare to itself
+	     */
+	    (number != number) ||
+	    /*
+	     * Fallback solution if NaN compares to NaN
+	     */
+	    ((number != 0.0) &&
+	     (fraction = modf(number, &integral),
+	      integral == fraction)));
+  
+# if defined(TRIO_PLATFORM_UNIX)
+  signal(SIGFPE, signal_handler);
+# endif
+  
+  return status;
+  
+#endif
+}
+
+/*************************************************************************
+ * trio_isinf
+ */
+int
+trio_isinf(volatile double number)
+{
+#if defined(TRIO_COMPILER_DECC)
+  /*
+   * DECC has an isinf() macro, but it works differently than that
+   * of C99, so we use the fp_class() function instead.
+   */
+  return ((fp_class(number) == FP_POS_INF)
+	  ? 1
+	  : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));
+  
+#elif defined(isinf)
+  /*
+   * C99 defines isinf() as a macro.
+   */
+  return isinf(number);
+  
+#elif defined(TRIO_COMPILER_MSVC)
+  /*
+   * MSVC has an _fpclass() function that can be used to detect infinity.
+   */
+  return ((_fpclass(number) == _FPCLASS_PINF)
+	  ? 1
+	  : ((_fpclass(number) == _FPCLASS_NINF) ? -1 : 0));
+
+#elif defined(USE_IEEE_754)
+  /*
+   * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
+   * pattern, and an empty mantissa.
+   */
+  int has_mantissa;
+  int is_special_quantity;
+
+  is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
+  
+  return (is_special_quantity && !has_mantissa)
+    ? ((number < 0.0) ? -1 : 1)
+    : 0;
+
+#else
+  /*
+   * Fallback solution.
+   */
+  int status;
+  
+# if defined(TRIO_PLATFORM_UNIX)
+  void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
+# endif
+  
+  double infinity = trio_pinf();
+  
+  status = ((number == infinity)
+	    ? 1
+	    : ((number == -infinity) ? -1 : 0));
+  
+# if defined(TRIO_PLATFORM_UNIX)
+  signal(SIGFPE, signal_handler);
+# endif
+  
+  return status;
+  
+#endif
+}
+
+/*************************************************************************
+ */
+#if defined(STANDALONE)
+# include <stdio.h>
+
+int main(void)
+{
+  double my_nan;
+  double my_pinf;
+  double my_ninf;
+# if defined(TRIO_PLATFORM_UNIX)
+  void (*signal_handler)(int);
+# endif
+
+  my_nan = trio_nan();
+  my_pinf = trio_pinf();
+  my_ninf = trio_ninf();
+
+  printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_nan,
+	 ((unsigned char *)&my_nan)[0],
+	 ((unsigned char *)&my_nan)[1],
+	 ((unsigned char *)&my_nan)[2],
+	 ((unsigned char *)&my_nan)[3],
+	 ((unsigned char *)&my_nan)[4],
+	 ((unsigned char *)&my_nan)[5],
+	 ((unsigned char *)&my_nan)[6],
+	 ((unsigned char *)&my_nan)[7],
+	 trio_isnan(my_nan), trio_isinf(my_nan));
+  printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_pinf,
+	 ((unsigned char *)&my_pinf)[0],
+	 ((unsigned char *)&my_pinf)[1],
+	 ((unsigned char *)&my_pinf)[2],
+	 ((unsigned char *)&my_pinf)[3],
+	 ((unsigned char *)&my_pinf)[4],
+	 ((unsigned char *)&my_pinf)[5],
+	 ((unsigned char *)&my_pinf)[6],
+	 ((unsigned char *)&my_pinf)[7],
+	 trio_isnan(my_pinf), trio_isinf(my_pinf));
+  printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_ninf,
+	 ((unsigned char *)&my_ninf)[0],
+	 ((unsigned char *)&my_ninf)[1],
+	 ((unsigned char *)&my_ninf)[2],
+	 ((unsigned char *)&my_ninf)[3],
+	 ((unsigned char *)&my_ninf)[4],
+	 ((unsigned char *)&my_ninf)[5],
+	 ((unsigned char *)&my_ninf)[6],
+	 ((unsigned char *)&my_ninf)[7],
+	 trio_isnan(my_ninf), trio_isinf(my_ninf));
+  
+# if defined(TRIO_PLATFORM_UNIX)
+  signal_handler = signal(SIGFPE, SIG_IGN);
+# endif
+  
+  my_pinf = DBL_MAX + DBL_MAX;
+  my_ninf = -my_pinf;
+  my_nan = my_pinf / my_pinf;
+
+# if defined(TRIO_PLATFORM_UNIX)
+  signal(SIGFPE, signal_handler);
+# endif
+  
+  printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_nan,
+	 ((unsigned char *)&my_nan)[0],
+	 ((unsigned char *)&my_nan)[1],
+	 ((unsigned char *)&my_nan)[2],
+	 ((unsigned char *)&my_nan)[3],
+	 ((unsigned char *)&my_nan)[4],
+	 ((unsigned char *)&my_nan)[5],
+	 ((unsigned char *)&my_nan)[6],
+	 ((unsigned char *)&my_nan)[7],
+	 trio_isnan(my_nan), trio_isinf(my_nan));
+  printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_pinf,
+	 ((unsigned char *)&my_pinf)[0],
+	 ((unsigned char *)&my_pinf)[1],
+	 ((unsigned char *)&my_pinf)[2],
+	 ((unsigned char *)&my_pinf)[3],
+	 ((unsigned char *)&my_pinf)[4],
+	 ((unsigned char *)&my_pinf)[5],
+	 ((unsigned char *)&my_pinf)[6],
+	 ((unsigned char *)&my_pinf)[7],
+	 trio_isnan(my_pinf), trio_isinf(my_pinf));
+  printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
+	 my_ninf,
+	 ((unsigned char *)&my_ninf)[0],
+	 ((unsigned char *)&my_ninf)[1],
+	 ((unsigned char *)&my_ninf)[2],
+	 ((unsigned char *)&my_ninf)[3],
+	 ((unsigned char *)&my_ninf)[4],
+	 ((unsigned char *)&my_ninf)[5],
+	 ((unsigned char *)&my_ninf)[6],
+	 ((unsigned char *)&my_ninf)[7],
+	 trio_isnan(my_ninf), trio_isinf(my_ninf));
+	 
+  return 0;
+}
+#endif