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  * MERCHANTABILITY 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 (amongst 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$";
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  * endianness by the compiler, which means that we can support any
127  * kind of endianness. 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)
133 static TRIO_CONST double internalEndianMagic = 7.949928895127362e-275;
134 #else
135 static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
136 #endif
137 
138 /* Mask for the exponent */
139 static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
140   0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
141 };
142 
143 /* Mask for the mantissa */
144 static 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 */
149 static 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 */
154 static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
155   0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
156 };
157 
158 /* Bit-pattern for infinity */
159 static 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 */
164 static 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  */
176 TRIO_PRIVATE double
177 trio_make_double
178 TRIO_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  */
193 TRIO_PRIVATE int
194 trio_is_special_quantity
195 TRIO_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  */
217 TRIO_PRIVATE int
218 trio_is_negative
219 TRIO_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 */
240 TRIO_PUBLIC double
trio_nzero(TRIO_NOARGS)241 trio_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 */
257 TRIO_PUBLIC double
trio_pinf(TRIO_NOARGS)258 trio_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 */
301 TRIO_PUBLIC double
trio_ninf(TRIO_NOARGS)302 trio_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 */
322 TRIO_PUBLIC double
trio_nan(TRIO_NOARGS)323 trio_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 */
370 TRIO_PUBLIC int
371 trio_isnan
372 TRIO_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 */
441 TRIO_PUBLIC int
442 trio_isinf
443 TRIO_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 */
519 TRIO_PUBLIC int
520 trio_isfinite
521 TRIO_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  */
559 TRIO_PUBLIC int
560 trio_fpclassify_and_signbit
561 TRIO_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 */
731 TRIO_PUBLIC int
732 trio_signbit
733 TRIO_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 */
750 TRIO_PUBLIC int
751 trio_fpclassify
752 TRIO_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 
775 static TRIO_CONST char *
776 getClassification
777 TRIO_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 
796 static void
797 print_class
798 TRIO_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 
main(TRIO_NOARGS)809 int 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