trionan.c

   1 /*************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Copyright (C) 2001 Bjorn Reese <breese@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