blob: 95baae1e022ca03aa06d5ace5286e5449abbb94e [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)
William M. Brack476cd962003-08-13 11:09:42 +000054# if defined(__linux__)
55# include <cpml.h>
56# else
57# include <fp_class.h>
58# endif
Bjorn Reese026d29f2002-01-19 15:40:18 +000059#endif
Bjorn Reese45029602001-08-21 09:23:53 +000060#include <assert.h>
61
Bjorn Reese026d29f2002-01-19 15:40:18 +000062#if defined(TRIO_DOCUMENTATION)
63# include "doc/doc_nan.h"
Bjorn Reese45029602001-08-21 09:23:53 +000064#endif
Bjorn Reese026d29f2002-01-19 15:40:18 +000065/** @addtogroup SpecialQuantities
66 @{
67*/
Bjorn Reese45029602001-08-21 09:23:53 +000068
69/*************************************************************************
70 * Definitions
71 */
72
Daniel Veillard21458c82002-03-27 16:12:22 +000073#define TRIO_TRUE (1 == 1)
74#define TRIO_FALSE (0 == 1)
75
Daniel Veillarda48ed3d2003-04-03 15:28:28 +000076/*
77 * We must enable IEEE floating-point on Alpha
78 */
Bjorn Reese45029602001-08-21 09:23:53 +000079#if defined(__alpha) && !defined(_IEEE_FP)
80# if defined(TRIO_COMPILER_DECC)
Bjorn Reese026d29f2002-01-19 15:40:18 +000081# if defined(TRIO_PLATFORM_VMS)
82# error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"
83# else
Daniel Veillardb7c29c32002-09-25 22:44:43 +000084# if !defined(_CFE)
85# error "Must be compiled with option -ieee"
86# endif
Bjorn Reese026d29f2002-01-19 15:40:18 +000087# endif
Bjorn Reese45029602001-08-21 09:23:53 +000088# elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))
89# error "Must be compiled with option -mieee"
90# endif
91#endif /* __alpha && ! _IEEE_FP */
92
93/*
94 * In ANSI/IEEE 754-1985 64-bits double format numbers have the
95 * following properties (amoungst others)
96 *
97 * o FLT_RADIX == 2: binary encoding
98 * o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used
99 * to indicate special numbers (e.g. NaN and Infinity), so the
100 * maximum exponent is 10 bits wide (2^10 == 1024).
101 * o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because
102 * numbers are normalized the initial binary 1 is represented
Daniel Veillardcbaf3992001-12-31 16:16:02 +0000103 * implicitly (the so-called "hidden bit"), which leaves us with
Bjorn Reese45029602001-08-21 09:23:53 +0000104 * the ability to represent 53 bits wide mantissa.
105 */
106#if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)
107# define USE_IEEE_754
108#endif
109
110
111/*************************************************************************
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000112 * Constants
Bjorn Reese45029602001-08-21 09:23:53 +0000113 */
114
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000115static TRIO_CONST char rcsid[] = "@(#)$Id$";
116
Bjorn Reese45029602001-08-21 09:23:53 +0000117#if defined(USE_IEEE_754)
118
119/*
120 * Endian-agnostic indexing macro.
121 *
122 * The value of internalEndianMagic, when converted into a 64-bit
Bjorn Reese026d29f2002-01-19 15:40:18 +0000123 * integer, becomes 0x0706050403020100 (we could have used a 64-bit
Bjorn Reese45029602001-08-21 09:23:53 +0000124 * integer value instead of a double, but not all platforms supports
125 * that type). The value is automatically encoded with the correct
126 * endianess by the compiler, which means that we can support any
127 * kind of endianess. The individual bytes are then used as an index
128 * for the IEEE 754 bit-patterns and masks.
129 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000130#define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])
Bjorn Reese45029602001-08-21 09:23:53 +0000131
Daniel Veillardda3fee42008-09-01 13:08:57 +0000132#if (defined(__BORLANDC__) && __BORLANDC__ >= 0x0590)
133static TRIO_CONST double internalEndianMagic = 7.949928895127362e-275;
134#else
Bjorn Reese026d29f2002-01-19 15:40:18 +0000135static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
Daniel Veillardda3fee42008-09-01 13:08:57 +0000136#endif
Bjorn Reese45029602001-08-21 09:23:53 +0000137
138/* Mask for the exponent */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000139static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000140 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
141};
142
143/* Mask for the mantissa */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000144static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000145 0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
146};
147
Daniel Veillard21458c82002-03-27 16:12:22 +0000148/* Mask for the sign bit */
149static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
150 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
151};
152
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000153/* Bit-pattern for negative zero */
154static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
155 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
156};
157
Bjorn Reese45029602001-08-21 09:23:53 +0000158/* Bit-pattern for infinity */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000159static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000160 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
161};
162
163/* Bit-pattern for quiet NaN */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000164static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
Bjorn Reese45029602001-08-21 09:23:53 +0000165 0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
166};
167
168
169/*************************************************************************
Bjorn Reese026d29f2002-01-19 15:40:18 +0000170 * Functions
171 */
172
173/*
Bjorn Reese45029602001-08-21 09:23:53 +0000174 * trio_make_double
175 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000176TRIO_PRIVATE double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000177trio_make_double
178TRIO_ARGS1((values),
179 TRIO_CONST unsigned char *values)
Bjorn Reese45029602001-08-21 09:23:53 +0000180{
Bjorn Reese026d29f2002-01-19 15:40:18 +0000181 TRIO_VOLATILE double result;
Bjorn Reese45029602001-08-21 09:23:53 +0000182 int i;
183
184 for (i = 0; i < (int)sizeof(double); i++) {
Bjorn Reese026d29f2002-01-19 15:40:18 +0000185 ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
Bjorn Reese45029602001-08-21 09:23:53 +0000186 }
187 return result;
188}
189
Bjorn Reese026d29f2002-01-19 15:40:18 +0000190/*
Daniel Veillard21458c82002-03-27 16:12:22 +0000191 * trio_is_special_quantity
Bjorn Reese45029602001-08-21 09:23:53 +0000192 */
Bjorn Reese026d29f2002-01-19 15:40:18 +0000193TRIO_PRIVATE int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000194trio_is_special_quantity
195TRIO_ARGS2((number, has_mantissa),
196 double number,
197 int *has_mantissa)
Bjorn Reese45029602001-08-21 09:23:53 +0000198{
199 unsigned int i;
200 unsigned char current;
Daniel Veillard21458c82002-03-27 16:12:22 +0000201 int is_special_quantity = TRIO_TRUE;
Bjorn Reese45029602001-08-21 09:23:53 +0000202
203 *has_mantissa = 0;
204
205 for (i = 0; i < (unsigned int)sizeof(double); i++) {
206 current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
207 is_special_quantity
208 &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
209 *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
210 }
211 return is_special_quantity;
212}
213
Daniel Veillard21458c82002-03-27 16:12:22 +0000214/*
215 * trio_is_negative
216 */
217TRIO_PRIVATE int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000218trio_is_negative
219TRIO_ARGS1((number),
220 double number)
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000221{
222 unsigned int i;
Daniel Veillard21458c82002-03-27 16:12:22 +0000223 int is_negative = TRIO_FALSE;
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000224
225 for (i = 0; i < (unsigned int)sizeof(double); i++) {
Daniel Veillard21458c82002-03-27 16:12:22 +0000226 is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
227 & ieee_754_sign_mask[i]);
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000228 }
Daniel Veillard21458c82002-03-27 16:12:22 +0000229 return is_negative;
Daniel Veillard5fc1f082002-03-27 09:05:40 +0000230}
231
Bjorn Reese45029602001-08-21 09:23:53 +0000232#endif /* USE_IEEE_754 */
233
234
Bjorn Reese026d29f2002-01-19 15:40:18 +0000235/**
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000236 Generate negative zero.
237
238 @return Floating-point representation of negative zero.
239*/
240TRIO_PUBLIC double
241trio_nzero(TRIO_NOARGS)
242{
243#if defined(USE_IEEE_754)
244 return trio_make_double(ieee_754_negzero_array);
245#else
246 TRIO_VOLATILE double zero = 0.0;
247
248 return -zero;
249#endif
250}
251
252/**
Bjorn Reese026d29f2002-01-19 15:40:18 +0000253 Generate positive infinity.
254
255 @return Floating-point representation of positive infinity.
256*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000257TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000258trio_pinf(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000259{
260 /* Cache the result */
261 static double result = 0.0;
262
263 if (result == 0.0) {
264
265#if defined(INFINITY) && defined(__STDC_IEC_559__)
266 result = (double)INFINITY;
267
268#elif defined(USE_IEEE_754)
269 result = trio_make_double(ieee_754_infinity_array);
270
271#else
272 /*
273 * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
274 * as infinity. Otherwise we have to resort to an overflow
275 * operation to generate infinity.
276 */
277# if defined(TRIO_PLATFORM_UNIX)
278 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
279# endif
280
281 result = HUGE_VAL;
282 if (HUGE_VAL == DBL_MAX) {
283 /* Force overflow */
284 result += HUGE_VAL;
285 }
286
287# if defined(TRIO_PLATFORM_UNIX)
288 signal(SIGFPE, signal_handler);
289# endif
290
291#endif
292 }
293 return result;
294}
295
Bjorn Reese026d29f2002-01-19 15:40:18 +0000296/**
297 Generate negative infinity.
298
299 @return Floating-point value of negative infinity.
300*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000301TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000302trio_ninf(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000303{
304 static double result = 0.0;
305
306 if (result == 0.0) {
307 /*
308 * Negative infinity is calculated by negating positive infinity,
309 * which can be done because it is legal to do calculations on
310 * infinity (for example, 1 / infinity == 0).
311 */
312 result = -trio_pinf();
313 }
314 return result;
315}
316
Bjorn Reese026d29f2002-01-19 15:40:18 +0000317/**
318 Generate NaN.
319
320 @return Floating-point representation of NaN.
321*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000322TRIO_PUBLIC double
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000323trio_nan(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000324{
325 /* Cache the result */
326 static double result = 0.0;
327
328 if (result == 0.0) {
329
330#if defined(TRIO_COMPILER_SUPPORTS_C99)
Bjorn Reese54d02fb2002-04-19 15:16:01 +0000331 result = nan("");
Bjorn Reese45029602001-08-21 09:23:53 +0000332
333#elif defined(NAN) && defined(__STDC_IEC_559__)
334 result = (double)NAN;
335
336#elif defined(USE_IEEE_754)
337 result = trio_make_double(ieee_754_qnan_array);
338
339#else
340 /*
341 * There are several ways to generate NaN. The one used here is
342 * to divide infinity by infinity. I would have preferred to add
343 * negative infinity to positive infinity, but that yields wrong
344 * result (infinity) on FreeBSD.
345 *
346 * This may fail if the hardware does not support NaN, or if
347 * the Invalid Operation floating-point exception is unmasked.
348 */
349# if defined(TRIO_PLATFORM_UNIX)
350 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
351# endif
352
353 result = trio_pinf() / trio_pinf();
354
355# if defined(TRIO_PLATFORM_UNIX)
356 signal(SIGFPE, signal_handler);
357# endif
358
359#endif
360 }
361 return result;
362}
363
Bjorn Reese026d29f2002-01-19 15:40:18 +0000364/**
365 Check for NaN.
366
367 @param number An arbitrary floating-point number.
368 @return Boolean value indicating whether or not the number is a NaN.
369*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000370TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000371trio_isnan
372TRIO_ARGS1((number),
373 double number)
Bjorn Reese45029602001-08-21 09:23:53 +0000374{
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000375#if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \
376 || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
Bjorn Reese45029602001-08-21 09:23:53 +0000377 /*
378 * C99 defines isnan() as a macro. UNIX95 defines isnan() as a
379 * function. This function was already present in XPG4, but this
380 * is a bit tricky to detect with compiler defines, so we choose
381 * the conservative approach and only use it for UNIX95.
382 */
383 return isnan(number);
384
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000385#elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
Bjorn Reese45029602001-08-21 09:23:53 +0000386 /*
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000387 * Microsoft Visual C++ and Borland C++ Builder have an _isnan()
388 * function.
Bjorn Reese45029602001-08-21 09:23:53 +0000389 */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000390 return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;
Bjorn Reese45029602001-08-21 09:23:53 +0000391
392#elif defined(USE_IEEE_754)
393 /*
394 * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
395 * pattern, and a non-empty mantissa.
396 */
397 int has_mantissa;
398 int is_special_quantity;
399
400 is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
401
402 return (is_special_quantity && has_mantissa);
403
404#else
405 /*
406 * Fallback solution
407 */
408 int status;
409 double integral, fraction;
410
411# if defined(TRIO_PLATFORM_UNIX)
412 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
413# endif
414
415 status = (/*
416 * NaN is the only number which does not compare to itself
417 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000418 ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
Bjorn Reese45029602001-08-21 09:23:53 +0000419 /*
420 * Fallback solution if NaN compares to NaN
421 */
422 ((number != 0.0) &&
423 (fraction = modf(number, &integral),
424 integral == fraction)));
425
426# if defined(TRIO_PLATFORM_UNIX)
427 signal(SIGFPE, signal_handler);
428# endif
429
430 return status;
431
432#endif
433}
434
Bjorn Reese026d29f2002-01-19 15:40:18 +0000435/**
436 Check for infinity.
437
438 @param number An arbitrary floating-point number.
439 @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.
440*/
Daniel Veillardcda96922001-08-21 10:56:31 +0000441TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000442trio_isinf
443TRIO_ARGS1((number),
444 double number)
Bjorn Reese45029602001-08-21 09:23:53 +0000445{
William M. Brack476cd962003-08-13 11:09:42 +0000446#if defined(TRIO_COMPILER_DECC) && !defined(__linux__)
Bjorn Reese45029602001-08-21 09:23:53 +0000447 /*
448 * DECC has an isinf() macro, but it works differently than that
449 * of C99, so we use the fp_class() function instead.
450 */
451 return ((fp_class(number) == FP_POS_INF)
452 ? 1
453 : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000454
Bjorn Reese45029602001-08-21 09:23:53 +0000455#elif defined(isinf)
456 /*
457 * C99 defines isinf() as a macro.
458 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000459 return isinf(number)
460 ? ((number > 0.0) ? 1 : -1)
461 : 0;
Bjorn Reese45029602001-08-21 09:23:53 +0000462
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000463#elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
Bjorn Reese45029602001-08-21 09:23:53 +0000464 /*
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000465 * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
466 * function that can be used to detect infinity.
Bjorn Reese45029602001-08-21 09:23:53 +0000467 */
468 return ((_fpclass(number) == _FPCLASS_PINF)
469 ? 1
470 : ((_fpclass(number) == _FPCLASS_NINF) ? -1 : 0));
471
472#elif defined(USE_IEEE_754)
473 /*
474 * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
475 * pattern, and an empty mantissa.
476 */
477 int has_mantissa;
478 int is_special_quantity;
479
480 is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
481
482 return (is_special_quantity && !has_mantissa)
483 ? ((number < 0.0) ? -1 : 1)
484 : 0;
485
486#else
487 /*
488 * Fallback solution.
489 */
490 int status;
491
492# if defined(TRIO_PLATFORM_UNIX)
493 void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
494# endif
495
496 double infinity = trio_pinf();
497
498 status = ((number == infinity)
499 ? 1
500 : ((number == -infinity) ? -1 : 0));
501
502# if defined(TRIO_PLATFORM_UNIX)
503 signal(SIGFPE, signal_handler);
504# endif
505
506 return status;
507
508#endif
509}
510
William M. Bracka71a8ef2003-08-06 04:43:55 +0000511#if 0
512 /* Temporary fix - this routine is not used anywhere */
Daniel Veillard21458c82002-03-27 16:12:22 +0000513/**
514 Check for finity.
515
516 @param number An arbitrary floating-point number.
517 @return Boolean value indicating whether or not the number is a finite.
518*/
519TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000520trio_isfinite
521TRIO_ARGS1((number),
522 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000523{
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000524#if defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isfinite)
Daniel Veillard21458c82002-03-27 16:12:22 +0000525 /*
526 * C99 defines isfinite() as a macro.
527 */
528 return isfinite(number);
529
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000530#elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
Daniel Veillard21458c82002-03-27 16:12:22 +0000531 /*
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000532 * Microsoft Visual C++ and Borland C++ Builder use _finite().
Daniel Veillard21458c82002-03-27 16:12:22 +0000533 */
534 return _finite(number);
535
536#elif defined(USE_IEEE_754)
537 /*
538 * Examine IEEE 754 bit-pattern. For finity we do not care about the
539 * mantissa.
540 */
541 int dummy;
542
543 return (! trio_is_special_quantity(number, &dummy));
544
545#else
546 /*
547 * Fallback solution.
548 */
549 return ((trio_isinf(number) == 0) && (trio_isnan(number) == 0));
550
551#endif
552}
553
William M. Bracka71a8ef2003-08-06 04:43:55 +0000554#endif
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000555
Daniel Veillard21458c82002-03-27 16:12:22 +0000556/*
557 * The sign of NaN is always false
558 */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000559TRIO_PUBLIC int
560trio_fpclassify_and_signbit
561TRIO_ARGS2((number, is_negative),
562 double number,
563 int *is_negative)
Daniel Veillard21458c82002-03-27 16:12:22 +0000564{
565#if defined(fpclassify) && defined(signbit)
566 /*
567 * C99 defines fpclassify() and signbit() as a macros
568 */
569 *is_negative = signbit(number);
570 switch (fpclassify(number)) {
571 case FP_NAN:
572 return TRIO_FP_NAN;
573 case FP_INFINITE:
574 return TRIO_FP_INFINITE;
575 case FP_SUBNORMAL:
576 return TRIO_FP_SUBNORMAL;
577 case FP_ZERO:
578 return TRIO_FP_ZERO;
579 default:
580 return TRIO_FP_NORMAL;
581 }
582
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000583#else
584# if defined(TRIO_COMPILER_DECC)
Daniel Veillard21458c82002-03-27 16:12:22 +0000585 /*
586 * DECC has an fp_class() function.
587 */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000588# define TRIO_FPCLASSIFY(n) fp_class(n)
589# define TRIO_QUIET_NAN FP_QNAN
590# define TRIO_SIGNALLING_NAN FP_SNAN
591# define TRIO_POSITIVE_INFINITY FP_POS_INF
592# define TRIO_NEGATIVE_INFINITY FP_NEG_INF
593# define TRIO_POSITIVE_SUBNORMAL FP_POS_DENORM
594# define TRIO_NEGATIVE_SUBNORMAL FP_NEG_DENORM
595# define TRIO_POSITIVE_ZERO FP_POS_ZERO
596# define TRIO_NEGATIVE_ZERO FP_NEG_ZERO
597# define TRIO_POSITIVE_NORMAL FP_POS_NORM
598# define TRIO_NEGATIVE_NORMAL FP_NEG_NORM
599
600# elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
Daniel Veillard21458c82002-03-27 16:12:22 +0000601 /*
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000602 * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
603 * function.
Daniel Veillard21458c82002-03-27 16:12:22 +0000604 */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000605# define TRIO_FPCLASSIFY(n) _fpclass(n)
606# define TRIO_QUIET_NAN _FPCLASS_QNAN
607# define TRIO_SIGNALLING_NAN _FPCLASS_SNAN
608# define TRIO_POSITIVE_INFINITY _FPCLASS_PINF
609# define TRIO_NEGATIVE_INFINITY _FPCLASS_NINF
610# define TRIO_POSITIVE_SUBNORMAL _FPCLASS_PD
611# define TRIO_NEGATIVE_SUBNORMAL _FPCLASS_ND
612# define TRIO_POSITIVE_ZERO _FPCLASS_PZ
613# define TRIO_NEGATIVE_ZERO _FPCLASS_NZ
614# define TRIO_POSITIVE_NORMAL _FPCLASS_PN
615# define TRIO_NEGATIVE_NORMAL _FPCLASS_NN
616
617# elif defined(FP_PLUS_NORM)
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000618 /*
619 * HP-UX 9.x and 10.x have an fpclassify() function, that is different
620 * from the C99 fpclassify() macro supported on HP-UX 11.x.
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000621 *
622 * AIX has class() for C, and _class() for C++, which returns the
623 * same values as the HP-UX fpclassify() function.
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000624 */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000625# if defined(TRIO_PLATFORM_AIX)
626# if defined(__cplusplus)
627# define TRIO_FPCLASSIFY(n) _class(n)
628# else
629# define TRIO_FPCLASSIFY(n) class(n)
630# endif
631# else
632# define TRIO_FPCLASSIFY(n) fpclassify(n)
633# endif
634# define TRIO_QUIET_NAN FP_QNAN
635# define TRIO_SIGNALLING_NAN FP_SNAN
636# define TRIO_POSITIVE_INFINITY FP_PLUS_INF
637# define TRIO_NEGATIVE_INFINITY FP_MINUS_INF
638# define TRIO_POSITIVE_SUBNORMAL FP_PLUS_DENORM
639# define TRIO_NEGATIVE_SUBNORMAL FP_MINUS_DENORM
640# define TRIO_POSITIVE_ZERO FP_PLUS_ZERO
641# define TRIO_NEGATIVE_ZERO FP_MINUS_ZERO
642# define TRIO_POSITIVE_NORMAL FP_PLUS_NORM
643# define TRIO_NEGATIVE_NORMAL FP_MINUS_NORM
644# endif
645
646# if defined(TRIO_FPCLASSIFY)
647 switch (TRIO_FPCLASSIFY(number)) {
648 case TRIO_QUIET_NAN:
649 case TRIO_SIGNALLING_NAN:
650 *is_negative = TRIO_FALSE; /* NaN has no sign */
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000651 return TRIO_FP_NAN;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000652 case TRIO_POSITIVE_INFINITY:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000653 *is_negative = TRIO_FALSE;
654 return TRIO_FP_INFINITE;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000655 case TRIO_NEGATIVE_INFINITY:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000656 *is_negative = TRIO_TRUE;
657 return TRIO_FP_INFINITE;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000658 case TRIO_POSITIVE_SUBNORMAL:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000659 *is_negative = TRIO_FALSE;
660 return TRIO_FP_SUBNORMAL;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000661 case TRIO_NEGATIVE_SUBNORMAL:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000662 *is_negative = TRIO_TRUE;
663 return TRIO_FP_SUBNORMAL;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000664 case TRIO_POSITIVE_ZERO:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000665 *is_negative = TRIO_FALSE;
666 return TRIO_FP_ZERO;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000667 case TRIO_NEGATIVE_ZERO:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000668 *is_negative = TRIO_TRUE;
669 return TRIO_FP_ZERO;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000670 case TRIO_POSITIVE_NORMAL:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000671 *is_negative = TRIO_FALSE;
672 return TRIO_FP_NORMAL;
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000673 case TRIO_NEGATIVE_NORMAL:
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000674 *is_negative = TRIO_TRUE;
675 return TRIO_FP_NORMAL;
676 default:
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000677 /* Just in case... */
678 *is_negative = (number < 0.0);
679 return TRIO_FP_NORMAL;
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000680 }
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000681
682# else
Daniel Veillard21458c82002-03-27 16:12:22 +0000683 /*
684 * Fallback solution.
685 */
686 int rc;
687
688 if (number == 0.0) {
689 /*
690 * In IEEE 754 the sign of zero is ignored in comparisons, so we
691 * have to handle this as a special case by examining the sign bit
692 * directly.
693 */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000694# if defined(USE_IEEE_754)
Daniel Veillard21458c82002-03-27 16:12:22 +0000695 *is_negative = trio_is_negative(number);
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000696# else
Daniel Veillard21458c82002-03-27 16:12:22 +0000697 *is_negative = TRIO_FALSE; /* FIXME */
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000698# endif
Daniel Veillard21458c82002-03-27 16:12:22 +0000699 return TRIO_FP_ZERO;
700 }
701 if (trio_isnan(number)) {
702 *is_negative = TRIO_FALSE;
703 return TRIO_FP_NAN;
704 }
705 if ((rc = trio_isinf(number))) {
706 *is_negative = (rc == -1);
707 return TRIO_FP_INFINITE;
708 }
709 if ((number > 0.0) && (number < DBL_MIN)) {
710 *is_negative = TRIO_FALSE;
711 return TRIO_FP_SUBNORMAL;
712 }
713 if ((number < 0.0) && (number > -DBL_MIN)) {
714 *is_negative = TRIO_TRUE;
715 return TRIO_FP_SUBNORMAL;
716 }
717 *is_negative = (number < 0.0);
718 return TRIO_FP_NORMAL;
719
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000720# endif
Daniel Veillard21458c82002-03-27 16:12:22 +0000721#endif
722}
723
724/**
725 Examine the sign of a number.
726
727 @param number An arbitrary floating-point number.
728 @return Boolean value indicating whether or not the number has the
729 sign bit set (i.e. is negative).
730*/
731TRIO_PUBLIC int
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000732trio_signbit
733TRIO_ARGS1((number),
734 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000735{
736 int is_negative;
737
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000738 (void)trio_fpclassify_and_signbit(number, &is_negative);
Daniel Veillard21458c82002-03-27 16:12:22 +0000739 return is_negative;
740}
741
William M. Bracka71a8ef2003-08-06 04:43:55 +0000742#if 0
743 /* Temporary fix - this routine is not used in libxml */
Daniel Veillard21458c82002-03-27 16:12:22 +0000744/**
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
William M. Bracka71a8ef2003-08-06 04:43:55 +0000760#endif
Daniel Veillard21458c82002-03-27 16:12:22 +0000761
Bjorn Reese026d29f2002-01-19 15:40:18 +0000762/** @} SpecialQuantities */
763
Bjorn Reese45029602001-08-21 09:23:53 +0000764/*************************************************************************
Bjorn Reese026d29f2002-01-19 15:40:18 +0000765 * For test purposes.
766 *
767 * Add the following compiler option to include this test code.
768 *
769 * Unix : -DSTANDALONE
770 * VMS : /DEFINE=(STANDALONE)
Bjorn Reese45029602001-08-21 09:23:53 +0000771 */
772#if defined(STANDALONE)
773# include <stdio.h>
774
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000775static TRIO_CONST char *
776getClassification
Daniel Veillarde645e8c2002-10-22 17:35:37 +0000777TRIO_ARGS1((type),
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000778 int type)
Daniel Veillard21458c82002-03-27 16:12:22 +0000779{
780 switch (type) {
781 case TRIO_FP_INFINITE:
782 return "FP_INFINITE";
783 case TRIO_FP_NAN:
784 return "FP_NAN";
785 case TRIO_FP_NORMAL:
786 return "FP_NORMAL";
787 case TRIO_FP_SUBNORMAL:
788 return "FP_SUBNORMAL";
789 case TRIO_FP_ZERO:
790 return "FP_ZERO";
791 default:
792 return "FP_UNKNOWN";
793 }
794}
795
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000796static void
797print_class
Daniel Veillarde645e8c2002-10-22 17:35:37 +0000798TRIO_ARGS2((prefix, number),
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000799 TRIO_CONST char *prefix,
800 double number)
Daniel Veillard21458c82002-03-27 16:12:22 +0000801{
802 printf("%-6s: %s %-15s %g\n",
803 prefix,
804 trio_signbit(number) ? "-" : "+",
Daniel Veillard9339b742003-10-15 08:18:00 +0000805 getClassification(TRIO_FPCLASSIFY(number)),
Daniel Veillard21458c82002-03-27 16:12:22 +0000806 number);
807}
808
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000809int main(TRIO_NOARGS)
Bjorn Reese45029602001-08-21 09:23:53 +0000810{
811 double my_nan;
812 double my_pinf;
813 double my_ninf;
814# if defined(TRIO_PLATFORM_UNIX)
Daniel Veillardb7c29c32002-09-25 22:44:43 +0000815 void (*signal_handler) TRIO_PROTO((int));
Bjorn Reese45029602001-08-21 09:23:53 +0000816# endif
817
818 my_nan = trio_nan();
819 my_pinf = trio_pinf();
820 my_ninf = trio_ninf();
821
Daniel Veillard21458c82002-03-27 16:12:22 +0000822 print_class("Nan", my_nan);
823 print_class("PInf", my_pinf);
824 print_class("NInf", my_ninf);
825 print_class("PZero", 0.0);
826 print_class("NZero", -0.0);
827 print_class("PNorm", 1.0);
828 print_class("NNorm", -1.0);
829 print_class("PSub", 1.01e-307 - 1.00e-307);
830 print_class("NSub", 1.00e-307 - 1.01e-307);
831
Bjorn Reese45029602001-08-21 09:23:53 +0000832 printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
833 my_nan,
834 ((unsigned char *)&my_nan)[0],
835 ((unsigned char *)&my_nan)[1],
836 ((unsigned char *)&my_nan)[2],
837 ((unsigned char *)&my_nan)[3],
838 ((unsigned char *)&my_nan)[4],
839 ((unsigned char *)&my_nan)[5],
840 ((unsigned char *)&my_nan)[6],
841 ((unsigned char *)&my_nan)[7],
842 trio_isnan(my_nan), trio_isinf(my_nan));
843 printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
844 my_pinf,
845 ((unsigned char *)&my_pinf)[0],
846 ((unsigned char *)&my_pinf)[1],
847 ((unsigned char *)&my_pinf)[2],
848 ((unsigned char *)&my_pinf)[3],
849 ((unsigned char *)&my_pinf)[4],
850 ((unsigned char *)&my_pinf)[5],
851 ((unsigned char *)&my_pinf)[6],
852 ((unsigned char *)&my_pinf)[7],
853 trio_isnan(my_pinf), trio_isinf(my_pinf));
854 printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
855 my_ninf,
856 ((unsigned char *)&my_ninf)[0],
857 ((unsigned char *)&my_ninf)[1],
858 ((unsigned char *)&my_ninf)[2],
859 ((unsigned char *)&my_ninf)[3],
860 ((unsigned char *)&my_ninf)[4],
861 ((unsigned char *)&my_ninf)[5],
862 ((unsigned char *)&my_ninf)[6],
863 ((unsigned char *)&my_ninf)[7],
864 trio_isnan(my_ninf), trio_isinf(my_ninf));
865
866# if defined(TRIO_PLATFORM_UNIX)
867 signal_handler = signal(SIGFPE, SIG_IGN);
868# endif
869
870 my_pinf = DBL_MAX + DBL_MAX;
871 my_ninf = -my_pinf;
872 my_nan = my_pinf / my_pinf;
873
874# if defined(TRIO_PLATFORM_UNIX)
875 signal(SIGFPE, signal_handler);
876# endif
877
878 printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
879 my_nan,
880 ((unsigned char *)&my_nan)[0],
881 ((unsigned char *)&my_nan)[1],
882 ((unsigned char *)&my_nan)[2],
883 ((unsigned char *)&my_nan)[3],
884 ((unsigned char *)&my_nan)[4],
885 ((unsigned char *)&my_nan)[5],
886 ((unsigned char *)&my_nan)[6],
887 ((unsigned char *)&my_nan)[7],
888 trio_isnan(my_nan), trio_isinf(my_nan));
889 printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
890 my_pinf,
891 ((unsigned char *)&my_pinf)[0],
892 ((unsigned char *)&my_pinf)[1],
893 ((unsigned char *)&my_pinf)[2],
894 ((unsigned char *)&my_pinf)[3],
895 ((unsigned char *)&my_pinf)[4],
896 ((unsigned char *)&my_pinf)[5],
897 ((unsigned char *)&my_pinf)[6],
898 ((unsigned char *)&my_pinf)[7],
899 trio_isnan(my_pinf), trio_isinf(my_pinf));
900 printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
901 my_ninf,
902 ((unsigned char *)&my_ninf)[0],
903 ((unsigned char *)&my_ninf)[1],
904 ((unsigned char *)&my_ninf)[2],
905 ((unsigned char *)&my_ninf)[3],
906 ((unsigned char *)&my_ninf)[4],
907 ((unsigned char *)&my_ninf)[5],
908 ((unsigned char *)&my_ninf)[6],
909 ((unsigned char *)&my_ninf)[7],
910 trio_isnan(my_ninf), trio_isinf(my_ninf));
Daniel Veillarda48ed3d2003-04-03 15:28:28 +0000911
Bjorn Reese45029602001-08-21 09:23:53 +0000912 return 0;
913}
914#endif