blob: 3a65986b8577cdacdc082b77273022ba01ec6dc6 [file] [log] [blame]
Bjorn Reese45029602001-08-21 09:23:53 +00001/*************************************************************************
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 * 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 *
Bjorn Reese45029602001-08-21 09:23:53 +000032 ************************************************************************/
33
Daniel Veillardb7c29c32002-09-25 22:44:43 +000034/*
35 * TODO:
36 * o Put all the magic into trio_fpclassify_and_signbit(), and use this from
37 * trio_isnan() etc.
38 */
Bjorn Reese45029602001-08-21 09:23:53 +000039
Bjorn Reese45029602001-08-21 09:23:53 +000040/*************************************************************************
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
Bjorn Reese026d29f2002-01-19 15:40:18 +000053#if defined(TRIO_COMPILER_DECC)
54# include <fp_class.h>
55#endif
Bjorn Reese45029602001-08-21 09:23:53 +000056#include <assert.h>
57
Bjorn Reese026d29f2002-01-19 15:40:18 +000058#if defined(TRIO_DOCUMENTATION)
59# include "doc/doc_nan.h"
Bjorn Reese45029602001-08-21 09:23:53 +000060#endif
Bjorn Reese026d29f2002-01-19 15:40:18 +000061/** @addtogroup SpecialQuantities
62 @{
63*/
Bjorn Reese45029602001-08-21 09:23:53 +000064
65/*************************************************************************
66 * Definitions
67 */
68
Daniel Veillard21458c82002-03-27 16:12:22 +000069#define TRIO_TRUE (1 == 1)
70#define TRIO_FALSE (0 == 1)
71
Bjorn Reese45029602001-08-21 09:23:53 +000072/* We must enable IEEE floating-point on Alpha */
73#if defined(__alpha) && !defined(_IEEE_FP)
74# if defined(TRIO_COMPILER_DECC)
Bjorn Reese026d29f2002-01-19 15:40:18 +000075# if defined(TRIO_PLATFORM_VMS)
76# error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"
77# else
Daniel Veillardb7c29c32002-09-25 22:44:43 +000078# if !defined(_CFE)
79# error "Must be compiled with option -ieee"
80# endif
Bjorn Reese026d29f2002-01-19 15:40:18 +000081# endif
Bjorn Reese45029602001-08-21 09:23:53 +000082# elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))
83# error "Must be compiled with option -mieee"
84# endif
85#endif /* __alpha && ! _IEEE_FP */
86
87/*
88 * In ANSI/IEEE 754-1985 64-bits double format numbers have the
89 * following properties (amoungst others)
90 *
91 * o FLT_RADIX == 2: binary encoding
92 * o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used
93 * to indicate special numbers (e.g. NaN and Infinity), so the
94 * maximum exponent is 10 bits wide (2^10 == 1024).
95 * o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because
96 * numbers are normalized the initial binary 1 is represented
Daniel Veillardcbaf3992001-12-31 16:16:02 +000097 * implicitly (the so-called "hidden bit"), which leaves us with
Bjorn Reese45029602001-08-21 09:23:53 +000098 * the ability to represent 53 bits wide mantissa.
99 */
100#if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)
101# define USE_IEEE_754
102#endif
103
104
105/*************************************************************************
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000106 * Constants
Bjorn Reese45029602001-08-21 09:23:53 +0000107 */
108
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000109static TRIO_CONST char rcsid[] = "@(#)$Id$";
110
Bjorn Reese45029602001-08-21 09:23:53 +0000111#if defined(USE_IEEE_754)
112
113/*
114 * Endian-agnostic indexing macro.
115 *
116 * The value of internalEndianMagic, when converted into a 64-bit
Bjorn Reese026d29f2002-01-19 15:40:18 +0000117 * integer, becomes 0x0706050403020100 (we could have used a 64-bit
Bjorn Reese45029602001-08-21 09:23:53 +0000118 * integer value instead of a double, but not all platforms supports
119 * that type). The value is automatically encoded with the correct
120 * endianess by the compiler, which means that we can support any
121 * kind of endianess. The individual bytes are then used as an index
122 * for the IEEE 754 bit-patterns and masks.
123 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000124#define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])
Bjorn Reese45029602001-08-21 09:23:53 +0000125
Bjorn Reese026d29f2002-01-19 15:40:18 +0000126static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
Bjorn Reese45029602001-08-21 09:23:53 +0000127
128/* Mask for the exponent */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000129static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000130 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
131};
132
133/* Mask for the mantissa */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000134static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000135 0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
136};
137
Daniel Veillard21458c82002-03-27 16:12:22 +0000138/* Mask for the sign bit */
139static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
140 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
141};
142
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000143/* Bit-pattern for negative zero */
144static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
145 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
146};
147
Bjorn Reese45029602001-08-21 09:23:53 +0000148/* Bit-pattern for infinity */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000149static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000150 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
151};
152
153/* Bit-pattern for quiet NaN */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000154static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000155 0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
156};
157
158
159/*************************************************************************
Bjorn Reese026d29f2002-01-19 15:40:18 +0000160 * Functions
161 */
162
163/*
Bjorn Reese45029602001-08-21 09:23:53 +0000164 * trio_make_double
165 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000166TRIO_PRIVATE double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000167trio_make_double
168TRIO_ARGS1((values),
169 TRIO_CONST unsigned char *values)
Bjorn Reese45029602001-08-21 09:23:53 +0000170{
Bjorn Reese026d29f2002-01-19 15:40:18 +0000171 TRIO_VOLATILE double result;
Bjorn Reese45029602001-08-21 09:23:53 +0000172 int i;
173
174 for (i = 0; i < (int)sizeof(double); i++) {
Bjorn Reese026d29f2002-01-19 15:40:18 +0000175 ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
Bjorn Reese45029602001-08-21 09:23:53 +0000176 }
177 return result;
178}
179
Bjorn Reese026d29f2002-01-19 15:40:18 +0000180/*
Daniel Veillard21458c82002-03-27 16:12:22 +0000181 * trio_is_special_quantity
Bjorn Reese45029602001-08-21 09:23:53 +0000182 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000183TRIO_PRIVATE int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000184trio_is_special_quantity
185TRIO_ARGS2((number, has_mantissa),
186 double number,
187 int *has_mantissa)
Bjorn Reese45029602001-08-21 09:23:53 +0000188{
189 unsigned int i;
190 unsigned char current;
Daniel Veillard21458c82002-03-27 16:12:22 +0000191 int is_special_quantity = TRIO_TRUE;
Bjorn Reese45029602001-08-21 09:23:53 +0000192
193 *has_mantissa = 0;
194
195 for (i = 0; i < (unsigned int)sizeof(double); i++) {
196 current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
197 is_special_quantity
198 &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
199 *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
200 }
201 return is_special_quantity;
202}
203
Daniel Veillard21458c82002-03-27 16:12:22 +0000204/*
205 * trio_is_negative
206 */
207TRIO_PRIVATE int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000208trio_is_negative
209TRIO_ARGS1((number),
210 double number)
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000211{
212 unsigned int i;
Daniel Veillard21458c82002-03-27 16:12:22 +0000213 int is_negative = TRIO_FALSE;
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000214
215 for (i = 0; i < (unsigned int)sizeof(double); i++) {
Daniel Veillard21458c82002-03-27 16:12:22 +0000216 is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
217 & ieee_754_sign_mask[i]);
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000218 }
Daniel Veillard21458c82002-03-27 16:12:22 +0000219 return is_negative;
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000220}
221
Bjorn Reese45029602001-08-21 09:23:53 +0000222#endif /* USE_IEEE_754 */
223
224
Bjorn Reese026d29f2002-01-19 15:40:18 +0000225/**
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000226 Generate negative zero.
227
228 @return Floating-point representation of negative zero.
229*/
230TRIO_PUBLIC double
231trio_nzero(TRIO_NOARGS)
232{
233#if defined(USE_IEEE_754)
234 return trio_make_double(ieee_754_negzero_array);
235#else
236 TRIO_VOLATILE double zero = 0.0;
237
238 return -zero;
239#endif
240}
241
242/**
Bjorn Reese026d29f2002-01-19 15:40:18 +0000243 Generate positive infinity.
244
245 @return Floating-point representation of positive infinity.
246*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000247TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000248trio_pinf(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000249{
250 /* Cache the result */
251 static double result = 0.0;
252
253 if (result == 0.0) {
254
255#if defined(INFINITY) && defined(__STDC_IEC_559__)
256 result = (double)INFINITY;
257
258#elif defined(USE_IEEE_754)
259 result = trio_make_double(ieee_754_infinity_array);
260
261#else
262 /*
263 * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
264 * as infinity. Otherwise we have to resort to an overflow
265 * operation to generate infinity.
266 */
267# if defined(TRIO_PLATFORM_UNIX)
268 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
269# endif
270
271 result = HUGE_VAL;
272 if (HUGE_VAL == DBL_MAX) {
273 /* Force overflow */
274 result += HUGE_VAL;
275 }
276
277# if defined(TRIO_PLATFORM_UNIX)
278 signal(SIGFPE, signal_handler);
279# endif
280
281#endif
282 }
283 return result;
284}
285
Bjorn Reese026d29f2002-01-19 15:40:18 +0000286/**
287 Generate negative infinity.
288
289 @return Floating-point value of negative infinity.
290*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000291TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000292trio_ninf(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000293{
294 static double result = 0.0;
295
296 if (result == 0.0) {
297 /*
298 * Negative infinity is calculated by negating positive infinity,
299 * which can be done because it is legal to do calculations on
300 * infinity (for example, 1 / infinity == 0).
301 */
302 result = -trio_pinf();
303 }
304 return result;
305}
306
Bjorn Reese026d29f2002-01-19 15:40:18 +0000307/**
308 Generate NaN.
309
310 @return Floating-point representation of NaN.
311*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000312TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000313trio_nan(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000314{
315 /* Cache the result */
316 static double result = 0.0;
317
318 if (result == 0.0) {
319
320#if defined(TRIO_COMPILER_SUPPORTS_C99)
Bjorn Reese54d02fb2002-04-19 15:16:01 +0000321 result = nan("");
Bjorn Reese45029602001-08-21 09:23:53 +0000322
323#elif defined(NAN) && defined(__STDC_IEC_559__)
324 result = (double)NAN;
325
326#elif defined(USE_IEEE_754)
327 result = trio_make_double(ieee_754_qnan_array);
328
329#else
330 /*
331 * There are several ways to generate NaN. The one used here is
332 * to divide infinity by infinity. I would have preferred to add
333 * negative infinity to positive infinity, but that yields wrong
334 * result (infinity) on FreeBSD.
335 *
336 * This may fail if the hardware does not support NaN, or if
337 * the Invalid Operation floating-point exception is unmasked.
338 */
339# if defined(TRIO_PLATFORM_UNIX)
340 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
341# endif
342
343 result = trio_pinf() / trio_pinf();
344
345# if defined(TRIO_PLATFORM_UNIX)
346 signal(SIGFPE, signal_handler);
347# endif
348
349#endif
350 }
351 return result;
352}
353
Bjorn Reese026d29f2002-01-19 15:40:18 +0000354/**
355 Check for NaN.
356
357 @param number An arbitrary floating-point number.
358 @return Boolean value indicating whether or not the number is a NaN.
359*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000360TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000361trio_isnan
362TRIO_ARGS1((number),
363 double number)
Bjorn Reese45029602001-08-21 09:23:53 +0000364{
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000365#if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \
366 || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
Bjorn Reese45029602001-08-21 09:23:53 +0000367 /*
368 * C99 defines isnan() as a macro. UNIX95 defines isnan() as a
369 * function. This function was already present in XPG4, but this
370 * is a bit tricky to detect with compiler defines, so we choose
371 * the conservative approach and only use it for UNIX95.
372 */
373 return isnan(number);
374
375#elif defined(TRIO_COMPILER_MSVC)
376 /*
Daniel Veillard21458c82002-03-27 16:12:22 +0000377 * MSVC has an _isnan() function
Bjorn Reese45029602001-08-21 09:23:53 +0000378 */
379 return _isnan(number);
380
381#elif defined(USE_IEEE_754)
382 /*
383 * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
384 * pattern, and a non-empty mantissa.
385 */
386 int has_mantissa;
387 int is_special_quantity;
388
389 is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
390
391 return (is_special_quantity && has_mantissa);
392
393#else
394 /*
395 * Fallback solution
396 */
397 int status;
398 double integral, fraction;
399
400# if defined(TRIO_PLATFORM_UNIX)
401 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
402# endif
403
404 status = (/*
405 * NaN is the only number which does not compare to itself
406 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000407 ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
Bjorn Reese45029602001-08-21 09:23:53 +0000408 /*
409 * Fallback solution if NaN compares to NaN
410 */
411 ((number != 0.0) &&
412 (fraction = modf(number, &integral),
413 integral == fraction)));
414
415# if defined(TRIO_PLATFORM_UNIX)
416 signal(SIGFPE, signal_handler);
417# endif
418
419 return status;
420
421#endif
422}
423
Bjorn Reese026d29f2002-01-19 15:40:18 +0000424/**
425 Check for infinity.
426
427 @param number An arbitrary floating-point number.
428 @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.
429*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000430TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000431trio_isinf
432TRIO_ARGS1((number),
433 double number)
Bjorn Reese45029602001-08-21 09:23:53 +0000434{
435#if defined(TRIO_COMPILER_DECC)
436 /*
437 * DECC has an isinf() macro, but it works differently than that
438 * of C99, so we use the fp_class() function instead.
439 */
440 return ((fp_class(number) == FP_POS_INF)
441 ? 1
442 : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000443
Bjorn Reese45029602001-08-21 09:23:53 +0000444#elif defined(isinf)
445 /*
446 * C99 defines isinf() as a macro.
447 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000448 return isinf(number)
449 ? ((number > 0.0) ? 1 : -1)
450 : 0;
Bjorn Reese45029602001-08-21 09:23:53 +0000451
452#elif defined(TRIO_COMPILER_MSVC)
453 /*
454 * MSVC has an _fpclass() function that can be used to detect infinity.
455 */
456 return ((_fpclass(number) == _FPCLASS_PINF)
457 ? 1
458 : ((_fpclass(number) == _FPCLASS_NINF) ? -1 : 0));
459
460#elif defined(USE_IEEE_754)
461 /*
462 * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
463 * pattern, and an empty mantissa.
464 */
465 int has_mantissa;
466 int is_special_quantity;
467
468 is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
469
470 return (is_special_quantity && !has_mantissa)
471 ? ((number < 0.0) ? -1 : 1)
472 : 0;
473
474#else
475 /*
476 * Fallback solution.
477 */
478 int status;
479
480# if defined(TRIO_PLATFORM_UNIX)
481 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
482# endif
483
484 double infinity = trio_pinf();
485
486 status = ((number == infinity)
487 ? 1
488 : ((number == -infinity) ? -1 : 0));
489
490# if defined(TRIO_PLATFORM_UNIX)
491 signal(SIGFPE, signal_handler);
492# endif
493
494 return status;
495
496#endif
497}
498
Daniel Veillard21458c82002-03-27 16:12:22 +0000499
500/**
501 Check for finity.
502
503 @param number An arbitrary floating-point number.
504 @return Boolean value indicating whether or not the number is a finite.
505*/
506TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000507trio_isfinite
508TRIO_ARGS1((number),
509 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000510{
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000511#if defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isfinite)
Daniel Veillard21458c82002-03-27 16:12:22 +0000512 /*
513 * C99 defines isfinite() as a macro.
514 */
515 return isfinite(number);
516
517#elif defined(TRIO_COMPILER_MSVC)
518 /*
519 * MSVC uses _finite().
520 */
521 return _finite(number);
522
523#elif defined(USE_IEEE_754)
524 /*
525 * Examine IEEE 754 bit-pattern. For finity we do not care about the
526 * mantissa.
527 */
528 int dummy;
529
530 return (! trio_is_special_quantity(number, &dummy));
531
532#else
533 /*
534 * Fallback solution.
535 */
536 return ((trio_isinf(number) == 0) && (trio_isnan(number) == 0));
537
538#endif
539}
540
541/*
542 * The sign of NaN is always false
543 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000544TRIO_PUBLIC int
545trio_fpclassify_and_signbit
546TRIO_ARGS2((number, is_negative),
547 double number,
548 int *is_negative)
Daniel Veillard21458c82002-03-27 16:12:22 +0000549{
550#if defined(fpclassify) && defined(signbit)
551 /*
552 * C99 defines fpclassify() and signbit() as a macros
553 */
554 *is_negative = signbit(number);
555 switch (fpclassify(number)) {
556 case FP_NAN:
557 return TRIO_FP_NAN;
558 case FP_INFINITE:
559 return TRIO_FP_INFINITE;
560 case FP_SUBNORMAL:
561 return TRIO_FP_SUBNORMAL;
562 case FP_ZERO:
563 return TRIO_FP_ZERO;
564 default:
565 return TRIO_FP_NORMAL;
566 }
567
568#elif defined(TRIO_COMPILER_DECC)
569 /*
570 * DECC has an fp_class() function.
571 */
572 switch (fp_class(number)) {
573 case FP_QNAN:
574 case FP_SNAN:
575 *is_negative = TRIO_FALSE; /* NaN has no sign */
576 return TRIO_FP_NAN;
577 case FP_POS_INF:
578 *is_negative = TRIO_FALSE;
579 return TRIO_FP_INFINITE;
580 case FP_NEG_INF:
581 *is_negative = TRIO_TRUE;
582 return TRIO_FP_INFINITE;
583 case FP_POS_DENORM:
584 *is_negative = TRIO_FALSE;
585 return TRIO_FP_SUBNORMAL;
586 case FP_NEG_DENORM:
587 *is_negative = TRIO_TRUE;
588 return TRIO_FP_SUBNORMAL;
589 case FP_POS_ZERO:
590 *is_negative = TRIO_FALSE;
591 return TRIO_FP_ZERO;
592 case FP_NEG_ZERO:
593 *is_negative = TRIO_TRUE;
594 return TRIO_FP_ZERO;
595 case FP_POS_NORM:
596 *is_negative = TRIO_FALSE;
597 return TRIO_FP_NORMAL;
598 case FP_NEG_NORM:
599 *is_negative = TRIO_TRUE;
600 return TRIO_FP_NORMAL;
601 default:
602 /* Just in case... */
603 *is_negative = (number < 0.0);
604 return TRIO_FP_NORMAL;
605 }
606
607#elif defined(TRIO_COMPILER_MSVC)
608 /*
609 * MSVC has an _fpclass() function.
610 */
611 switch (_fpclass(number)) {
612 case _FPCLASS_QNAN:
613 case _FPCLASS_SNAN:
614 *is_negative = TRIO_FALSE;
615 return TRIO_FP_NAN;
616 case _FPCLASS_PINF:
617 *is_negative = TRIO_FALSE;
618 return TRIO_FP_INFINITE;
619 case _FPCLASS_NINF:
620 *is_negative = TRIO_TRUE;
621 return TRIO_FP_INFINITE;
622 case _FPCLASS_PD:
623 *is_negative = TRIO_FALSE;
624 return TRIO_FP_SUBNORMAL;
625 case _FPCLASS_ND:
626 *is_negative = TRIO_TRUE;
627 return TRIO_FP_SUBNORMAL;
628 case _FPCLASS_PZ:
629 *is_negative = TRIO_FALSE;
630 return TRIO_FP_ZERO;
631 case _FPCLASS_NZ:
632 *is_negative = TRIO_TRUE;
633 return TRIO_FP_ZERO;
634 case _FPCLASS_PN:
635 *is_negative = TRIO_FALSE;
636 return TRIO_FP_NORMAL;
637 case _FPCLASS_NN:
638 *is_negative = TRIO_TRUE;
639 return TRIO_FP_NORMAL;
640 default:
641 /* Just in case... */
642 *is_negative = (number < 0.0);
643 return TRIO_FP_NORMAL;
644 }
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000645
646#elif defined(FP_PLUS_NORM) || defined(__hpux)
647
648 /*
649 * HP-UX 9.x and 10.x have an fpclassify() function, that is different
650 * from the C99 fpclassify() macro supported on HP-UX 11.x.
651 */
652 switch (fpclassify(number)) {
653 case FP_QNAN:
654 case FP_SNAN:
655 *is_negative = TRIO_FALSE;
656 return TRIO_FP_NAN;
657 case FP_PLUS_INF:
658 *is_negative = TRIO_FALSE;
659 return TRIO_FP_INFINITE;
660 case FP_MINUS_INF:
661 *is_negative = TRIO_TRUE;
662 return TRIO_FP_INFINITE;
663 case FP_PLUS_DENORM:
664 *is_negative = TRIO_FALSE;
665 return TRIO_FP_SUBNORMAL;
666 case FP_MINUS_DENORM:
667 *is_negative = TRIO_TRUE;
668 return TRIO_FP_SUBNORMAL;
669 case FP_PLUS_ZERO:
670 *is_negative = TRIO_FALSE;
671 return TRIO_FP_ZERO;
672 case FP_MINUS_ZERO:
673 *is_negative = TRIO_TRUE;
674 return TRIO_FP_ZERO;
675 case FP_PLUS_NORM:
676 *is_negative = TRIO_FALSE;
677 return TRIO_FP_NORMAL;
678 case FP_MINUS_NORM:
679 *is_negative = TRIO_TRUE;
680 return TRIO_FP_NORMAL;
681 default:
682 assert(0);
683 }
684
Daniel Veillard21458c82002-03-27 16:12:22 +0000685#else
686 /*
687 * Fallback solution.
688 */
689 int rc;
690
691 if (number == 0.0) {
692 /*
693 * In IEEE 754 the sign of zero is ignored in comparisons, so we
694 * have to handle this as a special case by examining the sign bit
695 * directly.
696 */
697#if defined(USE_IEEE_754)
698 *is_negative = trio_is_negative(number);
699#else
700 *is_negative = TRIO_FALSE; /* FIXME */
701#endif
702 return TRIO_FP_ZERO;
703 }
704 if (trio_isnan(number)) {
705 *is_negative = TRIO_FALSE;
706 return TRIO_FP_NAN;
707 }
708 if ((rc = trio_isinf(number))) {
709 *is_negative = (rc == -1);
710 return TRIO_FP_INFINITE;
711 }
712 if ((number > 0.0) && (number < DBL_MIN)) {
713 *is_negative = TRIO_FALSE;
714 return TRIO_FP_SUBNORMAL;
715 }
716 if ((number < 0.0) && (number > -DBL_MIN)) {
717 *is_negative = TRIO_TRUE;
718 return TRIO_FP_SUBNORMAL;
719 }
720 *is_negative = (number < 0.0);
721 return TRIO_FP_NORMAL;
722
723#endif
724}
725
726/**
727 Examine the sign of a number.
728
729 @param number An arbitrary floating-point number.
730 @return Boolean value indicating whether or not the number has the
731 sign bit set (i.e. is negative).
732*/
733TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000734trio_signbit
735TRIO_ARGS1((number),
736 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000737{
738 int is_negative;
739
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000740 (void)trio_fpclassify_and_signbit(number, &is_negative);
Daniel Veillard21458c82002-03-27 16:12:22 +0000741 return is_negative;
742}
743
744/**
745 Examine the class of a number.
746
747 @param number An arbitrary floating-point number.
748 @return Enumerable value indicating the class of @p number
749*/
750TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000751trio_fpclassify
752TRIO_ARGS1((number),
753 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000754{
755 int dummy;
756
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000757 return trio_fpclassify_and_signbit(number, &dummy);
Daniel Veillard21458c82002-03-27 16:12:22 +0000758}
759
760
Bjorn Reese026d29f2002-01-19 15:40:18 +0000761/** @} SpecialQuantities */
762
Bjorn Reese45029602001-08-21 09:23:53 +0000763/*************************************************************************
Bjorn Reese026d29f2002-01-19 15:40:18 +0000764 * For test purposes.
765 *
766 * Add the following compiler option to include this test code.
767 *
768 * Unix : -DSTANDALONE
769 * VMS : /DEFINE=(STANDALONE)
Bjorn Reese45029602001-08-21 09:23:53 +0000770 */
771#if defined(STANDALONE)
772# include <stdio.h>
773
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000774static TRIO_CONST char *
775getClassification
Daniel Veillarde645e8c2002-10-22 17:35:37 +0000776TRIO_ARGS1((type),
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000777 int type)
Daniel Veillard21458c82002-03-27 16:12:22 +0000778{
779 switch (type) {
780 case TRIO_FP_INFINITE:
781 return "FP_INFINITE";
782 case TRIO_FP_NAN:
783 return "FP_NAN";
784 case TRIO_FP_NORMAL:
785 return "FP_NORMAL";
786 case TRIO_FP_SUBNORMAL:
787 return "FP_SUBNORMAL";
788 case TRIO_FP_ZERO:
789 return "FP_ZERO";
790 default:
791 return "FP_UNKNOWN";
792 }
793}
794
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000795static void
796print_class
Daniel Veillarde645e8c2002-10-22 17:35:37 +0000797TRIO_ARGS2((prefix, number),
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000798 TRIO_CONST char *prefix,
799 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000800{
801 printf("%-6s: %s %-15s %g\n",
802 prefix,
803 trio_signbit(number) ? "-" : "+",
804 getClassification(trio_fpclassify(number)),
805 number);
806}
807
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000808int main(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000809{
810 double my_nan;
811 double my_pinf;
812 double my_ninf;
813# if defined(TRIO_PLATFORM_UNIX)
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000814 void (*signal_handler) TRIO_PROTO((int));
Bjorn Reese45029602001-08-21 09:23:53 +0000815# endif
816
817 my_nan = trio_nan();
818 my_pinf = trio_pinf();
819 my_ninf = trio_ninf();
820
Daniel Veillard21458c82002-03-27 16:12:22 +0000821 print_class("Nan", my_nan);
822 print_class("PInf", my_pinf);
823 print_class("NInf", my_ninf);
824 print_class("PZero", 0.0);
825 print_class("NZero", -0.0);
826 print_class("PNorm", 1.0);
827 print_class("NNorm", -1.0);
828 print_class("PSub", 1.01e-307 - 1.00e-307);
829 print_class("NSub", 1.00e-307 - 1.01e-307);
830
Bjorn Reese45029602001-08-21 09:23:53 +0000831 printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
832 my_nan,
833 ((unsigned char *)&my_nan)[0],
834 ((unsigned char *)&my_nan)[1],
835 ((unsigned char *)&my_nan)[2],
836 ((unsigned char *)&my_nan)[3],
837 ((unsigned char *)&my_nan)[4],
838 ((unsigned char *)&my_nan)[5],
839 ((unsigned char *)&my_nan)[6],
840 ((unsigned char *)&my_nan)[7],
841 trio_isnan(my_nan), trio_isinf(my_nan));
842 printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
843 my_pinf,
844 ((unsigned char *)&my_pinf)[0],
845 ((unsigned char *)&my_pinf)[1],
846 ((unsigned char *)&my_pinf)[2],
847 ((unsigned char *)&my_pinf)[3],
848 ((unsigned char *)&my_pinf)[4],
849 ((unsigned char *)&my_pinf)[5],
850 ((unsigned char *)&my_pinf)[6],
851 ((unsigned char *)&my_pinf)[7],
852 trio_isnan(my_pinf), trio_isinf(my_pinf));
853 printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
854 my_ninf,
855 ((unsigned char *)&my_ninf)[0],
856 ((unsigned char *)&my_ninf)[1],
857 ((unsigned char *)&my_ninf)[2],
858 ((unsigned char *)&my_ninf)[3],
859 ((unsigned char *)&my_ninf)[4],
860 ((unsigned char *)&my_ninf)[5],
861 ((unsigned char *)&my_ninf)[6],
862 ((unsigned char *)&my_ninf)[7],
863 trio_isnan(my_ninf), trio_isinf(my_ninf));
864
865# if defined(TRIO_PLATFORM_UNIX)
866 signal_handler = signal(SIGFPE, SIG_IGN);
867# endif
868
869 my_pinf = DBL_MAX + DBL_MAX;
870 my_ninf = -my_pinf;
871 my_nan = my_pinf / my_pinf;
872
873# if defined(TRIO_PLATFORM_UNIX)
874 signal(SIGFPE, signal_handler);
875# endif
876
877 printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
878 my_nan,
879 ((unsigned char *)&my_nan)[0],
880 ((unsigned char *)&my_nan)[1],
881 ((unsigned char *)&my_nan)[2],
882 ((unsigned char *)&my_nan)[3],
883 ((unsigned char *)&my_nan)[4],
884 ((unsigned char *)&my_nan)[5],
885 ((unsigned char *)&my_nan)[6],
886 ((unsigned char *)&my_nan)[7],
887 trio_isnan(my_nan), trio_isinf(my_nan));
888 printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
889 my_pinf,
890 ((unsigned char *)&my_pinf)[0],
891 ((unsigned char *)&my_pinf)[1],
892 ((unsigned char *)&my_pinf)[2],
893 ((unsigned char *)&my_pinf)[3],
894 ((unsigned char *)&my_pinf)[4],
895 ((unsigned char *)&my_pinf)[5],
896 ((unsigned char *)&my_pinf)[6],
897 ((unsigned char *)&my_pinf)[7],
898 trio_isnan(my_pinf), trio_isinf(my_pinf));
899 printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
900 my_ninf,
901 ((unsigned char *)&my_ninf)[0],
902 ((unsigned char *)&my_ninf)[1],
903 ((unsigned char *)&my_ninf)[2],
904 ((unsigned char *)&my_ninf)[3],
905 ((unsigned char *)&my_ninf)[4],
906 ((unsigned char *)&my_ninf)[5],
907 ((unsigned char *)&my_ninf)[6],
908 ((unsigned char *)&my_ninf)[7],
909 trio_isnan(my_ninf), trio_isinf(my_ninf));
910
911 return 0;
912}
913#endif