1/*************************************************************************
2 *
3 * $Id: trionan.c 3790 2008-09-01 13:08:57Z veillard $
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
115static TRIO_CONST char rcsid[] = "@(#)$Id: trionan.c 3790 2008-09-01 13:08:57Z veillard $";
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)
133static TRIO_CONST double internalEndianMagic = 7.949928895127362e-275;
134#else
135static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
136#endif
137
138/* Mask for the exponent */
139static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
140  0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
141};
142
143/* Mask for the mantissa */
144static 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 */
149static 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 */
154static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
155  0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
156};
157
158/* Bit-pattern for infinity */
159static 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 */
164static 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 */
176TRIO_PRIVATE double
177trio_make_double
178TRIO_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 */
193TRIO_PRIVATE int
194trio_is_special_quantity
195TRIO_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 */
217TRIO_PRIVATE int
218trio_is_negative
219TRIO_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*/
240TRIO_PUBLIC double
241trio_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*/
257TRIO_PUBLIC double
258trio_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*/
301TRIO_PUBLIC double
302trio_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*/
322TRIO_PUBLIC double
323trio_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*/
370TRIO_PUBLIC int
371trio_isnan
372TRIO_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*/
441TRIO_PUBLIC int
442trio_isinf
443TRIO_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*/
519TRIO_PUBLIC int
520trio_isfinite
521TRIO_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 */
559TRIO_PUBLIC int
560trio_fpclassify_and_signbit
561TRIO_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*/
731TRIO_PUBLIC int
732trio_signbit
733TRIO_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*/
750TRIO_PUBLIC int
751trio_fpclassify
752TRIO_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
775static TRIO_CONST char *
776getClassification
777TRIO_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
796static void
797print_class
798TRIO_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
809int main(TRIO_NOARGS)
810{
811  double my_nan;
812  double my_pinf;
813  double my_ninf;
814# if defined(TRIO_PLATFORM_UNIX)
815  void (*signal_handler) TRIO_PROTO((int));
816# endif
817
818  my_nan = trio_nan();
819  my_pinf = trio_pinf();
820  my_ninf = trio_ninf();
821
822  print_class("Nan", my_nan);
823  print_class("PInf", my_pinf);
824  print_class("NInf", my_ninf);
825  print_class("PZero", 0.0);
826  print_class("NZero", -0.0);
827  print_class("PNorm", 1.0);
828  print_class("NNorm", -1.0);
829  print_class("PSub", 1.01e-307 - 1.00e-307);
830  print_class("NSub", 1.00e-307 - 1.01e-307);
831
832  printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
833	 my_nan,
834	 ((unsigned char *)&my_nan)[0],
835	 ((unsigned char *)&my_nan)[1],
836	 ((unsigned char *)&my_nan)[2],
837	 ((unsigned char *)&my_nan)[3],
838	 ((unsigned char *)&my_nan)[4],
839	 ((unsigned char *)&my_nan)[5],
840	 ((unsigned char *)&my_nan)[6],
841	 ((unsigned char *)&my_nan)[7],
842	 trio_isnan(my_nan), trio_isinf(my_nan));
843  printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
844	 my_pinf,
845	 ((unsigned char *)&my_pinf)[0],
846	 ((unsigned char *)&my_pinf)[1],
847	 ((unsigned char *)&my_pinf)[2],
848	 ((unsigned char *)&my_pinf)[3],
849	 ((unsigned char *)&my_pinf)[4],
850	 ((unsigned char *)&my_pinf)[5],
851	 ((unsigned char *)&my_pinf)[6],
852	 ((unsigned char *)&my_pinf)[7],
853	 trio_isnan(my_pinf), trio_isinf(my_pinf));
854  printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
855	 my_ninf,
856	 ((unsigned char *)&my_ninf)[0],
857	 ((unsigned char *)&my_ninf)[1],
858	 ((unsigned char *)&my_ninf)[2],
859	 ((unsigned char *)&my_ninf)[3],
860	 ((unsigned char *)&my_ninf)[4],
861	 ((unsigned char *)&my_ninf)[5],
862	 ((unsigned char *)&my_ninf)[6],
863	 ((unsigned char *)&my_ninf)[7],
864	 trio_isnan(my_ninf), trio_isinf(my_ninf));
865
866# if defined(TRIO_PLATFORM_UNIX)
867  signal_handler = signal(SIGFPE, SIG_IGN);
868# endif
869
870  my_pinf = DBL_MAX + DBL_MAX;
871  my_ninf = -my_pinf;
872  my_nan = my_pinf / my_pinf;
873
874# if defined(TRIO_PLATFORM_UNIX)
875  signal(SIGFPE, signal_handler);
876# endif
877
878  printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
879	 my_nan,
880	 ((unsigned char *)&my_nan)[0],
881	 ((unsigned char *)&my_nan)[1],
882	 ((unsigned char *)&my_nan)[2],
883	 ((unsigned char *)&my_nan)[3],
884	 ((unsigned char *)&my_nan)[4],
885	 ((unsigned char *)&my_nan)[5],
886	 ((unsigned char *)&my_nan)[6],
887	 ((unsigned char *)&my_nan)[7],
888	 trio_isnan(my_nan), trio_isinf(my_nan));
889  printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
890	 my_pinf,
891	 ((unsigned char *)&my_pinf)[0],
892	 ((unsigned char *)&my_pinf)[1],
893	 ((unsigned char *)&my_pinf)[2],
894	 ((unsigned char *)&my_pinf)[3],
895	 ((unsigned char *)&my_pinf)[4],
896	 ((unsigned char *)&my_pinf)[5],
897	 ((unsigned char *)&my_pinf)[6],
898	 ((unsigned char *)&my_pinf)[7],
899	 trio_isnan(my_pinf), trio_isinf(my_pinf));
900  printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
901	 my_ninf,
902	 ((unsigned char *)&my_ninf)[0],
903	 ((unsigned char *)&my_ninf)[1],
904	 ((unsigned char *)&my_ninf)[2],
905	 ((unsigned char *)&my_ninf)[3],
906	 ((unsigned char *)&my_ninf)[4],
907	 ((unsigned char *)&my_ninf)[5],
908	 ((unsigned char *)&my_ninf)[6],
909	 ((unsigned char *)&my_ninf)[7],
910	 trio_isnan(my_ninf), trio_isinf(my_ninf));
911
912  return 0;
913}
914#endif
915