Re-worked NaN and Inf support
diff --git a/trionan.c b/trionan.c
new file mode 100644
index 0000000..d8bd99b
--- /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