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