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
trio_nzero(TRIO_NOARGS)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
trio_pinf(TRIO_NOARGS)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
trio_ninf(TRIO_NOARGS)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
trio_nan(TRIO_NOARGS)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
main(TRIO_NOARGS)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
1258