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