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

   1 /*************************************************************************
   2  *
   3  * $Id: trionan.c,v 1.33 2005/05/29 11:57:25 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  * Include files
  36  */
  37 #include "triodef.h"
  38 #include "trionan.h"
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include <limits.h>
  43 #if !defined(TRIO_PLATFORM_SYMBIAN)
  44 # include <float.h>
  45 #endif
  46 #if defined(TRIO_PLATFORM_UNIX)
  47 # include <signal.h>
  48 #endif
  49 #if defined(TRIO_COMPILER_DECC)
  50 # include <fp_class.h>
  51 #endif
  52 #include <assert.h>
  53
  54 #if defined(TRIO_DOCUMENTATION)
  55 # include "doc/doc_nan.h"
  56 #endif
  57 /** @addtogroup SpecialQuantities
  58     @{
  59 */
  60
  61 /*************************************************************************
  62  * Definitions
  63  */
  64
  65 #if !defined(TRIO_PUBLIC_NAN)
  66 # define TRIO_PUBLIC_NAN TRIO_PUBLIC
  67 #endif
  68 #if !defined(TRIO_PRIVATE_NAN)
  69 # define TRIO_PRIVATE_NAN TRIO_PRIVATE
  70 #endif
  71
  72 #define TRIO_TRUE (1 == 1)
  73 #define TRIO_FALSE (0 == 1)
  74
  75 /*
  76  * We must enable IEEE floating-point on Alpha
  77  */
  78 #if defined(__alpha) && !defined(_IEEE_FP)
  79 # if defined(TRIO_COMPILER_DECC)
  80 #  if defined(TRIO_PLATFORM_VMS)
  81 #   error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"
  82 #  else
  83 #   if !defined(_CFE)
  84 #    error "Must be compiled with option -ieee"
  85 #   endif
  86 #  endif
  87 # else
  88 #  if defined(TRIO_COMPILER_GCC)
  89 #   error "Must be compiled with option -mieee"
  90 #  endif
  91 # endif
  92 #endif /* __alpha && ! _IEEE_FP */
  93
  94 /*
  95  * In ANSI/IEEE 754-1985 64-bits double format numbers have the
  96  * following properties (amoungst others)
  97  *
  98  *   o FLT_RADIX == 2: binary encoding
  99  *   o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used
 100  *     to indicate special numbers (e.g. NaN and Infinity), so the
 101  *     maximum exponent is 10 bits wide (2^10 == 1024).
 102  *   o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because
 103  *     numbers are normalized the initial binary 1 is represented
 104  *     implicitly (the so-called "hidden bit"), which leaves us with
 105  *     the ability to represent 53 bits wide mantissa.
 106  */
 107 #if defined(__STDC_IEC_559__)
 108 # define TRIO_IEEE_754
 109 #else
 110 # if (FLT_RADIX - 0 == 2) && (DBL_MAX_EXP - 0 == 1024) && (DBL_MANT_DIG - 0 == 53)
 111 #  define TRIO_IEEE_754
 112 # endif
 113 #endif
 114
 115 /*
 116  * Determine which fpclassify_and_sign() function to use.
 117  */
 118 #if defined(TRIO_FUNC_FPCLASSIFY_AND_SIGNBIT)
 119 # if defined(PREDEF_STANDARD_C99) && defined(fpclassify)
 120 #  define TRIO_FUNC_C99_FPCLASSIFY_AND_SIGNBIT
 121 # else
 122 #  if defined(TRIO_COMPILER_DECC)
 123 #   define TRIO_FUNC_DECC_FPCLASSIFY_AND_SIGNBIT
 124 #  else
 125 #   if defined(TRIO_COMPILER_VISUALC) || defined(TRIO_COMPILER_BORLAND)
 126 #    define TRIO_FUNC_MS_FPCLASSIFY_AND_SIGNBIT
 127 #   else
 128 #    if defined(TRIO_COMPILER_HP) && defined(FP_PLUS_NORM)
 129 #     define TRIO_FUNC_HP_FPCLASSIFY_AND_SIGNBIT
 130 #    else
 131 #     if defined(TRIO_COMPILER_XLC) && defined(FP_PLUS_NORM)
 132 #      define TRIO_FUNC_XLC_FPCLASSIFY_AND_SIGNBIT
 133 #     else
 134 #      define TRIO_FUNC_INTERNAL_FPCLASSIFY_AND_SIGNBIT
 135 #     endif
 136 #    endif
 137 #   endif
 138 #  endif
 139 # endif
 140 #endif
 141
 142 /*
 143  * Determine how to generate negative zero.
 144  */
 145 #if defined(TRIO_FUNC_NZERO)
 146 # if defined(TRIO_IEEE_754)
 147 #  define TRIO_NZERO_IEEE_754
 148 # else
 149 #  define TRIO_NZERO_FALLBACK
 150 # endif
 151 #endif
 152
 153 /*
 154  * Determine how to generate positive infinity.
 155  */
 156 #if defined(TRIO_FUNC_PINF)
 157 # if defined(INFINITY) && defined(__STDC_IEC_559__)
 158 #  define TRIO_PINF_C99_MACRO
 159 # else
 160 #  if defined(TRIO_IEEE_754)
 161 #   define TRIO_PINF_IEEE_754
 162 #  else
 163 #   define TRIO_PINF_FALLBACK
 164 #  endif
 165 # endif
 166 #endif
 167
 168 /*
 169  * Determine how to generate NaN.
 170  */
 171 #if defined(TRIO_FUNC_NAN)
 172 # if defined(PREDEF_STANDARD_C99) && !defined(TRIO_COMPILER_DECC)
 173 #  define TRIO_NAN_C99_FUNCTION
 174 # else
 175 #  if defined(NAN) && defined(__STDC_IEC_559__)
 176 #   define TRIO_NAN_C99_MACRO
 177 #  else
 178 #   if defined(TRIO_IEEE_754)
 179 #    define TRIO_NAN_IEEE_754
 180 #   else
 181 #    define TRIO_NAN_FALLBACK
 182 #   endif
 183 #  endif
 184 # endif
 185 #endif
 186
 187 /*
 188  * Resolve internal dependencies.
 189  */
 190 #if defined(TRIO_FUNC_INTERNAL_FPCLASSIFY_AND_SIGNBIT)
 191 # define TRIO_FUNC_INTERNAL_ISNAN
 192 # define TRIO_FUNC_INTERNAL_ISINF
 193 # if defined(TRIO_IEEE_754)
 194 #  define TRIO_FUNC_INTERNAL_IS_SPECIAL_QUANTITY
 195 #  define TRIO_FUNC_INTERNAL_IS_NEGATIVE
 196 # endif
 197 #endif
 198
 199 #if defined(TRIO_NZERO_IEEE_754) \
 200  || defined(TRIO_PINF_IEEE_754) \
 201  || defined(TRIO_NAN_IEEE_754)
 202 # define TRIO_FUNC_INTERNAL_MAKE_DOUBLE
 203 #endif
 204
 205 #if defined(TRIO_FUNC_INTERNAL_ISNAN)
 206 # if defined(PREDEF_STANDARD_XPG3)
 207 #  define TRIO_INTERNAL_ISNAN_XPG3
 208 # else
 209 #  if defined(TRIO_IEEE_754)
 210 #   define TRIO_INTERNAL_ISNAN_IEEE_754
 211 #  else
 212 #   define TRIO_INTERNAL_ISNAN_FALLBACK
 213 #  endif
 214 # endif
 215 #endif
 216
 217 #if defined(TRIO_FUNC_INTERNAL_ISINF)
 218 # if defined(TRIO_IEEE_754)
 219 #  define TRIO_INTERNAL_ISINF_IEEE_754
 220 # else
 221 #  define TRIO_INTERNAL_ISINF_FALLBACK
 222 # endif
 223 #endif
 224
 225 /*************************************************************************
 226  * Constants
 227  */
 228
 229 #if !defined(TRIO_EMBED_NAN)
 230 static TRIO_CONST char rcsid[] = "@(#)$Id: trionan.c,v 1.33 2005/05/29 11:57:25 breese Exp $";
 231 #endif
 232
 233 #if defined(TRIO_FUNC_INTERNAL_MAKE_DOUBLE) \
 234  || defined(TRIO_FUNC_INTERNAL_IS_SPECIAL_QUANTITY) \
 235  || defined(TRIO_FUNC_INTERNAL_IS_NEGATIVE)
 236 /*
 237  * Endian-agnostic indexing macro.
 238  *
 239  * The value of internalEndianMagic, when converted into a 64-bit
 240  * integer, becomes 0x0706050403020100 (we could have used a 64-bit
 241  * integer value instead of a double, but not all platforms supports
 242  * that type). The value is automatically encoded with the correct
 243  * endianess by the compiler, which means that we can support any
 244  * kind of endianess. The individual bytes are then used as an index
 245  * for the IEEE 754 bit-patterns and masks.
 246  */
 247 #define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])
 248 static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
 249 #endif
 250
 251 #if defined(TRIO_FUNC_INTERNAL_IS_SPECIAL_QUANTITY)
 252 /* Mask for the exponent */
 253 static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
 254   0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
 255 };
 256
 257 /* Mask for the mantissa */
 258 static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
 259   0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
 260 };
 261 #endif
 262
 263 #if defined(TRIO_FUNC_INTERNAL_IS_NEGATIVE)
 264 /* Mask for the sign bit */
 265 static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
 266   0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
 267 };
 268 #endif
 269
 270 #if defined(TRIO_NZERO_IEEE_754)
 271 /* Bit-pattern for negative zero */
 272 static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
 273   0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
 274 };
 275 #endif
 276
 277 #if defined(TRIO_PINF_IEEE_754)
 278 /* Bit-pattern for infinity */
 279 static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
 280   0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
 281 };
 282 #endif
 283
 284 #if defined(TRIO_NAN_IEEE_754)
 285 /* Bit-pattern for quiet NaN */
 286 static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
 287   0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
 288 };
 289 #endif
 290
 291
 292 /*************************************************************************
 293  * Internal functions
 294  */
 295
 296 /*
 297  * internal_make_double
 298  */
 299 #if defined(TRIO_FUNC_INTERNAL_MAKE_DOUBLE)
 300
 301 TRIO_PRIVATE_NAN double
 302 internal_make_double
 303 TRIO_ARGS1((values),
 304            TRIO_CONST unsigned char *values)
 305 {
 306   TRIO_VOLATILE double result;
 307   int i;
 308
 309   for (i = 0; i < (int)sizeof(double); i++) {
 310     ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
 311   }
 312   return result;
 313 }
 314
 315 #endif
 316
 317 /*
 318  * internal_is_special_quantity
 319  */
 320 #if defined(TRIO_FUNC_INTERNAL_IS_SPECIAL_QUANTITY)
 321
 322 TRIO_PRIVATE_NAN int
 323 internal_is_special_quantity
 324 TRIO_ARGS2((number, has_mantissa),
 325            double number,
 326            int *has_mantissa)
 327 {
 328   unsigned int i;
 329   unsigned char current;
 330   int is_special_quantity = TRIO_TRUE;
 331
 332   *has_mantissa = 0;
 333
 334   for (i = 0; i < (unsigned int)sizeof(double); i++) {
 335     current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
 336     is_special_quantity
 337       &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
 338     *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
 339   }
 340   return is_special_quantity;
 341 }
 342
 343 #endif
 344
 345 /*
 346  * internal_is_negative
 347  */
 348 #if defined(TRIO_FUNC_INTERNAL_IS_NEGATIVE)
 349
 350 TRIO_PRIVATE_NAN int
 351 internal_is_negative
 352 TRIO_ARGS1((number),
 353            double number)
 354 {
 355   unsigned int i;
 356   int is_negative = TRIO_FALSE;
 357
 358   for (i = 0; i < (unsigned int)sizeof(double); i++) {
 359     is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
 360                     & ieee_754_sign_mask[i]);
 361   }
 362   return is_negative;
 363 }
 364
 365 #endif
 366
 367 #if defined(TRIO_FUNC_C99_FPCLASSIFY_AND_SIGNBIT)
 368
 369 TRIO_PRIVATE_NAN TRIO_INLINE int
 370 c99_fpclassify_and_signbit
 371 TRIO_ARGS2((number, is_negative),
 372            double number,
 373            int *is_negative)
 374 {
 375   *is_negative = signbit(number);
 376   switch (fpclassify(number)) {
 377   case FP_NAN:
 378     return TRIO_FP_NAN;
 379   case FP_INFINITE:
 380     return TRIO_FP_INFINITE;
 381   case FP_SUBNORMAL:
 382     return TRIO_FP_SUBNORMAL;
 383   case FP_ZERO:
 384     return TRIO_FP_ZERO;
 385   default:
 386     return TRIO_FP_NORMAL;
 387   }
 388 }
 389
 390 #endif /* TRIO_FUNC_C99_FPCLASSIFY_AND_SIGNBIT */
 391
 392 #if defined(TRIO_FUNC_DECC_FPCLASSIFY_AND_SIGNBIT)
 393
 394 TRIO_PRIVATE_NAN TRIO_INLINE int
 395 decc_fpclassify_and_signbit
 396 TRIO_ARGS2((number, is_negative),
 397           double number,
 398           int *is_negative)
 399 {
 400   switch (fp_class(number)) {
 401   case FP_QNAN:
 402   case FP_SNAN:
 403     *is_negative = TRIO_FALSE; /* NaN has no sign */
 404     return TRIO_FP_NAN;
 405   case FP_POS_INF:
 406     *is_negative = TRIO_FALSE;
 407     return TRIO_FP_INFINITE;
 408   case FP_NEG_INF:
 409     *is_negative = TRIO_TRUE;
 410     return TRIO_FP_INFINITE;
 411   case FP_POS_DENORM:
 412     *is_negative = TRIO_FALSE;
 413     return TRIO_FP_SUBNORMAL;
 414   case FP_NEG_DENORM:
 415     *is_negative = TRIO_TRUE;
 416     return TRIO_FP_SUBNORMAL;
 417   case FP_POS_ZERO:
 418     *is_negative = TRIO_FALSE;
 419     return TRIO_FP_ZERO;
 420   case FP_NEG_ZERO:
 421     *is_negative = TRIO_TRUE;
 422     return TRIO_FP_ZERO;
 423   case FP_POS_NORM:
 424     *is_negative = TRIO_FALSE;
 425     return TRIO_FP_NORMAL;
 426   case FP_NEG_NORM:
 427     *is_negative = TRIO_TRUE;
 428     return TRIO_FP_NORMAL;
 429   default:
 430     *is_negative = (number < 0.0);
 431     return TRIO_FP_NORMAL;
 432   }
 433 }
 434
 435 #endif /* TRIO_FUNC_DECC_FPCLASSIFY_AND_SIGNBIT */
 436
 437 #if defined(TRIO_FUNC_MS_FPCLASSIFY_AND_SIGNBIT)
 438
 439 TRIO_PRIVATE_NAN int
 440 ms_fpclassify_and_signbit
 441 TRIO_ARGS2((number, is_negative),
 442           double number,
 443           int *is_negative)
 444 {
 445   int result;
 446 # if defined(TRIO_COMPILER_BORLAND)
 447   /*
 448    * The floating-point precision may be changed by the Borland _fpclass()
 449    * function, so we have to save and restore the floating-point control mask.
 450    */
 451   unsigned int mask;
 452   /* Remember the old mask */
 453   mask = _control87(0, 0);
 454 # endif
 455
 456   switch (_fpclass(number)) {
 457   case _FPCLASS_QNAN:
 458   case _FPCLASS_SNAN:
 459     *is_negative = TRIO_FALSE; /* NaN has no sign */
 460     result = TRIO_FP_NAN;
 461     break;
 462   case _FPCLASS_PINF:
 463     *is_negative = TRIO_FALSE;
 464     result = TRIO_FP_INFINITE;
 465     break;
 466   case _FPCLASS_NINF:
 467     *is_negative = TRIO_TRUE;
 468     result = TRIO_FP_INFINITE;
 469     break;
 470   case _FPCLASS_PD:
 471     *is_negative = TRIO_FALSE;
 472     result = TRIO_FP_SUBNORMAL;
 473     break;
 474   case _FPCLASS_ND:
 475     *is_negative = TRIO_TRUE;
 476     result = TRIO_FP_SUBNORMAL;
 477     break;
 478   case _FPCLASS_PZ:
 479     *is_negative = TRIO_FALSE;
 480     result = TRIO_FP_ZERO;
 481     break;
 482   case _FPCLASS_NZ:
 483     *is_negative = TRIO_TRUE;
 484     result = TRIO_FP_ZERO;
 485     break;
 486   case _FPCLASS_PN:
 487     *is_negative = TRIO_FALSE;
 488     result = TRIO_FP_NORMAL;
 489     break;
 490   case _FPCLASS_NN:
 491     *is_negative = TRIO_TRUE;
 492     result = TRIO_FP_NORMAL;
 493     break;
 494   default:
 495     *is_negative = (number < 0.0);
 496     result = TRIO_FP_NORMAL;
 497     break;
 498   }
 499
 500 # if defined(TRIO_COMPILER_BORLAND)
 501   /* Restore the old precision */
 502   (void)_control87(mask, MCW_PC);
 503 # endif
 504
 505   return result;
 506 }
 507
 508 #endif /* TRIO_FUNC_MS_FPCLASSIFY_AND_SIGNBIT */
 509
 510 #if defined(TRIO_FUNC_HP_FPCLASSIFY_AND_SIGNBIT)
 511
 512 TRIO_PRIVATE_NAN TRIO_INLINE int
 513 hp_fpclassify_and_signbit
 514 TRIO_ARGS2((number, is_negative),
 515           double number,
 516           int *is_negative)
 517 {
 518   /*
 519    * HP-UX 9.x and 10.x have an fpclassify() function, that is different
 520    * from the C99 fpclassify() macro supported on HP-UX 11.x.
 521    */
 522   switch (fpclassify(number)) {
 523   case FP_QNAN:
 524   case FP_SNAN:
 525     *is_negative = TRIO_FALSE; /* NaN has no sign */
 526     return TRIO_FP_NAN;
 527   case FP_PLUS_INF:
 528     *is_negative = TRIO_FALSE;
 529     return TRIO_FP_INFINITE;
 530   case FP_MINUS_INF:
 531     *is_negative = TRIO_TRUE;
 532     return TRIO_FP_INFINITE;
 533   case FP_PLUS_DENORM:
 534     *is_negative = TRIO_FALSE;
 535     return TRIO_FP_SUBNORMAL;
 536   case FP_MINUS_DENORM:
 537     *is_negative = TRIO_TRUE;
 538     return TRIO_FP_SUBNORMAL;
 539   case FP_PLUS_ZERO:
 540     *is_negative = TRIO_FALSE;
 541     return TRIO_FP_ZERO;
 542   case FP_MINUS_ZERO:
 543     *is_negative = TRIO_TRUE;
 544     return TRIO_FP_ZERO;
 545   case FP_PLUS_NORM:
 546     *is_negative = TRIO_FALSE;
 547     return TRIO_FP_NORMAL;
 548   case FP_MINUS_NORM:
 549     *is_negative = TRIO_TRUE;
 550     return TRIO_FP_NORMAL;
 551   default:
 552     *is_negative = (number < 0.0);
 553     return TRIO_FP_NORMAL;
 554   }
 555 }
 556
 557 #endif /* TRIO_FUNC_HP_FPCLASSIFY_AND_SIGNBIT */
 558
 559 #if defined(TRIO_FUNC_XLC_FPCLASSIFY_AND_SIGNBIT)
 560
 561 TRIO_PRIVATE_NAN TRIO_INLINE int
 562 xlc_fpclassify_and_signbit
 563 TRIO_ARGS2((number, is_negative),
 564           double number,
 565           int *is_negative)
 566 {
 567   /*
 568    * AIX has class() for C, and _class() for C++
 569    */
 570 # if defined(__cplusplus)
 571 #  define AIX_CLASS(n) _class(n)
 572 # else
 573 #  define AIX_CLASS(n) class(n)
 574 # endif
 575
 576   switch (AIX_CLASS(number)) {
 577   case FP_QNAN:
 578   case FP_SNAN:
 579     *is_negative = TRIO_FALSE; /* NaN has no sign */
 580     return TRIO_FP_NAN;
 581   case FP_PLUS_INF:
 582     *is_negative = TRIO_FALSE;
 583     return TRIO_FP_INFINITE;
 584   case FP_MINUS_INF:
 585     *is_negative = TRIO_TRUE;
 586     return TRIO_FP_INFINITE;
 587   case FP_PLUS_DENORM:
 588     *is_negative = TRIO_FALSE;
 589     return TRIO_FP_SUBNORMAL;
 590   case FP_MINUS_DENORM:
 591     *is_negative = TRIO_TRUE;
 592     return TRIO_FP_SUBNORMAL;
 593   case FP_PLUS_ZERO:
 594     *is_negative = TRIO_FALSE;
 595     return TRIO_FP_ZERO;
 596   case FP_MINUS_ZERO:
 597     *is_negative = TRIO_TRUE;
 598     return TRIO_FP_ZERO;
 599   case FP_PLUS_NORM:
 600     *is_negative = TRIO_FALSE;
 601     return TRIO_FP_NORMAL;
 602   case FP_MINUS_NORM:
 603     *is_negative = TRIO_TRUE;
 604     return TRIO_FP_NORMAL;
 605   default:
 606     *is_negative = (number < 0.0);
 607     return TRIO_FP_NORMAL;
 608   }
 609 }
 610
 611 #endif /* TRIO_FUNC_XLC_FPCLASSIFY_AND_SIGNBIT */
 612
 613 #if defined(TRIO_FUNC_INTERNAL_ISNAN)
 614
 615 TRIO_PRIVATE_NAN TRIO_INLINE int
 616 internal_isnan
 617 TRIO_ARGS1((number),
 618            double number)
 619 {
 620 # if defined(TRIO_INTERNAL_ISNAN_XPG3) || defined(TRIO_PLATFORM_SYMBIAN)
 621   /*
 622    * XPG3 defines isnan() as a function.
 623    */
 624   return isnan(number);
 625
 626 # endif
 627
 628 # if defined(TRIO_INTERNAL_ISNAN_IEEE_754)
 629
 630   /*
 631    * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
 632    * pattern, and a non-empty mantissa.
 633    */
 634   int has_mantissa;
 635   int is_special_quantity;
 636
 637   is_special_quantity = internal_is_special_quantity(number, &has_mantissa);
 638
 639   return (is_special_quantity && has_mantissa);
 640
 641 # endif
 642
 643 # if defined(TRIO_INTERNAL_ISNAN_FALLBACK)
 644
 645   /*
 646    * Fallback solution
 647    */
 648   int status;
 649   double integral, fraction;
 650
 651 #  if defined(TRIO_PLATFORM_UNIX)
 652   void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
 653 #  endif
 654
 655   status = (/*
 656              * NaN is the only number which does not compare to itself
 657              */
 658             ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
 659             /*
 660              * Fallback solution if NaN compares to NaN
 661              */
 662             ((number != 0.0) &&
 663              (fraction = modf(number, &integral),
 664               integral == fraction)));
 665
 666 #  if defined(TRIO_PLATFORM_UNIX)
 667   signal(SIGFPE, signal_handler);
 668 #  endif
 669
 670   return status;
 671
 672 # endif
 673 }
 674
 675 #endif /* TRIO_FUNC_INTERNAL_ISNAN */
 676
 677 #if defined(TRIO_FUNC_INTERNAL_ISINF)
 678
 679 TRIO_PRIVATE_NAN TRIO_INLINE int
 680 internal_isinf
 681 TRIO_ARGS1((number),
 682            double number)
 683 {
 684 # if defined(TRIO_PLATFORM_SYMBIAN)
 685
 686   return isinf(number);
 687
 688 # endif
 689
 690 # if defined(TRIO_INTERNAL_ISINF_IEEE_754)
 691   /*
 692    * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
 693    * pattern, and an empty mantissa.
 694    */
 695   int has_mantissa;
 696   int is_special_quantity;
 697
 698   is_special_quantity = internal_is_special_quantity(number, &has_mantissa);
 699
 700   return (is_special_quantity && !has_mantissa)
 701     ? ((number < 0.0) ? -1 : 1)
 702     : 0;
 703
 704 # endif
 705
 706 # if defined(TRIO_INTERNAL_ISINF_FALLBACK)
 707
 708   /*
 709    * Fallback solution.
 710    */
 711   int status;
 712
 713 #  if defined(TRIO_PLATFORM_UNIX)
 714   void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
 715 #  endif
 716
 717   double infinity = trio_pinf();
 718
 719   status = ((number == infinity)
 720             ? 1
 721             : ((number == -infinity) ? -1 : 0));
 722
 723 #  if defined(TRIO_PLATFORM_UNIX)
 724   signal(SIGFPE, signal_handler);
 725 #  endif
 726
 727   return status;
 728
 729 # endif
 730 }
 731
 732 #endif /* TRIO_FUNC_INTERNAL_ISINF */
 733
 734 /*************************************************************************
 735  * Public functions
 736  */
 737
 738 #if defined(TRIO_FUNC_FPCLASSIFY_AND_SIGNBIT)
 739
 740 TRIO_PUBLIC_NAN int
 741 trio_fpclassify_and_signbit
 742 TRIO_ARGS2((number, is_negative),
 743            double number,
 744            int *is_negative)
 745 {
 746   /* The TRIO_FUNC_xxx_FPCLASSIFY_AND_SIGNBIT macros are mutually exclusive */
 747
 748 #if defined(TRIO_FUNC_C99_FPCLASSIFY_AND_SIGNBIT)
 749
 750   return c99_fpclassify_and_signbit(number, is_negative);
 751
 752 #endif
 753
 754 #if defined(TRIO_FUNC_DECC_FPCLASSIFY_AND_SIGNBIT)
 755
 756   return decc_fpclassify_and_signbit(number, is_negative);
 757
 758 #endif
 759
 760 #if defined(TRIO_FUNC_MS_FPCLASSIFY_AND_SIGNBIT)
 761
 762   return ms_fpclassify_and_signbit(number, is_negative);
 763
 764 #endif
 765
 766 #if defined(TRIO_FUNC_HP_FPCLASSIFY_AND_SIGNBIT)
 767
 768   return hp_fpclassify_and_signbit(number, is_negative);
 769
 770 #endif
 771
 772 #if defined(TRIO_FUNC_XLC_FPCLASSIFY_AND_SIGNBIT)
 773
 774   return xlc_fpclassify_and_signbit(number, is_negative);
 775
 776 #endif
 777
 778 #if defined(TRIO_FUNC_INTERNAL_FPCLASSIFY_AND_SIGNBIT)
 779
 780   /*
 781    * Fallback solution.
 782    */
 783   int rc;
 784
 785   if (number == 0.0) {
 786     /*
 787      * In IEEE 754 the sign of zero is ignored in comparisons, so we
 788      * have to handle this as a special case by examining the sign bit
 789      * directly.
 790      */
 791 # if defined(TRIO_IEEE_754)
 792     *is_negative = internal_is_negative(number);
 793 # else
 794     *is_negative = TRIO_FALSE; /* FIXME */
 795 # endif
 796     return TRIO_FP_ZERO;
 797   }
 798   if (internal_isnan(number)) {
 799     *is_negative = TRIO_FALSE;
 800     return TRIO_FP_NAN;
 801   }
 802   rc = internal_isinf(number);
 803   if (rc != 0) {
 804     *is_negative = (rc == -1);
 805     return TRIO_FP_INFINITE;
 806   }
 807   if ((number > 0.0) && (number < DBL_MIN)) {
 808     *is_negative = TRIO_FALSE;
 809     return TRIO_FP_SUBNORMAL;
 810   }
 811   if ((number < 0.0) && (number > -DBL_MIN)) {
 812     *is_negative = TRIO_TRUE;
 813     return TRIO_FP_SUBNORMAL;
 814   }
 815   *is_negative = (number < 0.0);
 816   return TRIO_FP_NORMAL;
 817
 818 #endif
 819 }
 820
 821 #endif
 822
 823 /**
 824    Check for NaN.
 825
 826    @param number An arbitrary floating-point number.
 827    @return Boolean value indicating whether or not the number is a NaN.
 828 */
 829 #if defined(TRIO_FUNC_ISNAN)
 830
 831 TRIO_PUBLIC_NAN int
 832 trio_isnan
 833 TRIO_ARGS1((number),
 834            double number)
 835 {
 836   int dummy;
 837
 838   return (trio_fpclassify_and_signbit(number, &dummy) == TRIO_FP_NAN);
 839 }
 840
 841 #endif
 842
 843 /**
 844    Check for infinity.
 845
 846    @param number An arbitrary floating-point number.
 847    @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.
 848 */
 849 #if defined(TRIO_FUNC_ISINF)
 850
 851 TRIO_PUBLIC_NAN int
 852 trio_isinf
 853 TRIO_ARGS1((number),
 854            double number)
 855 {
 856   int is_negative;
 857
 858   if (trio_fpclassify_and_signbit(number, &is_negative) == TRIO_FP_INFINITE)
 859     {
 860       return (is_negative) ? -1 : 1;
 861     }
 862   else
 863     {
 864       return 0;
 865     }
 866 }
 867
 868 #endif
 869
 870 /**
 871    Check for finity.
 872
 873    @param number An arbitrary floating-point number.
 874    @return Boolean value indicating whether or not the number is a finite.
 875 */
 876 #if defined(TRIO_FUNC_ISFINITE)
 877
 878 TRIO_PUBLIC_NAN int
 879 trio_isfinite
 880 TRIO_ARGS1((number),
 881            double number)
 882 {
 883   int dummy;
 884
 885   switch (trio_fpclassify_and_signbit(number, &dummy))
 886     {
 887     case TRIO_FP_INFINITE:
 888     case TRIO_FP_NAN:
 889       return 0;
 890     default:
 891       return 1;
 892     }
 893 }
 894
 895 #endif
 896
 897 /**
 898    Examine the sign of a number.
 899
 900    @param number An arbitrary floating-point number.
 901    @return Boolean value indicating whether or not the number has the
 902    sign bit set (i.e. is negative).
 903 */
 904 #if defined(TRIO_FUNC_SIGNBIT)
 905
 906 TRIO_PUBLIC_NAN int
 907 trio_signbit
 908 TRIO_ARGS1((number),
 909            double number)
 910 {
 911   int is_negative;
 912
 913   (void)trio_fpclassify_and_signbit(number, &is_negative);
 914   return is_negative;
 915 }
 916
 917 #endif
 918
 919 /**
 920    Examine the class of a number.
 921
 922    @param number An arbitrary floating-point number.
 923    @return Enumerable value indicating the class of @p number
 924 */
 925 #if defined(TRIO_FUNC_FPCLASSIFY)
 926
 927 TRIO_PUBLIC_NAN int
 928 trio_fpclassify
 929 TRIO_ARGS1((number),
 930            double number)
 931 {
 932   int dummy;
 933
 934   return trio_fpclassify_and_signbit(number, &dummy);
 935 }
 936
 937 #endif
 938
 939 /**
 940    Generate negative zero.
 941
 942    @return Floating-point representation of negative zero.
 943 */
 944 #if defined(TRIO_FUNC_NZERO)
 945
 946 TRIO_PUBLIC_NAN double
 947 trio_nzero(TRIO_NOARGS)
 948 {
 949 # if defined(TRIO_NZERO_IEEE_754)
 950
 951   return internal_make_double(ieee_754_negzero_array);
 952
 953 # endif
 954
 955 # if defined(TRIO_NZERO_FALLBACK)
 956
 957   TRIO_VOLATILE double zero = 0.0;
 958
 959   return -zero;
 960
 961 # endif
 962 }
 963
 964 #endif
 965
 966 /**
 967    Generate positive infinity.
 968
 969    @return Floating-point representation of positive infinity.
 970 */
 971 #if defined(TRIO_FUNC_PINF)
 972
 973 TRIO_PUBLIC_NAN double
 974 trio_pinf(TRIO_NOARGS)
 975 {
 976   /* Cache the result */
 977   static double pinf_value = 0.0;
 978
 979   if (pinf_value == 0.0) {
 980
 981 # if defined(TRIO_PINF_C99_MACRO)
 982
 983     pinf_value = (double)INFINITY;
 984
 985 # endif
 986
 987 # if defined(TRIO_PINF_IEEE_754)
 988
 989     pinf_value = internal_make_double(ieee_754_infinity_array);
 990
 991 # endif
 992
 993 # if defined(TRIO_PINF_FALLBACK)
 994     /*
 995      * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
 996      * as infinity. Otherwise we have to resort to an overflow
 997      * operation to generate infinity.
 998      */
 999 #  if defined(TRIO_PLATFORM_UNIX)
1000     void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
1001 #  endif
1002
1003     pinf_value = HUGE_VAL;
1004     if (HUGE_VAL == DBL_MAX) {
1005       /* Force overflow */
1006       pinf_value += HUGE_VAL;
1007     }
1008
1009 #  if defined(TRIO_PLATFORM_UNIX)
1010     signal(SIGFPE, signal_handler);
1011 #  endif
1012
1013 # endif
1014   }
1015   return pinf_value;
1016 }
1017
1018 #endif
1019
1020 /**
1021    Generate negative infinity.
1022
1023    @return Floating-point value of negative infinity.
1024 */
1025 #if defined(TRIO_FUNC_NINF)
1026
1027 TRIO_PUBLIC_NAN double
1028 trio_ninf(TRIO_NOARGS)
1029 {
1030   static double ninf_value = 0.0;
1031
1032   if (ninf_value == 0.0) {
1033     /*
1034      * Negative infinity is calculated by negating positive infinity,
1035      * which can be done because it is legal to do calculations on
1036      * infinity (for example,  1 / infinity == 0).
1037      */
1038     ninf_value = -trio_pinf();
1039   }
1040   return ninf_value;
1041 }
1042
1043 #endif
1044
1045 /**
1046    Generate NaN.
1047
1048    @return Floating-point representation of NaN.
1049 */
1050 #if defined(TRIO_FUNC_NAN)
1051
1052 TRIO_PUBLIC_NAN double
1053 trio_nan(TRIO_NOARGS)
1054 {
1055   /* Cache the result */
1056   static double nan_value = 0.0;
1057
1058   if (nan_value == 0.0) {
1059
1060 # if defined(TRIO_NAN_C99_FUNCTION) || defined(TRIO_PLATFORM_SYMBIAN)
1061
1062     nan_value = nan("");
1063
1064 # endif
1065
1066 # if defined(TRIO_NAN_C99_MACRO)
1067
1068     nan_value = (double)NAN;
1069
1070 # endif
1071
1072 # if defined(TRIO_NAN_IEEE_754)
1073
1074     nan_value = internal_make_double(ieee_754_qnan_array);
1075
1076 # endif
1077
1078 # if defined(TRIO_NAN_FALLBACK)
1079     /*
1080      * There are several ways to generate NaN. The one used here is
1081      * to divide infinity by infinity. I would have preferred to add
1082      * negative infinity to positive infinity, but that yields wrong
1083      * result (infinity) on FreeBSD.
1084      *
1085      * This may fail if the hardware does not support NaN, or if
1086      * the Invalid Operation floating-point exception is unmasked.
1087      */
1088 #  if defined(TRIO_PLATFORM_UNIX)
1089     void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
1090 #  endif
1091
1092     nan_value = trio_pinf() / trio_pinf();
1093
1094 #  if defined(TRIO_PLATFORM_UNIX)
1095     signal(SIGFPE, signal_handler);
1096 #  endif
1097
1098 # endif
1099   }
1100   return nan_value;
1101 }
1102
1103 #endif
1104
1105 /** @} SpecialQuantities */
1106
1107 /*************************************************************************
1108  * For test purposes.
1109  *
1110  * Add the following compiler option to include this test code.
1111  *
1112  *  Unix : -DSTANDALONE
1113  *  VMS  : /DEFINE=(STANDALONE)
1114  */
1115 #if defined(STANDALONE)
1116 # include <stdio.h>
1117
1118 static TRIO_CONST char *
1119 getClassification
1120 TRIO_ARGS1((type),
1121            int type)
1122 {
1123   switch (type) {
1124   case TRIO_FP_INFINITE:
1125     return "FP_INFINITE";
1126   case TRIO_FP_NAN:
1127     return "FP_NAN";
1128   case TRIO_FP_NORMAL:
1129     return "FP_NORMAL";
1130   case TRIO_FP_SUBNORMAL:
1131     return "FP_SUBNORMAL";
1132   case TRIO_FP_ZERO:
1133     return "FP_ZERO";
1134   default:
1135     return "FP_UNKNOWN";
1136   }
1137 }
1138
1139 static void
1140 print_class
1141 TRIO_ARGS2((prefix, number),
1142            TRIO_CONST char *prefix,
1143            double number)
1144 {
1145   printf("%-6s: %s %-15s %g\n",
1146          prefix,
1147          trio_signbit(number) ? "-" : "+",
1148          getClassification(trio_fpclassify(number)),
1149          number);
1150 }
1151
1152 int main(TRIO_NOARGS)
1153 {
1154   double my_nan;
1155   double my_pinf;
1156   double my_ninf;
1157 # if defined(TRIO_PLATFORM_UNIX)
1158   void (*signal_handler) TRIO_PROTO((int));
1159 # endif
1160
1161   my_nan = trio_nan();
1162   my_pinf = trio_pinf();
1163   my_ninf = trio_ninf();
1164
1165   print_class("Nan", my_nan);
1166   print_class("PInf", my_pinf);
1167   print_class("NInf", my_ninf);
1168   print_class("PZero", 0.0);
1169   print_class("NZero", -0.0);
1170   print_class("PNorm", 1.0);
1171   print_class("NNorm", -1.0);
1172   print_class("PSub", 1.01e-307 - 1.00e-307);
1173   print_class("NSub", 1.00e-307 - 1.01e-307);
1174
1175   printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1176          my_nan,
1177          ((unsigned char *)&my_nan)[0],
1178          ((unsigned char *)&my_nan)[1],
1179          ((unsigned char *)&my_nan)[2],
1180          ((unsigned char *)&my_nan)[3],
1181          ((unsigned char *)&my_nan)[4],
1182          ((unsigned char *)&my_nan)[5],
1183          ((unsigned char *)&my_nan)[6],
1184          ((unsigned char *)&my_nan)[7],
1185          trio_isnan(my_nan), trio_isinf(my_nan), trio_isfinite(my_nan));
1186   printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1187          my_pinf,
1188          ((unsigned char *)&my_pinf)[0],
1189          ((unsigned char *)&my_pinf)[1],
1190          ((unsigned char *)&my_pinf)[2],
1191          ((unsigned char *)&my_pinf)[3],
1192          ((unsigned char *)&my_pinf)[4],
1193          ((unsigned char *)&my_pinf)[5],
1194          ((unsigned char *)&my_pinf)[6],
1195          ((unsigned char *)&my_pinf)[7],
1196          trio_isnan(my_pinf), trio_isinf(my_pinf), trio_isfinite(my_pinf));
1197   printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1198          my_ninf,
1199          ((unsigned char *)&my_ninf)[0],
1200          ((unsigned char *)&my_ninf)[1],
1201          ((unsigned char *)&my_ninf)[2],
1202          ((unsigned char *)&my_ninf)[3],
1203          ((unsigned char *)&my_ninf)[4],
1204          ((unsigned char *)&my_ninf)[5],
1205          ((unsigned char *)&my_ninf)[6],
1206          ((unsigned char *)&my_ninf)[7],
1207          trio_isnan(my_ninf), trio_isinf(my_ninf), trio_isfinite(my_ninf));
1208
1209 # if defined(TRIO_PLATFORM_UNIX)
1210   signal_handler = signal(SIGFPE, SIG_IGN);
1211 # endif
1212
1213   my_pinf = DBL_MAX + DBL_MAX;
1214   my_ninf = -my_pinf;
1215   my_nan = my_pinf / my_pinf;
1216
1217 # if defined(TRIO_PLATFORM_UNIX)
1218   signal(SIGFPE, signal_handler);
1219 # endif
1220
1221   printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1222          my_nan,
1223          ((unsigned char *)&my_nan)[0],
1224          ((unsigned char *)&my_nan)[1],
1225          ((unsigned char *)&my_nan)[2],
1226          ((unsigned char *)&my_nan)[3],
1227          ((unsigned char *)&my_nan)[4],
1228          ((unsigned char *)&my_nan)[5],
1229          ((unsigned char *)&my_nan)[6],
1230          ((unsigned char *)&my_nan)[7],
1231          trio_isnan(my_nan), trio_isinf(my_nan), trio_isfinite(my_nan));
1232   printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1233          my_pinf,
1234          ((unsigned char *)&my_pinf)[0],
1235          ((unsigned char *)&my_pinf)[1],
1236          ((unsigned char *)&my_pinf)[2],
1237          ((unsigned char *)&my_pinf)[3],
1238          ((unsigned char *)&my_pinf)[4],
1239          ((unsigned char *)&my_pinf)[5],
1240          ((unsigned char *)&my_pinf)[6],
1241          ((unsigned char *)&my_pinf)[7],
1242          trio_isnan(my_pinf), trio_isinf(my_pinf), trio_isfinite(my_pinf));
1243   printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d, %2d)\n",
1244          my_ninf,
1245          ((unsigned char *)&my_ninf)[0],
1246          ((unsigned char *)&my_ninf)[1],
1247          ((unsigned char *)&my_ninf)[2],
1248          ((unsigned char *)&my_ninf)[3],
1249          ((unsigned char *)&my_ninf)[4],
1250          ((unsigned char *)&my_ninf)[5],
1251          ((unsigned char *)&my_ninf)[6],
1252          ((unsigned char *)&my_ninf)[7],
1253          trio_isnan(my_ninf), trio_isinf(my_ninf), trio_isfinite(my_ninf));
1254
1255   return 0;
1256 }
1257 #endif