/*************************************************************************

 *

 *

 * 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.

 *

 ************************************************************************/

/*

 * TODO:

 *  o Put all the magic into trio_fpclassify_and_signbit(), and use this from

 *    trio_isnan() etc.

 */

/*************************************************************************

 * Include files

 */

#include "XmlEngtriodef.h"

#include "XmlEngtrionan.h"

#include <math.h>

#include <string.h>

#include <limits.h>

// TODO: XE Cleanup

#if defined(TRIO_PLATFORM_UNIX)

# include <signal.h>

#endif

#if defined(TRIO_COMPILER_DECC)

#  if defined(__linux__)

#   include <cpml.h>

#  else

#   include <fp_class.h>

#  endif

#endif

#include <assert.h>

#if defined(TRIO_DOCUMENTATION)

# include "doc/doc_nan.h"

#endif

/** @addtogroup SpecialQuantities

    @{

*/

/*************************************************************************

 * Definitions

 */

#define TRIO_TRUE (1 == 1)

#define TRIO_FALSE (0 == 1)

/*

 * We must enable IEEE floating-point on Alpha

 */

 // TODO: XE Cleanup

#if defined(__alpha) && !defined(_IEEE_FP)

# if defined(TRIO_COMPILER_DECC)

#  if defined(TRIO_PLATFORM_VMS)

#   error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"

#  else

#   if !defined(_CFE)

#    error "Must be compiled with option -ieee"

#   endif

#  endif

# 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

 *     implicitly (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

#ifdef WIN32

# define USE_IEEE_754

#endif

#ifndef DBL_MAX

#define DBL_MAX 1.79769313486231500e+308

#endif

#ifndef DBL_MIN

#define DBL_MIN 2.22507385850720200e-308

#endif

/*************************************************************************

 * Constants

 */

#if defined(USE_IEEE_754)

/*

 * Endian-agnostic indexing macro.

 *

 * The value of internalEndianMagic, when converted into a 64-bit

 * integer, becomes 0x0706050403020100 (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)[7-(x)])

static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;

/* Mask for the exponent */

static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {

  0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00

};

/* Mask for the mantissa */

static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {

  0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF

};

/* Mask for the sign bit */

static TRIO_CONST unsigned char ieee_754_sign_mask[] = {

  0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00

};

/* Bit-pattern for negative zero */

static TRIO_CONST unsigned char ieee_754_negzero_array[] = {

  0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00

};

/* Bit-pattern for infinity */

static TRIO_CONST unsigned char ieee_754_infinity_array[] = {

  0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00

};

/* Bit-pattern for quiet NaN */

static TRIO_CONST unsigned char ieee_754_qnan_array[] = {

  0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00

};

/*************************************************************************

 * Functions

 */

/*

 * trio_make_double

 */

TRIO_PRIVATE double

trio_make_double

TRIO_ARGS1((values),

	   TRIO_CONST unsigned char *values)

{

  TRIO_VOLATILE double result;

  int i;

  for (i = 0; i < (int)sizeof(double); i++) {

    ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];

  }

  return result;

}

/*

 * trio_is_special_quantity

 */

TRIO_PRIVATE int

trio_is_special_quantity

TRIO_ARGS2((number, has_mantissa),

	   double number,

	   int *has_mantissa)

{

  unsigned int i;

  unsigned char current;

  int is_special_quantity = TRIO_TRUE;

  *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;

}

/*

 * trio_is_negative

 */

TRIO_PRIVATE int

trio_is_negative

TRIO_ARGS1((number),

	   double number)

{

  unsigned int i;

  int is_negative = TRIO_FALSE;

  for (i = 0; i < (unsigned int)sizeof(double); i++) {

    is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]

		    & ieee_754_sign_mask[i]);

  }

  return is_negative;

}

#endif /* USE_IEEE_754 */

/**

   Generate negative zero.

   @return Floating-point representation of negative zero.

*/

TRIO_PUBLIC double

trio_nzero(TRIO_NOARGS)

{

#if defined(USE_IEEE_754)

  return trio_make_double(ieee_754_negzero_array);

#else

  TRIO_VOLATILE double zero = 0.0;

  return -zero;

#endif

}

/**

   Generate positive infinity.

   @return Floating-point representation of positive infinity.

*/

TRIO_PUBLIC double

trio_pinf(TRIO_NOARGS)

{

  /* Cache the result */

	//FIXIT

//  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;

}

/**

   Generate negative infinity.

   @return Floating-point value of negative infinity.

*/

TRIO_PUBLIC double

trio_ninf(TRIO_NOARGS)

{

	//FIXIT

  //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;

}

/**

   Generate NaN.

   @return Floating-point representation of NaN.

*/

TRIO_PUBLIC double

trio_nan(TRIO_NOARGS)

{

  /* Cache the result */

	//FIXIT

  //static 

	double result = 0.0;

  if (result == 0.0) {

#if defined(TRIO_COMPILER_SUPPORTS_C99)

    result = nan("");

#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;

}

/**

   Check for NaN.

   @param number An arbitrary floating-point number.

   @return Boolean value indicating whether or not the number is a NaN.

*/

TRIO_PUBLIC int

trio_isnan

TRIO_ARGS1((number),

	   double number)

{

#if (defined(TRIO_COMPILER_SUPPORTS_C99) && 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) || defined(TRIO_COMPILER_BCB)

  /*

   * Microsoft Visual C++ and Borland C++ Builder have an _isnan()

   * function.

   */

  return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;

#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

	     */

	    ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)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

}

/**

   Check for infinity.

   @param number An arbitrary floating-point number.

   @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.

*/

TRIO_PUBLIC int

trio_isinf

TRIO_ARGS1((number),

	   double number)

{

#if defined(TRIO_COMPILER_DECC) && !defined(__linux__)

  /*

   * 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) && !defined(WIN32)

  /*

   * C99 defines isinf() as a macro.

   */

  return isinf(number)

    ? ((number > 0.0) ? 1 : -1)

    : 0;

#elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)

  /*

   * Microsoft Visual C++ and Borland C++ Builder have 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 0

	/* Temporary fix - this routine is not used anywhere */

/**

   Check for finity.