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