Home | History | Annotate | Download | only in src
      1 /*************************************************************************
      2  *
      3  * $Id$
      4  *
      5  * Copyright (C) 2001 Bjorn Reese <breese (at) users.sourceforge.net>
      6  *
      7  * Permission to use, copy, modify, and distribute this software for any
      8  * purpose with or without fee is hereby granted, provided that the above
      9  * copyright notice and this permission notice appear in all copies.
     10  *
     11  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
     12  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
     13  * MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE AUTHORS AND
     14  * CONTRIBUTORS ACCEPT NO RESPONSIBILITY IN ANY CONCEIVABLE MANNER.
     15  *
     16  ************************************************************************
     17  *
     18  * Functions to handle special quantities in floating-point numbers
     19  * (that is, NaNs and infinity). They provide the capability to detect
     20  * and fabricate special quantities.
     21  *
     22  * Although written to be as portable as possible, it can never be
     23  * guaranteed to work on all platforms, as not all hardware supports
     24  * special quantities.
     25  *
     26  * The approach used here (approximately) is to:
     27  *
     28  *   1. Use C99 functionality when available.
     29  *   2. Use IEEE 754 bit-patterns if possible.
     30  *   3. Use platform-specific techniques.
     31  *
     32  ************************************************************************/
     33 
     34 /*
     35  * TODO:
     36  *  o Put all the magic into trio_fpclassify_and_signbit(), and use this from
     37  *    trio_isnan() etc.
     38  */
     39 
     40 /*************************************************************************
     41  * Include files
     42  */
     43 #include "triodef.h"
     44 #include "trionan.h"
     45 
     46 #include <math.h>
     47 #include <string.h>
     48 #include <limits.h>
     49 #include <float.h>
     50 #if defined(TRIO_PLATFORM_UNIX)
     51 # include <signal.h>
     52 #endif
     53 #if defined(TRIO_COMPILER_DECC)
     54 #  if defined(__linux__)
     55 #   include <cpml.h>
     56 #  else
     57 #   include <fp_class.h>
     58 #  endif
     59 #endif
     60 #include <assert.h>
     61 
     62 #if defined(TRIO_DOCUMENTATION)
     63 # include "doc/doc_nan.h"
     64 #endif
     65 /** @addtogroup SpecialQuantities
     66     @{
     67 */
     68 
     69 /*************************************************************************
     70  * Definitions
     71  */
     72 
     73 #define TRIO_TRUE (1 == 1)
     74 #define TRIO_FALSE (0 == 1)
     75 
     76 /*
     77  * We must enable IEEE floating-point on Alpha
     78  */
     79 #if defined(__alpha) && !defined(_IEEE_FP)
     80 # if defined(TRIO_COMPILER_DECC)
     81 #  if defined(TRIO_PLATFORM_VMS)
     82 #   error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"
     83 #  else
     84 #   if !defined(_CFE)
     85 #    error "Must be compiled with option -ieee"
     86 #   endif
     87 #  endif
     88 # elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))
     89 #  error "Must be compiled with option -mieee"
     90 # endif
     91 #endif /* __alpha && ! _IEEE_FP */
     92 
     93 /*
     94  * In ANSI/IEEE 754-1985 64-bits double format numbers have the
     95  * following properties (amoungst others)
     96  *
     97  *   o FLT_RADIX == 2: binary encoding
     98  *   o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used
     99  *     to indicate special numbers (e.g. NaN and Infinity), so the
    100  *     maximum exponent is 10 bits wide (2^10 == 1024).
    101  *   o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because
    102  *     numbers are normalized the initial binary 1 is represented
    103  *     implicitly (the so-called "hidden bit"), which leaves us with
    104  *     the ability to represent 53 bits wide mantissa.
    105  */
    106 #if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)
    107 # define USE_IEEE_754
    108 #endif
    109 
    110 
    111 /*************************************************************************
    112  * Constants
    113  */
    114 
    115 static TRIO_CONST char rcsid[] = "@(#)$Id$";
    116 
    117 #if defined(USE_IEEE_754)
    118 
    119 /*
    120  * Endian-agnostic indexing macro.
    121  *
    122  * The value of internalEndianMagic, when converted into a 64-bit
    123  * integer, becomes 0x0706050403020100 (we could have used a 64-bit
    124  * integer value instead of a double, but not all platforms supports
    125  * that type). The value is automatically encoded with the correct
    126  * endianess by the compiler, which means that we can support any
    127  * kind of endianess. The individual bytes are then used as an index
    128  * for the IEEE 754 bit-patterns and masks.
    129  */
    130 #define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])
    131 
    132 #if (defined(__BORLANDC__) && __BORLANDC__ >= 0x0590)
    133 static TRIO_CONST double internalEndianMagic = 7.949928895127362e-275;
    134 #else
    135 static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
    136 #endif
    137 
    138 /* Mask for the exponent */
    139 static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
    140   0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
    141 };
    142 
    143 /* Mask for the mantissa */
    144 static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
    145   0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
    146 };
    147 
    148 /* Mask for the sign bit */
    149 static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
    150   0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
    151 };
    152 
    153 /* Bit-pattern for negative zero */
    154 static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
    155   0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
    156 };
    157 
    158 /* Bit-pattern for infinity */
    159 static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
    160   0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
    161 };
    162 
    163 /* Bit-pattern for quiet NaN */
    164 static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
    165   0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
    166 };
    167 
    168 
    169 /*************************************************************************
    170  * Functions
    171  */
    172 
    173 /*
    174  * trio_make_double
    175  */
    176 TRIO_PRIVATE double
    177 trio_make_double
    178 TRIO_ARGS1((values),
    179 	   TRIO_CONST unsigned char *values)
    180 {
    181   TRIO_VOLATILE double result;
    182   int i;
    183 
    184   for (i = 0; i < (int)sizeof(double); i++) {
    185     ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
    186   }
    187   return result;
    188 }
    189 
    190 /*
    191  * trio_is_special_quantity
    192  */
    193 TRIO_PRIVATE int
    194 trio_is_special_quantity
    195 TRIO_ARGS2((number, has_mantissa),
    196 	   double number,
    197 	   int *has_mantissa)
    198 {
    199   unsigned int i;
    200   unsigned char current;
    201   int is_special_quantity = TRIO_TRUE;
    202 
    203   *has_mantissa = 0;
    204 
    205   for (i = 0; i < (unsigned int)sizeof(double); i++) {
    206     current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
    207     is_special_quantity
    208       &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
    209     *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
    210   }
    211   return is_special_quantity;
    212 }
    213 
    214 /*
    215  * trio_is_negative
    216  */
    217 TRIO_PRIVATE int
    218 trio_is_negative
    219 TRIO_ARGS1((number),
    220 	   double number)
    221 {
    222   unsigned int i;
    223   int is_negative = TRIO_FALSE;
    224 
    225   for (i = 0; i < (unsigned int)sizeof(double); i++) {
    226     is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
    227 		    & ieee_754_sign_mask[i]);
    228   }
    229   return is_negative;
    230 }
    231 
    232 #endif /* USE_IEEE_754 */
    233 
    234 
    235 /**
    236    Generate negative zero.
    237 
    238    @return Floating-point representation of negative zero.
    239 */
    240 TRIO_PUBLIC double
    241 trio_nzero(TRIO_NOARGS)
    242 {
    243 #if defined(USE_IEEE_754)
    244   return trio_make_double(ieee_754_negzero_array);
    245 #else
    246   TRIO_VOLATILE double zero = 0.0;
    247 
    248   return -zero;
    249 #endif
    250 }
    251 
    252 /**
    253    Generate positive infinity.
    254 
    255    @return Floating-point representation of positive infinity.
    256 */
    257 TRIO_PUBLIC double
    258 trio_pinf(TRIO_NOARGS)
    259 {
    260   /* Cache the result */
    261   static double result = 0.0;
    262 
    263   if (result == 0.0) {
    264 
    265 #if defined(INFINITY) && defined(__STDC_IEC_559__)
    266     result = (double)INFINITY;
    267 
    268 #elif defined(USE_IEEE_754)
    269     result = trio_make_double(ieee_754_infinity_array);
    270 
    271 #else
    272     /*
    273      * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
    274      * as infinity. Otherwise we have to resort to an overflow
    275      * operation to generate infinity.
    276      */
    277 # if defined(TRIO_PLATFORM_UNIX)
    278     void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
    279 # endif
    280 
    281     result = HUGE_VAL;
    282     if (HUGE_VAL == DBL_MAX) {
    283       /* Force overflow */
    284       result += HUGE_VAL;
    285     }
    286 
    287 # if defined(TRIO_PLATFORM_UNIX)
    288     signal(SIGFPE, signal_handler);
    289 # endif
    290 
    291 #endif
    292   }
    293   return result;
    294 }
    295 
    296 /**
    297    Generate negative infinity.
    298 
    299    @return Floating-point value of negative infinity.
    300 */
    301 TRIO_PUBLIC double
    302 trio_ninf(TRIO_NOARGS)
    303 {
    304   static double result = 0.0;
    305 
    306   if (result == 0.0) {
    307     /*
    308      * Negative infinity is calculated by negating positive infinity,
    309      * which can be done because it is legal to do calculations on
    310      * infinity (for example,  1 / infinity == 0).
    311      */
    312     result = -trio_pinf();
    313   }
    314   return result;
    315 }
    316 
    317 /**
    318    Generate NaN.
    319 
    320    @return Floating-point representation of NaN.
    321 */
    322 TRIO_PUBLIC double
    323 trio_nan(TRIO_NOARGS)
    324 {
    325   /* Cache the result */
    326   static double result = 0.0;
    327 
    328   if (result == 0.0) {
    329 
    330 #if defined(TRIO_COMPILER_SUPPORTS_C99)
    331     result = nan("");
    332 
    333 #elif defined(NAN) && defined(__STDC_IEC_559__)
    334     result = (double)NAN;
    335 
    336 #elif defined(USE_IEEE_754)
    337     result = trio_make_double(ieee_754_qnan_array);
    338 
    339 #else
    340     /*
    341      * There are several ways to generate NaN. The one used here is
    342      * to divide infinity by infinity. I would have preferred to add
    343      * negative infinity to positive infinity, but that yields wrong
    344      * result (infinity) on FreeBSD.
    345      *
    346      * This may fail if the hardware does not support NaN, or if
    347      * the Invalid Operation floating-point exception is unmasked.
    348      */
    349 # if defined(TRIO_PLATFORM_UNIX)
    350     void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
    351 # endif
    352 
    353     result = trio_pinf() / trio_pinf();
    354 
    355 # if defined(TRIO_PLATFORM_UNIX)
    356     signal(SIGFPE, signal_handler);
    357 # endif
    358 
    359 #endif
    360   }
    361   return result;
    362 }
    363 
    364 /**
    365    Check for NaN.
    366 
    367    @param number An arbitrary floating-point number.
    368    @return Boolean value indicating whether or not the number is a NaN.
    369 */
    370 TRIO_PUBLIC int
    371 trio_isnan
    372 TRIO_ARGS1((number),
    373 	   double number)
    374 {
    375 #if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \
    376  || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
    377   /*
    378    * C99 defines isnan() as a macro. UNIX95 defines isnan() as a
    379    * function. This function was already present in XPG4, but this
    380    * is a bit tricky to detect with compiler defines, so we choose
    381    * the conservative approach and only use it for UNIX95.
    382    */
    383   return isnan(number);
    384 
    385 #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
    386   /*
    387    * Microsoft Visual C++ and Borland C++ Builder have an _isnan()
    388    * function.
    389    */
    390   return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;
    391 
    392 #elif defined(USE_IEEE_754)
    393   /*
    394    * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
    395    * pattern, and a non-empty mantissa.
    396    */
    397   int has_mantissa;
    398   int is_special_quantity;
    399 
    400   is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
    401 
    402   return (is_special_quantity && has_mantissa);
    403 
    404 #else
    405   /*
    406    * Fallback solution
    407    */
    408   int status;
    409   double integral, fraction;
    410 
    411 # if defined(TRIO_PLATFORM_UNIX)
    412   void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
    413 # endif
    414 
    415   status = (/*
    416 	     * NaN is the only number which does not compare to itself
    417 	     */
    418 	    ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
    419 	    /*
    420 	     * Fallback solution if NaN compares to NaN
    421 	     */
    422 	    ((number != 0.0) &&
    423 	     (fraction = modf(number, &integral),
    424 	      integral == fraction)));
    425 
    426 # if defined(TRIO_PLATFORM_UNIX)
    427   signal(SIGFPE, signal_handler);
    428 # endif
    429 
    430   return status;
    431 
    432 #endif
    433 }
    434 
    435 /**
    436    Check for infinity.
    437 
    438    @param number An arbitrary floating-point number.
    439    @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.
    440 */
    441 TRIO_PUBLIC int
    442 trio_isinf
    443 TRIO_ARGS1((number),
    444 	   double number)
    445 {
    446 #if defined(TRIO_COMPILER_DECC) && !defined(__linux__)
    447   /*
    448    * DECC has an isinf() macro, but it works differently than that
    449    * of C99, so we use the fp_class() function instead.
    450    */
    451   return ((fp_class(number) == FP_POS_INF)
    452 	  ? 1
    453 	  : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));
    454 
    455 #elif defined(isinf)
    456   /*
    457    * C99 defines isinf() as a macro.
    458    */
    459   return isinf(number)
    460     ? ((number > 0.0) ? 1 : -1)
    461     : 0;
    462 
    463 #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
    464   /*
    465    * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
    466    * function that can be used to detect infinity.
    467    */
    468   return ((_fpclass(number) == _FPCLASS_PINF)
    469 	  ? 1
    470 	  : ((_fpclass(number) == _FPCLASS_NINF) ? -1 : 0));
    471 
    472 #elif defined(USE_IEEE_754)
    473   /*
    474    * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
    475    * pattern, and an empty mantissa.
    476    */
    477   int has_mantissa;
    478   int is_special_quantity;
    479 
    480   is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
    481 
    482   return (is_special_quantity && !has_mantissa)
    483     ? ((number < 0.0) ? -1 : 1)
    484     : 0;
    485 
    486 #else
    487   /*
    488    * Fallback solution.
    489    */
    490   int status;
    491 
    492 # if defined(TRIO_PLATFORM_UNIX)
    493   void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
    494 # endif
    495 
    496   double infinity = trio_pinf();
    497 
    498   status = ((number == infinity)
    499 	    ? 1
    500 	    : ((number == -infinity) ? -1 : 0));
    501 
    502 # if defined(TRIO_PLATFORM_UNIX)
    503   signal(SIGFPE, signal_handler);
    504 # endif
    505 
    506   return status;
    507 
    508 #endif
    509 }
    510 
    511 #if 0
    512 	/* Temporary fix - this routine is not used anywhere */
    513 /**
    514    Check for finity.
    515 
    516    @param number An arbitrary floating-point number.
    517    @return Boolean value indicating whether or not the number is a finite.
    518 */
    519 TRIO_PUBLIC int
    520 trio_isfinite
    521 TRIO_ARGS1((number),
    522 	   double number)
    523 {
    524 #if defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isfinite)
    525   /*
    526    * C99 defines isfinite() as a macro.
    527    */
    528   return isfinite(number);
    529 
    530 #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
    531   /*
    532    * Microsoft Visual C++ and Borland C++ Builder use _finite().
    533    */
    534   return _finite(number);
    535 
    536 #elif defined(USE_IEEE_754)
    537   /*
    538    * Examine IEEE 754 bit-pattern. For finity we do not care about the
    539    * mantissa.
    540    */
    541   int dummy;
    542 
    543   return (! trio_is_special_quantity(number, &dummy));
    544 
    545 #else
    546   /*
    547    * Fallback solution.
    548    */
    549   return ((trio_isinf(number) == 0) && (trio_isnan(number) == 0));
    550 
    551 #endif
    552 }
    553 
    554 #endif
    555 
    556 /*
    557  * The sign of NaN is always false
    558  */
    559 TRIO_PUBLIC int
    560 trio_fpclassify_and_signbit
    561 TRIO_ARGS2((number, is_negative),
    562 	   double number,
    563 	   int *is_negative)
    564 {
    565 #if defined(fpclassify) && defined(signbit)
    566   /*
    567    * C99 defines fpclassify() and signbit() as a macros
    568    */
    569   *is_negative = signbit(number);
    570   switch (fpclassify(number)) {
    571   case FP_NAN:
    572     return TRIO_FP_NAN;
    573   case FP_INFINITE:
    574     return TRIO_FP_INFINITE;
    575   case FP_SUBNORMAL:
    576     return TRIO_FP_SUBNORMAL;
    577   case FP_ZERO:
    578     return TRIO_FP_ZERO;
    579   default:
    580     return TRIO_FP_NORMAL;
    581   }
    582 
    583 #else
    584 # if defined(TRIO_COMPILER_DECC)
    585   /*
    586    * DECC has an fp_class() function.
    587    */
    588 #  define TRIO_FPCLASSIFY(n) fp_class(n)
    589 #  define TRIO_QUIET_NAN FP_QNAN
    590 #  define TRIO_SIGNALLING_NAN FP_SNAN
    591 #  define TRIO_POSITIVE_INFINITY FP_POS_INF
    592 #  define TRIO_NEGATIVE_INFINITY FP_NEG_INF
    593 #  define TRIO_POSITIVE_SUBNORMAL FP_POS_DENORM
    594 #  define TRIO_NEGATIVE_SUBNORMAL FP_NEG_DENORM
    595 #  define TRIO_POSITIVE_ZERO FP_POS_ZERO
    596 #  define TRIO_NEGATIVE_ZERO FP_NEG_ZERO
    597 #  define TRIO_POSITIVE_NORMAL FP_POS_NORM
    598 #  define TRIO_NEGATIVE_NORMAL FP_NEG_NORM
    599 
    600 # elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
    601   /*
    602    * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
    603    * function.
    604    */
    605 #  define TRIO_FPCLASSIFY(n) _fpclass(n)
    606 #  define TRIO_QUIET_NAN _FPCLASS_QNAN
    607 #  define TRIO_SIGNALLING_NAN _FPCLASS_SNAN
    608 #  define TRIO_POSITIVE_INFINITY _FPCLASS_PINF
    609 #  define TRIO_NEGATIVE_INFINITY _FPCLASS_NINF
    610 #  define TRIO_POSITIVE_SUBNORMAL _FPCLASS_PD
    611 #  define TRIO_NEGATIVE_SUBNORMAL _FPCLASS_ND
    612 #  define TRIO_POSITIVE_ZERO _FPCLASS_PZ
    613 #  define TRIO_NEGATIVE_ZERO _FPCLASS_NZ
    614 #  define TRIO_POSITIVE_NORMAL _FPCLASS_PN
    615 #  define TRIO_NEGATIVE_NORMAL _FPCLASS_NN
    616 
    617 # elif defined(FP_PLUS_NORM)
    618   /*
    619    * HP-UX 9.x and 10.x have an fpclassify() function, that is different
    620    * from the C99 fpclassify() macro supported on HP-UX 11.x.
    621    *
    622    * AIX has class() for C, and _class() for C++, which returns the
    623    * same values as the HP-UX fpclassify() function.
    624    */
    625 #  if defined(TRIO_PLATFORM_AIX)
    626 #   if defined(__cplusplus)
    627 #    define TRIO_FPCLASSIFY(n) _class(n)
    628 #   else
    629 #    define TRIO_FPCLASSIFY(n) class(n)
    630 #   endif
    631 #  else
    632 #   define TRIO_FPCLASSIFY(n) fpclassify(n)
    633 #  endif
    634 #  define TRIO_QUIET_NAN FP_QNAN
    635 #  define TRIO_SIGNALLING_NAN FP_SNAN
    636 #  define TRIO_POSITIVE_INFINITY FP_PLUS_INF
    637 #  define TRIO_NEGATIVE_INFINITY FP_MINUS_INF
    638 #  define TRIO_POSITIVE_SUBNORMAL FP_PLUS_DENORM
    639 #  define TRIO_NEGATIVE_SUBNORMAL FP_MINUS_DENORM
    640 #  define TRIO_POSITIVE_ZERO FP_PLUS_ZERO
    641 #  define TRIO_NEGATIVE_ZERO FP_MINUS_ZERO
    642 #  define TRIO_POSITIVE_NORMAL FP_PLUS_NORM
    643 #  define TRIO_NEGATIVE_NORMAL FP_MINUS_NORM
    644 # endif
    645 
    646 # if defined(TRIO_FPCLASSIFY)
    647   switch (TRIO_FPCLASSIFY(number)) {
    648   case TRIO_QUIET_NAN:
    649   case TRIO_SIGNALLING_NAN:
    650     *is_negative = TRIO_FALSE; /* NaN has no sign */
    651     return TRIO_FP_NAN;
    652   case TRIO_POSITIVE_INFINITY:
    653     *is_negative = TRIO_FALSE;
    654     return TRIO_FP_INFINITE;
    655   case TRIO_NEGATIVE_INFINITY:
    656     *is_negative = TRIO_TRUE;
    657     return TRIO_FP_INFINITE;
    658   case TRIO_POSITIVE_SUBNORMAL:
    659     *is_negative = TRIO_FALSE;
    660     return TRIO_FP_SUBNORMAL;
    661   case TRIO_NEGATIVE_SUBNORMAL:
    662     *is_negative = TRIO_TRUE;
    663     return TRIO_FP_SUBNORMAL;
    664   case TRIO_POSITIVE_ZERO:
    665     *is_negative = TRIO_FALSE;
    666     return TRIO_FP_ZERO;
    667   case TRIO_NEGATIVE_ZERO:
    668     *is_negative = TRIO_TRUE;
    669     return TRIO_FP_ZERO;
    670   case TRIO_POSITIVE_NORMAL:
    671     *is_negative = TRIO_FALSE;
    672     return TRIO_FP_NORMAL;
    673   case TRIO_NEGATIVE_NORMAL:
    674     *is_negative = TRIO_TRUE;
    675     return TRIO_FP_NORMAL;
    676   default:
    677     /* Just in case... */
    678     *is_negative = (number < 0.0);
    679     return TRIO_FP_NORMAL;
    680   }
    681 
    682 # else
    683   /*
    684    * Fallback solution.
    685    */
    686   int rc;
    687 
    688   if (number == 0.0) {
    689     /*
    690      * In IEEE 754 the sign of zero is ignored in comparisons, so we
    691      * have to handle this as a special case by examining the sign bit
    692      * directly.
    693      */
    694 #  if defined(USE_IEEE_754)
    695     *is_negative = trio_is_negative(number);
    696 #  else
    697     *is_negative = TRIO_FALSE; /* FIXME */
    698 #  endif
    699     return TRIO_FP_ZERO;
    700   }
    701   if (trio_isnan(number)) {
    702     *is_negative = TRIO_FALSE;
    703     return TRIO_FP_NAN;
    704   }
    705   if ((rc = trio_isinf(number))) {
    706     *is_negative = (rc == -1);
    707     return TRIO_FP_INFINITE;
    708   }
    709   if ((number > 0.0) && (number < DBL_MIN)) {
    710     *is_negative = TRIO_FALSE;
    711     return TRIO_FP_SUBNORMAL;
    712   }
    713   if ((number < 0.0) && (number > -DBL_MIN)) {
    714     *is_negative = TRIO_TRUE;
    715     return TRIO_FP_SUBNORMAL;
    716   }
    717   *is_negative = (number < 0.0);
    718   return TRIO_FP_NORMAL;
    719 
    720 # endif
    721 #endif
    722 }
    723 
    724 /**
    725    Examine the sign of a number.
    726 
    727    @param number An arbitrary floating-point number.
    728    @return Boolean value indicating whether or not the number has the
    729    sign bit set (i.e. is negative).
    730 */
    731 TRIO_PUBLIC int
    732 trio_signbit
    733 TRIO_ARGS1((number),
    734 	   double number)
    735 {
    736   int is_negative;
    737 
    738   (void)trio_fpclassify_and_signbit(number, &is_negative);
    739   return is_negative;
    740 }
    741 
    742 #if 0
    743 	/* Temporary fix - this routine is not used in libxml */
    744 /**
    745    Examine the class of a number.
    746 
    747    @param number An arbitrary floating-point number.
    748    @return Enumerable value indicating the class of @p number
    749 */
    750 TRIO_PUBLIC int
    751 trio_fpclassify
    752 TRIO_ARGS1((number),
    753 	   double number)
    754 {
    755   int dummy;
    756 
    757   return trio_fpclassify_and_signbit(number, &dummy);
    758 }
    759 
    760 #endif
    761 
    762 /** @} SpecialQuantities */
    763 
    764 /*************************************************************************
    765  * For test purposes.
    766  *
    767  * Add the following compiler option to include this test code.
    768  *
    769  *  Unix : -DSTANDALONE
    770  *  VMS  : /DEFINE=(STANDALONE)
    771  */
    772 #if defined(STANDALONE)
    773 # include <stdio.h>
    774 
    775 static TRIO_CONST char *
    776 getClassification
    777 TRIO_ARGS1((type),
    778 	   int type)
    779 {
    780   switch (type) {
    781   case TRIO_FP_INFINITE:
    782     return "FP_INFINITE";
    783   case TRIO_FP_NAN:
    784     return "FP_NAN";
    785   case TRIO_FP_NORMAL:
    786     return "FP_NORMAL";
    787   case TRIO_FP_SUBNORMAL:
    788     return "FP_SUBNORMAL";
    789   case TRIO_FP_ZERO:
    790     return "FP_ZERO";
    791   default:
    792     return "FP_UNKNOWN";
    793   }
    794 }
    795 
    796 static void
    797 print_class
    798 TRIO_ARGS2((prefix, number),
    799 	   TRIO_CONST char *prefix,
    800 	   double number)
    801 {
    802   printf("%-6s: %s %-15s %g\n",
    803 	 prefix,
    804 	 trio_signbit(number) ? "-" : "+",
    805 	 getClassification(TRIO_FPCLASSIFY(number)),
    806 	 number);
    807 }
    808 
    809 int main(TRIO_NOARGS)
    810 {
    811   double my_nan;
    812   double my_pinf;
    813   double my_ninf;
    814 # if defined(TRIO_PLATFORM_UNIX)
    815   void (*signal_handler) TRIO_PROTO((int));
    816 # endif
    817 
    818   my_nan = trio_nan();
    819   my_pinf = trio_pinf();
    820   my_ninf = trio_ninf();
    821 
    822   print_class("Nan", my_nan);
    823   print_class("PInf", my_pinf);
    824   print_class("NInf", my_ninf);
    825   print_class("PZero", 0.0);
    826   print_class("NZero", -0.0);
    827   print_class("PNorm", 1.0);
    828   print_class("NNorm", -1.0);
    829   print_class("PSub", 1.01e-307 - 1.00e-307);
    830   print_class("NSub", 1.00e-307 - 1.01e-307);
    831 
    832   printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    833 	 my_nan,
    834 	 ((unsigned char *)&my_nan)[0],
    835 	 ((unsigned char *)&my_nan)[1],
    836 	 ((unsigned char *)&my_nan)[2],
    837 	 ((unsigned char *)&my_nan)[3],
    838 	 ((unsigned char *)&my_nan)[4],
    839 	 ((unsigned char *)&my_nan)[5],
    840 	 ((unsigned char *)&my_nan)[6],
    841 	 ((unsigned char *)&my_nan)[7],
    842 	 trio_isnan(my_nan), trio_isinf(my_nan));
    843   printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    844 	 my_pinf,
    845 	 ((unsigned char *)&my_pinf)[0],
    846 	 ((unsigned char *)&my_pinf)[1],
    847 	 ((unsigned char *)&my_pinf)[2],
    848 	 ((unsigned char *)&my_pinf)[3],
    849 	 ((unsigned char *)&my_pinf)[4],
    850 	 ((unsigned char *)&my_pinf)[5],
    851 	 ((unsigned char *)&my_pinf)[6],
    852 	 ((unsigned char *)&my_pinf)[7],
    853 	 trio_isnan(my_pinf), trio_isinf(my_pinf));
    854   printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    855 	 my_ninf,
    856 	 ((unsigned char *)&my_ninf)[0],
    857 	 ((unsigned char *)&my_ninf)[1],
    858 	 ((unsigned char *)&my_ninf)[2],
    859 	 ((unsigned char *)&my_ninf)[3],
    860 	 ((unsigned char *)&my_ninf)[4],
    861 	 ((unsigned char *)&my_ninf)[5],
    862 	 ((unsigned char *)&my_ninf)[6],
    863 	 ((unsigned char *)&my_ninf)[7],
    864 	 trio_isnan(my_ninf), trio_isinf(my_ninf));
    865 
    866 # if defined(TRIO_PLATFORM_UNIX)
    867   signal_handler = signal(SIGFPE, SIG_IGN);
    868 # endif
    869 
    870   my_pinf = DBL_MAX + DBL_MAX;
    871   my_ninf = -my_pinf;
    872   my_nan = my_pinf / my_pinf;
    873 
    874 # if defined(TRIO_PLATFORM_UNIX)
    875   signal(SIGFPE, signal_handler);
    876 # endif
    877 
    878   printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    879 	 my_nan,
    880 	 ((unsigned char *)&my_nan)[0],
    881 	 ((unsigned char *)&my_nan)[1],
    882 	 ((unsigned char *)&my_nan)[2],
    883 	 ((unsigned char *)&my_nan)[3],
    884 	 ((unsigned char *)&my_nan)[4],
    885 	 ((unsigned char *)&my_nan)[5],
    886 	 ((unsigned char *)&my_nan)[6],
    887 	 ((unsigned char *)&my_nan)[7],
    888 	 trio_isnan(my_nan), trio_isinf(my_nan));
    889   printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    890 	 my_pinf,
    891 	 ((unsigned char *)&my_pinf)[0],
    892 	 ((unsigned char *)&my_pinf)[1],
    893 	 ((unsigned char *)&my_pinf)[2],
    894 	 ((unsigned char *)&my_pinf)[3],
    895 	 ((unsigned char *)&my_pinf)[4],
    896 	 ((unsigned char *)&my_pinf)[5],
    897 	 ((unsigned char *)&my_pinf)[6],
    898 	 ((unsigned char *)&my_pinf)[7],
    899 	 trio_isnan(my_pinf), trio_isinf(my_pinf));
    900   printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
    901 	 my_ninf,
    902 	 ((unsigned char *)&my_ninf)[0],
    903 	 ((unsigned char *)&my_ninf)[1],
    904 	 ((unsigned char *)&my_ninf)[2],
    905 	 ((unsigned char *)&my_ninf)[3],
    906 	 ((unsigned char *)&my_ninf)[4],
    907 	 ((unsigned char *)&my_ninf)[5],
    908 	 ((unsigned char *)&my_ninf)[6],
    909 	 ((unsigned char *)&my_ninf)[7],
    910 	 trio_isnan(my_ninf), trio_isinf(my_ninf));
    911 
    912   return 0;
    913 }
    914 #endif
    915