Pileus Git - ~andy/fetchmail/blob - trio/trionan.c

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