1 | /*************************************************************************
|
---|
2 | *
|
---|
3 | * $Id: trionan.c 2219 2003-10-15 08:18:00Z veillard $
|
---|
4 | *
|
---|
5 | * Copyright (C) 2001 Bjorn Reese <[email protected]>
|
---|
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 | * TODO:
|
---|
36 | * o Put all the magic into trio_fpclassify_and_signbit(), and use this from
|
---|
37 | * trio_isnan() etc.
|
---|
38 | */
|
---|
39 |
|
---|
40 | /*************************************************************************
|
---|
41 | * Include files
|
---|
42 | */
|
---|
43 | #include "triodef.h"
|
---|
44 | #include "trionan.h"
|
---|
45 |
|
---|
46 | #include <math.h>
|
---|
47 | #include <string.h>
|
---|
48 | #include <limits.h>
|
---|
49 | #include <float.h>
|
---|
50 | #if defined(TRIO_PLATFORM_UNIX)
|
---|
51 | # include <signal.h>
|
---|
52 | #endif
|
---|
53 | #if defined(TRIO_COMPILER_DECC)
|
---|
54 | # if defined(__linux__)
|
---|
55 | # include <cpml.h>
|
---|
56 | # else
|
---|
57 | # include <fp_class.h>
|
---|
58 | # endif
|
---|
59 | #endif
|
---|
60 | #include <assert.h>
|
---|
61 |
|
---|
62 | #if defined(TRIO_DOCUMENTATION)
|
---|
63 | # include "doc/doc_nan.h"
|
---|
64 | #endif
|
---|
65 | /** @addtogroup SpecialQuantities
|
---|
66 | @{
|
---|
67 | */
|
---|
68 |
|
---|
69 | /*************************************************************************
|
---|
70 | * Definitions
|
---|
71 | */
|
---|
72 |
|
---|
73 | #define TRIO_TRUE (1 == 1)
|
---|
74 | #define TRIO_FALSE (0 == 1)
|
---|
75 |
|
---|
76 | /*
|
---|
77 | * We must enable IEEE floating-point on Alpha
|
---|
78 | */
|
---|
79 | #if defined(__alpha) && !defined(_IEEE_FP)
|
---|
80 | # if defined(TRIO_COMPILER_DECC)
|
---|
81 | # if defined(TRIO_PLATFORM_VMS)
|
---|
82 | # error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"
|
---|
83 | # else
|
---|
84 | # if !defined(_CFE)
|
---|
85 | # error "Must be compiled with option -ieee"
|
---|
86 | # endif
|
---|
87 | # endif
|
---|
88 | # elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))
|
---|
89 | # error "Must be compiled with option -mieee"
|
---|
90 | # endif
|
---|
91 | #endif /* __alpha && ! _IEEE_FP */
|
---|
92 |
|
---|
93 | /*
|
---|
94 | * In ANSI/IEEE 754-1985 64-bits double format numbers have the
|
---|
95 | * following properties (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
|
---|
103 | * implicitly (the so-called "hidden bit"), which leaves us with
|
---|
104 | * the ability to represent 53 bits wide mantissa.
|
---|
105 | */
|
---|
106 | #if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)
|
---|
107 | # define USE_IEEE_754
|
---|
108 | #endif
|
---|
109 |
|
---|
110 |
|
---|
111 | /*************************************************************************
|
---|
112 | * Constants
|
---|
113 | */
|
---|
114 |
|
---|
115 | static TRIO_CONST char rcsid[] = "@(#)$Id: trionan.c 2219 2003-10-15 08:18:00Z veillard $";
|
---|
116 |
|
---|
117 | #if defined(USE_IEEE_754)
|
---|
118 |
|
---|
119 | /*
|
---|
120 | * Endian-agnostic indexing macro.
|
---|
121 | *
|
---|
122 | * The value of internalEndianMagic, when converted into a 64-bit
|
---|
123 | * integer, becomes 0x0706050403020100 (we could have used a 64-bit
|
---|
124 | * integer value instead of a double, but not all platforms supports
|
---|
125 | * that type). The value is automatically encoded with the correct
|
---|
126 | * 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 | */
|
---|
130 | #define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])
|
---|
131 |
|
---|
132 | static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
|
---|
133 |
|
---|
134 | /* Mask for the exponent */
|
---|
135 | static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
|
---|
136 | 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
137 | };
|
---|
138 |
|
---|
139 | /* Mask for the mantissa */
|
---|
140 | static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
|
---|
141 | 0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
|
---|
142 | };
|
---|
143 |
|
---|
144 | /* Mask for the sign bit */
|
---|
145 | static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
|
---|
146 | 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
147 | };
|
---|
148 |
|
---|
149 | /* Bit-pattern for negative zero */
|
---|
150 | static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
|
---|
151 | 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
152 | };
|
---|
153 |
|
---|
154 | /* Bit-pattern for infinity */
|
---|
155 | static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
|
---|
156 | 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
157 | };
|
---|
158 |
|
---|
159 | /* Bit-pattern for quiet NaN */
|
---|
160 | static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
|
---|
161 | 0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
162 | };
|
---|
163 |
|
---|
164 |
|
---|
165 | /*************************************************************************
|
---|
166 | * Functions
|
---|
167 | */
|
---|
168 |
|
---|
169 | /*
|
---|
170 | * trio_make_double
|
---|
171 | */
|
---|
172 | TRIO_PRIVATE double
|
---|
173 | trio_make_double
|
---|
174 | TRIO_ARGS1((values),
|
---|
175 | TRIO_CONST unsigned char *values)
|
---|
176 | {
|
---|
177 | TRIO_VOLATILE double result;
|
---|
178 | int i;
|
---|
179 |
|
---|
180 | for (i = 0; i < (int)sizeof(double); i++) {
|
---|
181 | ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
|
---|
182 | }
|
---|
183 | return result;
|
---|
184 | }
|
---|
185 |
|
---|
186 | /*
|
---|
187 | * trio_is_special_quantity
|
---|
188 | */
|
---|
189 | TRIO_PRIVATE int
|
---|
190 | trio_is_special_quantity
|
---|
191 | TRIO_ARGS2((number, has_mantissa),
|
---|
192 | double number,
|
---|
193 | int *has_mantissa)
|
---|
194 | {
|
---|
195 | unsigned int i;
|
---|
196 | unsigned char current;
|
---|
197 | int is_special_quantity = TRIO_TRUE;
|
---|
198 |
|
---|
199 | *has_mantissa = 0;
|
---|
200 |
|
---|
201 | for (i = 0; i < (unsigned int)sizeof(double); i++) {
|
---|
202 | current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)];
|
---|
203 | is_special_quantity
|
---|
204 | &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]);
|
---|
205 | *has_mantissa |= (current & ieee_754_mantissa_mask[i]);
|
---|
206 | }
|
---|
207 | return is_special_quantity;
|
---|
208 | }
|
---|
209 |
|
---|
210 | /*
|
---|
211 | * trio_is_negative
|
---|
212 | */
|
---|
213 | TRIO_PRIVATE int
|
---|
214 | trio_is_negative
|
---|
215 | TRIO_ARGS1((number),
|
---|
216 | double number)
|
---|
217 | {
|
---|
218 | unsigned int i;
|
---|
219 | int is_negative = TRIO_FALSE;
|
---|
220 |
|
---|
221 | for (i = 0; i < (unsigned int)sizeof(double); i++) {
|
---|
222 | is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
|
---|
223 | & ieee_754_sign_mask[i]);
|
---|
224 | }
|
---|
225 | return is_negative;
|
---|
226 | }
|
---|
227 |
|
---|
228 | #endif /* USE_IEEE_754 */
|
---|
229 |
|
---|
230 |
|
---|
231 | /**
|
---|
232 | Generate negative zero.
|
---|
233 |
|
---|
234 | @return Floating-point representation of negative zero.
|
---|
235 | */
|
---|
236 | TRIO_PUBLIC double
|
---|
237 | trio_nzero(TRIO_NOARGS)
|
---|
238 | {
|
---|
239 | #if defined(USE_IEEE_754)
|
---|
240 | return trio_make_double(ieee_754_negzero_array);
|
---|
241 | #else
|
---|
242 | TRIO_VOLATILE double zero = 0.0;
|
---|
243 |
|
---|
244 | return -zero;
|
---|
245 | #endif
|
---|
246 | }
|
---|
247 |
|
---|
248 | /**
|
---|
249 | Generate positive infinity.
|
---|
250 |
|
---|
251 | @return Floating-point representation of positive infinity.
|
---|
252 | */
|
---|
253 | TRIO_PUBLIC double
|
---|
254 | trio_pinf(TRIO_NOARGS)
|
---|
255 | {
|
---|
256 | /* Cache the result */
|
---|
257 | static double result = 0.0;
|
---|
258 |
|
---|
259 | if (result == 0.0) {
|
---|
260 |
|
---|
261 | #if defined(INFINITY) && defined(__STDC_IEC_559__)
|
---|
262 | result = (double)INFINITY;
|
---|
263 |
|
---|
264 | #elif defined(USE_IEEE_754)
|
---|
265 | result = trio_make_double(ieee_754_infinity_array);
|
---|
266 |
|
---|
267 | #else
|
---|
268 | /*
|
---|
269 | * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used
|
---|
270 | * as infinity. Otherwise we have to resort to an overflow
|
---|
271 | * operation to generate infinity.
|
---|
272 | */
|
---|
273 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
274 | void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
|
---|
275 | # endif
|
---|
276 |
|
---|
277 | result = HUGE_VAL;
|
---|
278 | if (HUGE_VAL == DBL_MAX) {
|
---|
279 | /* Force overflow */
|
---|
280 | result += HUGE_VAL;
|
---|
281 | }
|
---|
282 |
|
---|
283 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
284 | signal(SIGFPE, signal_handler);
|
---|
285 | # endif
|
---|
286 |
|
---|
287 | #endif
|
---|
288 | }
|
---|
289 | return result;
|
---|
290 | }
|
---|
291 |
|
---|
292 | /**
|
---|
293 | Generate negative infinity.
|
---|
294 |
|
---|
295 | @return Floating-point value of negative infinity.
|
---|
296 | */
|
---|
297 | TRIO_PUBLIC double
|
---|
298 | trio_ninf(TRIO_NOARGS)
|
---|
299 | {
|
---|
300 | static double result = 0.0;
|
---|
301 |
|
---|
302 | if (result == 0.0) {
|
---|
303 | /*
|
---|
304 | * Negative infinity is calculated by negating positive infinity,
|
---|
305 | * which can be done because it is legal to do calculations on
|
---|
306 | * infinity (for example, 1 / infinity == 0).
|
---|
307 | */
|
---|
308 | result = -trio_pinf();
|
---|
309 | }
|
---|
310 | return result;
|
---|
311 | }
|
---|
312 |
|
---|
313 | /**
|
---|
314 | Generate NaN.
|
---|
315 |
|
---|
316 | @return Floating-point representation of NaN.
|
---|
317 | */
|
---|
318 | TRIO_PUBLIC double
|
---|
319 | trio_nan(TRIO_NOARGS)
|
---|
320 | {
|
---|
321 | /* Cache the result */
|
---|
322 | static double result = 0.0;
|
---|
323 |
|
---|
324 | if (result == 0.0) {
|
---|
325 |
|
---|
326 | #if defined(TRIO_COMPILER_SUPPORTS_C99)
|
---|
327 | result = nan("");
|
---|
328 |
|
---|
329 | #elif defined(NAN) && defined(__STDC_IEC_559__)
|
---|
330 | result = (double)NAN;
|
---|
331 |
|
---|
332 | #elif defined(USE_IEEE_754)
|
---|
333 | result = trio_make_double(ieee_754_qnan_array);
|
---|
334 |
|
---|
335 | #else
|
---|
336 | /*
|
---|
337 | * There are several ways to generate NaN. The one used here is
|
---|
338 | * to divide infinity by infinity. I would have preferred to add
|
---|
339 | * negative infinity to positive infinity, but that yields wrong
|
---|
340 | * result (infinity) on FreeBSD.
|
---|
341 | *
|
---|
342 | * This may fail if the hardware does not support NaN, or if
|
---|
343 | * the Invalid Operation floating-point exception is unmasked.
|
---|
344 | */
|
---|
345 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
346 | void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
|
---|
347 | # endif
|
---|
348 |
|
---|
349 | result = trio_pinf() / trio_pinf();
|
---|
350 |
|
---|
351 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
352 | signal(SIGFPE, signal_handler);
|
---|
353 | # endif
|
---|
354 |
|
---|
355 | #endif
|
---|
356 | }
|
---|
357 | return result;
|
---|
358 | }
|
---|
359 |
|
---|
360 | /**
|
---|
361 | Check for NaN.
|
---|
362 |
|
---|
363 | @param number An arbitrary floating-point number.
|
---|
364 | @return Boolean value indicating whether or not the number is a NaN.
|
---|
365 | */
|
---|
366 | TRIO_PUBLIC int
|
---|
367 | trio_isnan
|
---|
368 | TRIO_ARGS1((number),
|
---|
369 | double number)
|
---|
370 | {
|
---|
371 | #if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \
|
---|
372 | || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
|
---|
373 | /*
|
---|
374 | * C99 defines isnan() as a macro. UNIX95 defines isnan() as a
|
---|
375 | * function. This function was already present in XPG4, but this
|
---|
376 | * is a bit tricky to detect with compiler defines, so we choose
|
---|
377 | * the conservative approach and only use it for UNIX95.
|
---|
378 | */
|
---|
379 | return isnan(number);
|
---|
380 |
|
---|
381 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
382 | /*
|
---|
383 | * Microsoft Visual C++ and Borland C++ Builder have an _isnan()
|
---|
384 | * function.
|
---|
385 | */
|
---|
386 | return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;
|
---|
387 |
|
---|
388 | #elif defined(USE_IEEE_754)
|
---|
389 | /*
|
---|
390 | * Examine IEEE 754 bit-pattern. A NaN must have a special exponent
|
---|
391 | * pattern, and a non-empty mantissa.
|
---|
392 | */
|
---|
393 | int has_mantissa;
|
---|
394 | int is_special_quantity;
|
---|
395 |
|
---|
396 | is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
|
---|
397 |
|
---|
398 | return (is_special_quantity && has_mantissa);
|
---|
399 |
|
---|
400 | #else
|
---|
401 | /*
|
---|
402 | * Fallback solution
|
---|
403 | */
|
---|
404 | int status;
|
---|
405 | double integral, fraction;
|
---|
406 |
|
---|
407 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
408 | void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
|
---|
409 | # endif
|
---|
410 |
|
---|
411 | status = (/*
|
---|
412 | * NaN is the only number which does not compare to itself
|
---|
413 | */
|
---|
414 | ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
|
---|
415 | /*
|
---|
416 | * Fallback solution if NaN compares to NaN
|
---|
417 | */
|
---|
418 | ((number != 0.0) &&
|
---|
419 | (fraction = modf(number, &integral),
|
---|
420 | integral == fraction)));
|
---|
421 |
|
---|
422 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
423 | signal(SIGFPE, signal_handler);
|
---|
424 | # endif
|
---|
425 |
|
---|
426 | return status;
|
---|
427 |
|
---|
428 | #endif
|
---|
429 | }
|
---|
430 |
|
---|
431 | /**
|
---|
432 | Check for infinity.
|
---|
433 |
|
---|
434 | @param number An arbitrary floating-point number.
|
---|
435 | @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.
|
---|
436 | */
|
---|
437 | TRIO_PUBLIC int
|
---|
438 | trio_isinf
|
---|
439 | TRIO_ARGS1((number),
|
---|
440 | double number)
|
---|
441 | {
|
---|
442 | #if defined(TRIO_COMPILER_DECC) && !defined(__linux__)
|
---|
443 | /*
|
---|
444 | * DECC has an isinf() macro, but it works differently than that
|
---|
445 | * of C99, so we use the fp_class() function instead.
|
---|
446 | */
|
---|
447 | return ((fp_class(number) == FP_POS_INF)
|
---|
448 | ? 1
|
---|
449 | : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));
|
---|
450 |
|
---|
451 | #elif defined(isinf)
|
---|
452 | /*
|
---|
453 | * C99 defines isinf() as a macro.
|
---|
454 | */
|
---|
455 | return isinf(number)
|
---|
456 | ? ((number > 0.0) ? 1 : -1)
|
---|
457 | : 0;
|
---|
458 |
|
---|
459 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
460 | /*
|
---|
461 | * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
|
---|
462 | * function that can be used to detect infinity.
|
---|
463 | */
|
---|
464 | return ((_fpclass(number) == _FPCLASS_PINF)
|
---|
465 | ? 1
|
---|
466 | : ((_fpclass(number) == _FPCLASS_NINF) ? -1 : 0));
|
---|
467 |
|
---|
468 | #elif defined(USE_IEEE_754)
|
---|
469 | /*
|
---|
470 | * Examine IEEE 754 bit-pattern. Infinity must have a special exponent
|
---|
471 | * pattern, and an empty mantissa.
|
---|
472 | */
|
---|
473 | int has_mantissa;
|
---|
474 | int is_special_quantity;
|
---|
475 |
|
---|
476 | is_special_quantity = trio_is_special_quantity(number, &has_mantissa);
|
---|
477 |
|
---|
478 | return (is_special_quantity && !has_mantissa)
|
---|
479 | ? ((number < 0.0) ? -1 : 1)
|
---|
480 | : 0;
|
---|
481 |
|
---|
482 | #else
|
---|
483 | /*
|
---|
484 | * Fallback solution.
|
---|
485 | */
|
---|
486 | int status;
|
---|
487 |
|
---|
488 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
489 | void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);
|
---|
490 | # endif
|
---|
491 |
|
---|
492 | double infinity = trio_pinf();
|
---|
493 |
|
---|
494 | status = ((number == infinity)
|
---|
495 | ? 1
|
---|
496 | : ((number == -infinity) ? -1 : 0));
|
---|
497 |
|
---|
498 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
499 | signal(SIGFPE, signal_handler);
|
---|
500 | # endif
|
---|
501 |
|
---|
502 | return status;
|
---|
503 |
|
---|
504 | #endif
|
---|
505 | }
|
---|
506 |
|
---|
507 | #if 0
|
---|
508 | /* Temporary fix - this routine is not used anywhere */
|
---|
509 | /**
|
---|
510 | Check for finity.
|
---|
511 |
|
---|
512 | @param number An arbitrary floating-point number.
|
---|
513 | @return Boolean value indicating whether or not the number is a finite.
|
---|
514 | */
|
---|
515 | TRIO_PUBLIC int
|
---|
516 | trio_isfinite
|
---|
517 | TRIO_ARGS1((number),
|
---|
518 | double number)
|
---|
519 | {
|
---|
520 | #if defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isfinite)
|
---|
521 | /*
|
---|
522 | * C99 defines isfinite() as a macro.
|
---|
523 | */
|
---|
524 | return isfinite(number);
|
---|
525 |
|
---|
526 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
527 | /*
|
---|
528 | * Microsoft Visual C++ and Borland C++ Builder use _finite().
|
---|
529 | */
|
---|
530 | return _finite(number);
|
---|
531 |
|
---|
532 | #elif defined(USE_IEEE_754)
|
---|
533 | /*
|
---|
534 | * Examine IEEE 754 bit-pattern. For finity we do not care about the
|
---|
535 | * mantissa.
|
---|
536 | */
|
---|
537 | int dummy;
|
---|
538 |
|
---|
539 | return (! trio_is_special_quantity(number, &dummy));
|
---|
540 |
|
---|
541 | #else
|
---|
542 | /*
|
---|
543 | * Fallback solution.
|
---|
544 | */
|
---|
545 | return ((trio_isinf(number) == 0) && (trio_isnan(number) == 0));
|
---|
546 |
|
---|
547 | #endif
|
---|
548 | }
|
---|
549 |
|
---|
550 | #endif
|
---|
551 |
|
---|
552 | /*
|
---|
553 | * The sign of NaN is always false
|
---|
554 | */
|
---|
555 | TRIO_PUBLIC int
|
---|
556 | trio_fpclassify_and_signbit
|
---|
557 | TRIO_ARGS2((number, is_negative),
|
---|
558 | double number,
|
---|
559 | int *is_negative)
|
---|
560 | {
|
---|
561 | #if defined(fpclassify) && defined(signbit)
|
---|
562 | /*
|
---|
563 | * C99 defines fpclassify() and signbit() as a macros
|
---|
564 | */
|
---|
565 | *is_negative = signbit(number);
|
---|
566 | switch (fpclassify(number)) {
|
---|
567 | case FP_NAN:
|
---|
568 | return TRIO_FP_NAN;
|
---|
569 | case FP_INFINITE:
|
---|
570 | return TRIO_FP_INFINITE;
|
---|
571 | case FP_SUBNORMAL:
|
---|
572 | return TRIO_FP_SUBNORMAL;
|
---|
573 | case FP_ZERO:
|
---|
574 | return TRIO_FP_ZERO;
|
---|
575 | default:
|
---|
576 | return TRIO_FP_NORMAL;
|
---|
577 | }
|
---|
578 |
|
---|
579 | #else
|
---|
580 | # if defined(TRIO_COMPILER_DECC)
|
---|
581 | /*
|
---|
582 | * DECC has an fp_class() function.
|
---|
583 | */
|
---|
584 | # define TRIO_FPCLASSIFY(n) fp_class(n)
|
---|
585 | # define TRIO_QUIET_NAN FP_QNAN
|
---|
586 | # define TRIO_SIGNALLING_NAN FP_SNAN
|
---|
587 | # define TRIO_POSITIVE_INFINITY FP_POS_INF
|
---|
588 | # define TRIO_NEGATIVE_INFINITY FP_NEG_INF
|
---|
589 | # define TRIO_POSITIVE_SUBNORMAL FP_POS_DENORM
|
---|
590 | # define TRIO_NEGATIVE_SUBNORMAL FP_NEG_DENORM
|
---|
591 | # define TRIO_POSITIVE_ZERO FP_POS_ZERO
|
---|
592 | # define TRIO_NEGATIVE_ZERO FP_NEG_ZERO
|
---|
593 | # define TRIO_POSITIVE_NORMAL FP_POS_NORM
|
---|
594 | # define TRIO_NEGATIVE_NORMAL FP_NEG_NORM
|
---|
595 |
|
---|
596 | # elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
597 | /*
|
---|
598 | * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
|
---|
599 | * function.
|
---|
600 | */
|
---|
601 | # define TRIO_FPCLASSIFY(n) _fpclass(n)
|
---|
602 | # define TRIO_QUIET_NAN _FPCLASS_QNAN
|
---|
603 | # define TRIO_SIGNALLING_NAN _FPCLASS_SNAN
|
---|
604 | # define TRIO_POSITIVE_INFINITY _FPCLASS_PINF
|
---|
605 | # define TRIO_NEGATIVE_INFINITY _FPCLASS_NINF
|
---|
606 | # define TRIO_POSITIVE_SUBNORMAL _FPCLASS_PD
|
---|
607 | # define TRIO_NEGATIVE_SUBNORMAL _FPCLASS_ND
|
---|
608 | # define TRIO_POSITIVE_ZERO _FPCLASS_PZ
|
---|
609 | # define TRIO_NEGATIVE_ZERO _FPCLASS_NZ
|
---|
610 | # define TRIO_POSITIVE_NORMAL _FPCLASS_PN
|
---|
611 | # define TRIO_NEGATIVE_NORMAL _FPCLASS_NN
|
---|
612 |
|
---|
613 | # elif defined(FP_PLUS_NORM)
|
---|
614 | /*
|
---|
615 | * HP-UX 9.x and 10.x have an fpclassify() function, that is different
|
---|
616 | * from the C99 fpclassify() macro supported on HP-UX 11.x.
|
---|
617 | *
|
---|
618 | * AIX has class() for C, and _class() for C++, which returns the
|
---|
619 | * same values as the HP-UX fpclassify() function.
|
---|
620 | */
|
---|
621 | # if defined(TRIO_PLATFORM_AIX)
|
---|
622 | # if defined(__cplusplus)
|
---|
623 | # define TRIO_FPCLASSIFY(n) _class(n)
|
---|
624 | # else
|
---|
625 | # define TRIO_FPCLASSIFY(n) class(n)
|
---|
626 | # endif
|
---|
627 | # else
|
---|
628 | # define TRIO_FPCLASSIFY(n) fpclassify(n)
|
---|
629 | # endif
|
---|
630 | # define TRIO_QUIET_NAN FP_QNAN
|
---|
631 | # define TRIO_SIGNALLING_NAN FP_SNAN
|
---|
632 | # define TRIO_POSITIVE_INFINITY FP_PLUS_INF
|
---|
633 | # define TRIO_NEGATIVE_INFINITY FP_MINUS_INF
|
---|
634 | # define TRIO_POSITIVE_SUBNORMAL FP_PLUS_DENORM
|
---|
635 | # define TRIO_NEGATIVE_SUBNORMAL FP_MINUS_DENORM
|
---|
636 | # define TRIO_POSITIVE_ZERO FP_PLUS_ZERO
|
---|
637 | # define TRIO_NEGATIVE_ZERO FP_MINUS_ZERO
|
---|
638 | # define TRIO_POSITIVE_NORMAL FP_PLUS_NORM
|
---|
639 | # define TRIO_NEGATIVE_NORMAL FP_MINUS_NORM
|
---|
640 | # endif
|
---|
641 |
|
---|
642 | # if defined(TRIO_FPCLASSIFY)
|
---|
643 | switch (TRIO_FPCLASSIFY(number)) {
|
---|
644 | case TRIO_QUIET_NAN:
|
---|
645 | case TRIO_SIGNALLING_NAN:
|
---|
646 | *is_negative = TRIO_FALSE; /* NaN has no sign */
|
---|
647 | return TRIO_FP_NAN;
|
---|
648 | case TRIO_POSITIVE_INFINITY:
|
---|
649 | *is_negative = TRIO_FALSE;
|
---|
650 | return TRIO_FP_INFINITE;
|
---|
651 | case TRIO_NEGATIVE_INFINITY:
|
---|
652 | *is_negative = TRIO_TRUE;
|
---|
653 | return TRIO_FP_INFINITE;
|
---|
654 | case TRIO_POSITIVE_SUBNORMAL:
|
---|
655 | *is_negative = TRIO_FALSE;
|
---|
656 | return TRIO_FP_SUBNORMAL;
|
---|
657 | case TRIO_NEGATIVE_SUBNORMAL:
|
---|
658 | *is_negative = TRIO_TRUE;
|
---|
659 | return TRIO_FP_SUBNORMAL;
|
---|
660 | case TRIO_POSITIVE_ZERO:
|
---|
661 | *is_negative = TRIO_FALSE;
|
---|
662 | return TRIO_FP_ZERO;
|
---|
663 | case TRIO_NEGATIVE_ZERO:
|
---|
664 | *is_negative = TRIO_TRUE;
|
---|
665 | return TRIO_FP_ZERO;
|
---|
666 | case TRIO_POSITIVE_NORMAL:
|
---|
667 | *is_negative = TRIO_FALSE;
|
---|
668 | return TRIO_FP_NORMAL;
|
---|
669 | case TRIO_NEGATIVE_NORMAL:
|
---|
670 | *is_negative = TRIO_TRUE;
|
---|
671 | return TRIO_FP_NORMAL;
|
---|
672 | default:
|
---|
673 | /* Just in case... */
|
---|
674 | *is_negative = (number < 0.0);
|
---|
675 | return TRIO_FP_NORMAL;
|
---|
676 | }
|
---|
677 |
|
---|
678 | # else
|
---|
679 | /*
|
---|
680 | * Fallback solution.
|
---|
681 | */
|
---|
682 | int rc;
|
---|
683 |
|
---|
684 | if (number == 0.0) {
|
---|
685 | /*
|
---|
686 | * In IEEE 754 the sign of zero is ignored in comparisons, so we
|
---|
687 | * have to handle this as a special case by examining the sign bit
|
---|
688 | * directly.
|
---|
689 | */
|
---|
690 | # if defined(USE_IEEE_754)
|
---|
691 | *is_negative = trio_is_negative(number);
|
---|
692 | # else
|
---|
693 | *is_negative = TRIO_FALSE; /* FIXME */
|
---|
694 | # endif
|
---|
695 | return TRIO_FP_ZERO;
|
---|
696 | }
|
---|
697 | if (trio_isnan(number)) {
|
---|
698 | *is_negative = TRIO_FALSE;
|
---|
699 | return TRIO_FP_NAN;
|
---|
700 | }
|
---|
701 | if ((rc = trio_isinf(number))) {
|
---|
702 | *is_negative = (rc == -1);
|
---|
703 | return TRIO_FP_INFINITE;
|
---|
704 | }
|
---|
705 | if ((number > 0.0) && (number < DBL_MIN)) {
|
---|
706 | *is_negative = TRIO_FALSE;
|
---|
707 | return TRIO_FP_SUBNORMAL;
|
---|
708 | }
|
---|
709 | if ((number < 0.0) && (number > -DBL_MIN)) {
|
---|
710 | *is_negative = TRIO_TRUE;
|
---|
711 | return TRIO_FP_SUBNORMAL;
|
---|
712 | }
|
---|
713 | *is_negative = (number < 0.0);
|
---|
714 | return TRIO_FP_NORMAL;
|
---|
715 |
|
---|
716 | # endif
|
---|
717 | #endif
|
---|
718 | }
|
---|
719 |
|
---|
720 | /**
|
---|
721 | Examine the sign of a number.
|
---|
722 |
|
---|
723 | @param number An arbitrary floating-point number.
|
---|
724 | @return Boolean value indicating whether or not the number has the
|
---|
725 | sign bit set (i.e. is negative).
|
---|
726 | */
|
---|
727 | TRIO_PUBLIC int
|
---|
728 | trio_signbit
|
---|
729 | TRIO_ARGS1((number),
|
---|
730 | double number)
|
---|
731 | {
|
---|
732 | int is_negative;
|
---|
733 |
|
---|
734 | (void)trio_fpclassify_and_signbit(number, &is_negative);
|
---|
735 | return is_negative;
|
---|
736 | }
|
---|
737 |
|
---|
738 | #if 0
|
---|
739 | /* Temporary fix - this routine is not used in libxml */
|
---|
740 | /**
|
---|
741 | Examine the class of a number.
|
---|
742 |
|
---|
743 | @param number An arbitrary floating-point number.
|
---|
744 | @return Enumerable value indicating the class of @p number
|
---|
745 | */
|
---|
746 | TRIO_PUBLIC int
|
---|
747 | trio_fpclassify
|
---|
748 | TRIO_ARGS1((number),
|
---|
749 | double number)
|
---|
750 | {
|
---|
751 | int dummy;
|
---|
752 |
|
---|
753 | return trio_fpclassify_and_signbit(number, &dummy);
|
---|
754 | }
|
---|
755 |
|
---|
756 | #endif
|
---|
757 |
|
---|
758 | /** @} SpecialQuantities */
|
---|
759 |
|
---|
760 | /*************************************************************************
|
---|
761 | * For test purposes.
|
---|
762 | *
|
---|
763 | * Add the following compiler option to include this test code.
|
---|
764 | *
|
---|
765 | * Unix : -DSTANDALONE
|
---|
766 | * VMS : /DEFINE=(STANDALONE)
|
---|
767 | */
|
---|
768 | #if defined(STANDALONE)
|
---|
769 | # include <stdio.h>
|
---|
770 |
|
---|
771 | static TRIO_CONST char *
|
---|
772 | getClassification
|
---|
773 | TRIO_ARGS1((type),
|
---|
774 | int type)
|
---|
775 | {
|
---|
776 | switch (type) {
|
---|
777 | case TRIO_FP_INFINITE:
|
---|
778 | return "FP_INFINITE";
|
---|
779 | case TRIO_FP_NAN:
|
---|
780 | return "FP_NAN";
|
---|
781 | case TRIO_FP_NORMAL:
|
---|
782 | return "FP_NORMAL";
|
---|
783 | case TRIO_FP_SUBNORMAL:
|
---|
784 | return "FP_SUBNORMAL";
|
---|
785 | case TRIO_FP_ZERO:
|
---|
786 | return "FP_ZERO";
|
---|
787 | default:
|
---|
788 | return "FP_UNKNOWN";
|
---|
789 | }
|
---|
790 | }
|
---|
791 |
|
---|
792 | static void
|
---|
793 | print_class
|
---|
794 | TRIO_ARGS2((prefix, number),
|
---|
795 | TRIO_CONST char *prefix,
|
---|
796 | double number)
|
---|
797 | {
|
---|
798 | printf("%-6s: %s %-15s %g\n",
|
---|
799 | prefix,
|
---|
800 | trio_signbit(number) ? "-" : "+",
|
---|
801 | getClassification(TRIO_FPCLASSIFY(number)),
|
---|
802 | number);
|
---|
803 | }
|
---|
804 |
|
---|
805 | int main(TRIO_NOARGS)
|
---|
806 | {
|
---|
807 | double my_nan;
|
---|
808 | double my_pinf;
|
---|
809 | double my_ninf;
|
---|
810 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
811 | void (*signal_handler) TRIO_PROTO((int));
|
---|
812 | # endif
|
---|
813 |
|
---|
814 | my_nan = trio_nan();
|
---|
815 | my_pinf = trio_pinf();
|
---|
816 | my_ninf = trio_ninf();
|
---|
817 |
|
---|
818 | print_class("Nan", my_nan);
|
---|
819 | print_class("PInf", my_pinf);
|
---|
820 | print_class("NInf", my_ninf);
|
---|
821 | print_class("PZero", 0.0);
|
---|
822 | print_class("NZero", -0.0);
|
---|
823 | print_class("PNorm", 1.0);
|
---|
824 | print_class("NNorm", -1.0);
|
---|
825 | print_class("PSub", 1.01e-307 - 1.00e-307);
|
---|
826 | print_class("NSub", 1.00e-307 - 1.01e-307);
|
---|
827 |
|
---|
828 | printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
829 | my_nan,
|
---|
830 | ((unsigned char *)&my_nan)[0],
|
---|
831 | ((unsigned char *)&my_nan)[1],
|
---|
832 | ((unsigned char *)&my_nan)[2],
|
---|
833 | ((unsigned char *)&my_nan)[3],
|
---|
834 | ((unsigned char *)&my_nan)[4],
|
---|
835 | ((unsigned char *)&my_nan)[5],
|
---|
836 | ((unsigned char *)&my_nan)[6],
|
---|
837 | ((unsigned char *)&my_nan)[7],
|
---|
838 | trio_isnan(my_nan), trio_isinf(my_nan));
|
---|
839 | printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
840 | my_pinf,
|
---|
841 | ((unsigned char *)&my_pinf)[0],
|
---|
842 | ((unsigned char *)&my_pinf)[1],
|
---|
843 | ((unsigned char *)&my_pinf)[2],
|
---|
844 | ((unsigned char *)&my_pinf)[3],
|
---|
845 | ((unsigned char *)&my_pinf)[4],
|
---|
846 | ((unsigned char *)&my_pinf)[5],
|
---|
847 | ((unsigned char *)&my_pinf)[6],
|
---|
848 | ((unsigned char *)&my_pinf)[7],
|
---|
849 | trio_isnan(my_pinf), trio_isinf(my_pinf));
|
---|
850 | printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
851 | my_ninf,
|
---|
852 | ((unsigned char *)&my_ninf)[0],
|
---|
853 | ((unsigned char *)&my_ninf)[1],
|
---|
854 | ((unsigned char *)&my_ninf)[2],
|
---|
855 | ((unsigned char *)&my_ninf)[3],
|
---|
856 | ((unsigned char *)&my_ninf)[4],
|
---|
857 | ((unsigned char *)&my_ninf)[5],
|
---|
858 | ((unsigned char *)&my_ninf)[6],
|
---|
859 | ((unsigned char *)&my_ninf)[7],
|
---|
860 | trio_isnan(my_ninf), trio_isinf(my_ninf));
|
---|
861 |
|
---|
862 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
863 | signal_handler = signal(SIGFPE, SIG_IGN);
|
---|
864 | # endif
|
---|
865 |
|
---|
866 | my_pinf = DBL_MAX + DBL_MAX;
|
---|
867 | my_ninf = -my_pinf;
|
---|
868 | my_nan = my_pinf / my_pinf;
|
---|
869 |
|
---|
870 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
871 | signal(SIGFPE, signal_handler);
|
---|
872 | # endif
|
---|
873 |
|
---|
874 | printf("NaN : %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
875 | my_nan,
|
---|
876 | ((unsigned char *)&my_nan)[0],
|
---|
877 | ((unsigned char *)&my_nan)[1],
|
---|
878 | ((unsigned char *)&my_nan)[2],
|
---|
879 | ((unsigned char *)&my_nan)[3],
|
---|
880 | ((unsigned char *)&my_nan)[4],
|
---|
881 | ((unsigned char *)&my_nan)[5],
|
---|
882 | ((unsigned char *)&my_nan)[6],
|
---|
883 | ((unsigned char *)&my_nan)[7],
|
---|
884 | trio_isnan(my_nan), trio_isinf(my_nan));
|
---|
885 | printf("PInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
886 | my_pinf,
|
---|
887 | ((unsigned char *)&my_pinf)[0],
|
---|
888 | ((unsigned char *)&my_pinf)[1],
|
---|
889 | ((unsigned char *)&my_pinf)[2],
|
---|
890 | ((unsigned char *)&my_pinf)[3],
|
---|
891 | ((unsigned char *)&my_pinf)[4],
|
---|
892 | ((unsigned char *)&my_pinf)[5],
|
---|
893 | ((unsigned char *)&my_pinf)[6],
|
---|
894 | ((unsigned char *)&my_pinf)[7],
|
---|
895 | trio_isnan(my_pinf), trio_isinf(my_pinf));
|
---|
896 | printf("NInf: %4g 0x%02x%02x%02x%02x%02x%02x%02x%02x (%2d, %2d)\n",
|
---|
897 | my_ninf,
|
---|
898 | ((unsigned char *)&my_ninf)[0],
|
---|
899 | ((unsigned char *)&my_ninf)[1],
|
---|
900 | ((unsigned char *)&my_ninf)[2],
|
---|
901 | ((unsigned char *)&my_ninf)[3],
|
---|
902 | ((unsigned char *)&my_ninf)[4],
|
---|
903 | ((unsigned char *)&my_ninf)[5],
|
---|
904 | ((unsigned char *)&my_ninf)[6],
|
---|
905 | ((unsigned char *)&my_ninf)[7],
|
---|
906 | trio_isnan(my_ninf), trio_isinf(my_ninf));
|
---|
907 |
|
---|
908 | return 0;
|
---|
909 | }
|
---|
910 | #endif
|
---|