1 | /*************************************************************************
|
---|
2 | *
|
---|
3 | * $Id$
|
---|
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 | * MERCHANTABILITY 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 (amongst 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$";
|
---|
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 | * endianness by the compiler, which means that we can support any
|
---|
127 | * kind of endianness. 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 | #if (defined(__BORLANDC__) && __BORLANDC__ >= 0x0590)
|
---|
133 | static TRIO_CONST double internalEndianMagic = 7.949928895127362e-275;
|
---|
134 | #else
|
---|
135 | static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;
|
---|
136 | #endif
|
---|
137 |
|
---|
138 | /* Mask for the exponent */
|
---|
139 | static TRIO_CONST unsigned char ieee_754_exponent_mask[] = {
|
---|
140 | 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
141 | };
|
---|
142 |
|
---|
143 | /* Mask for the mantissa */
|
---|
144 | static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = {
|
---|
145 | 0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
|
---|
146 | };
|
---|
147 |
|
---|
148 | /* Mask for the sign bit */
|
---|
149 | static TRIO_CONST unsigned char ieee_754_sign_mask[] = {
|
---|
150 | 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
151 | };
|
---|
152 |
|
---|
153 | /* Bit-pattern for negative zero */
|
---|
154 | static TRIO_CONST unsigned char ieee_754_negzero_array[] = {
|
---|
155 | 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
156 | };
|
---|
157 |
|
---|
158 | /* Bit-pattern for infinity */
|
---|
159 | static TRIO_CONST unsigned char ieee_754_infinity_array[] = {
|
---|
160 | 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
161 | };
|
---|
162 |
|
---|
163 | /* Bit-pattern for quiet NaN */
|
---|
164 | static TRIO_CONST unsigned char ieee_754_qnan_array[] = {
|
---|
165 | 0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
|
---|
166 | };
|
---|
167 |
|
---|
168 |
|
---|
169 | /*************************************************************************
|
---|
170 | * Functions
|
---|
171 | */
|
---|
172 |
|
---|
173 | /*
|
---|
174 | * trio_make_double
|
---|
175 | */
|
---|
176 | TRIO_PRIVATE double
|
---|
177 | trio_make_double
|
---|
178 | TRIO_ARGS1((values),
|
---|
179 | TRIO_CONST unsigned char *values)
|
---|
180 | {
|
---|
181 | TRIO_VOLATILE double result;
|
---|
182 | int i;
|
---|
183 |
|
---|
184 | for (i = 0; i < (int)sizeof(double); i++) {
|
---|
185 | ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i];
|
---|
186 | }
|
---|
187 | return result;
|
---|
188 | }
|
---|
189 |
|
---|
190 | /*
|
---|
191 | * trio_is_special_quantity
|
---|
192 | */
|
---|
193 | TRIO_PRIVATE int
|
---|
194 | trio_is_special_quantity
|
---|
195 | TRIO_ARGS2((number, has_mantissa),
|
---|
196 | double number,
|
---|
197 | int *has_mantissa)
|
---|
198 | {
|
---|
199 | unsigned int i;
|
---|
200 | unsigned char current;
|
---|
201 | int is_special_quantity = TRIO_TRUE;
|
---|
202 |
|
---|
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 |
|
---|
214 | /*
|
---|
215 | * trio_is_negative
|
---|
216 | */
|
---|
217 | TRIO_PRIVATE int
|
---|
218 | trio_is_negative
|
---|
219 | TRIO_ARGS1((number),
|
---|
220 | double number)
|
---|
221 | {
|
---|
222 | unsigned int i;
|
---|
223 | int is_negative = TRIO_FALSE;
|
---|
224 |
|
---|
225 | for (i = 0; i < (unsigned int)sizeof(double); i++) {
|
---|
226 | is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]
|
---|
227 | & ieee_754_sign_mask[i]);
|
---|
228 | }
|
---|
229 | return is_negative;
|
---|
230 | }
|
---|
231 |
|
---|
232 | #endif /* USE_IEEE_754 */
|
---|
233 |
|
---|
234 |
|
---|
235 | /**
|
---|
236 | Generate negative zero.
|
---|
237 |
|
---|
238 | @return Floating-point representation of negative zero.
|
---|
239 | */
|
---|
240 | TRIO_PUBLIC double
|
---|
241 | trio_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 | /**
|
---|
253 | Generate positive infinity.
|
---|
254 |
|
---|
255 | @return Floating-point representation of positive infinity.
|
---|
256 | */
|
---|
257 | TRIO_PUBLIC double
|
---|
258 | trio_pinf(TRIO_NOARGS)
|
---|
259 | {
|
---|
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 |
|
---|
296 | /**
|
---|
297 | Generate negative infinity.
|
---|
298 |
|
---|
299 | @return Floating-point value of negative infinity.
|
---|
300 | */
|
---|
301 | TRIO_PUBLIC double
|
---|
302 | trio_ninf(TRIO_NOARGS)
|
---|
303 | {
|
---|
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 |
|
---|
317 | /**
|
---|
318 | Generate NaN.
|
---|
319 |
|
---|
320 | @return Floating-point representation of NaN.
|
---|
321 | */
|
---|
322 | TRIO_PUBLIC double
|
---|
323 | trio_nan(TRIO_NOARGS)
|
---|
324 | {
|
---|
325 | /* Cache the result */
|
---|
326 | static double result = 0.0;
|
---|
327 |
|
---|
328 | if (result == 0.0) {
|
---|
329 |
|
---|
330 | #if defined(TRIO_COMPILER_SUPPORTS_C99)
|
---|
331 | result = nan("");
|
---|
332 |
|
---|
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 |
|
---|
364 | /**
|
---|
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 | */
|
---|
370 | TRIO_PUBLIC int
|
---|
371 | trio_isnan
|
---|
372 | TRIO_ARGS1((number),
|
---|
373 | double number)
|
---|
374 | {
|
---|
375 | #if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \
|
---|
376 | || defined(TRIO_COMPILER_SUPPORTS_UNIX95)
|
---|
377 | /*
|
---|
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 |
|
---|
385 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
386 | /*
|
---|
387 | * Microsoft Visual C++ and Borland C++ Builder have an _isnan()
|
---|
388 | * function.
|
---|
389 | */
|
---|
390 | return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;
|
---|
391 |
|
---|
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 | */
|
---|
418 | ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) ||
|
---|
419 | /*
|
---|
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 |
|
---|
435 | /**
|
---|
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 | */
|
---|
441 | TRIO_PUBLIC int
|
---|
442 | trio_isinf
|
---|
443 | TRIO_ARGS1((number),
|
---|
444 | double number)
|
---|
445 | {
|
---|
446 | #if defined(TRIO_COMPILER_DECC) && !defined(__linux__)
|
---|
447 | /*
|
---|
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));
|
---|
454 |
|
---|
455 | #elif defined(isinf)
|
---|
456 | /*
|
---|
457 | * C99 defines isinf() as a macro.
|
---|
458 | */
|
---|
459 | return isinf(number)
|
---|
460 | ? ((number > 0.0) ? 1 : -1)
|
---|
461 | : 0;
|
---|
462 |
|
---|
463 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
464 | /*
|
---|
465 | * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
|
---|
466 | * function that can be used to detect infinity.
|
---|
467 | */
|
---|
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 |
|
---|
511 | #if 0
|
---|
512 | /* Temporary fix - this routine is not used anywhere */
|
---|
513 | /**
|
---|
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 | */
|
---|
519 | TRIO_PUBLIC int
|
---|
520 | trio_isfinite
|
---|
521 | TRIO_ARGS1((number),
|
---|
522 | double number)
|
---|
523 | {
|
---|
524 | #if defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isfinite)
|
---|
525 | /*
|
---|
526 | * C99 defines isfinite() as a macro.
|
---|
527 | */
|
---|
528 | return isfinite(number);
|
---|
529 |
|
---|
530 | #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB)
|
---|
531 | /*
|
---|
532 | * Microsoft Visual C++ and Borland C++ Builder use _finite().
|
---|
533 | */
|
---|
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 |
|
---|
554 | #endif
|
---|
555 |
|
---|
556 | /*
|
---|
557 | * The sign of NaN is always false
|
---|
558 | */
|
---|
559 | TRIO_PUBLIC int
|
---|
560 | trio_fpclassify_and_signbit
|
---|
561 | TRIO_ARGS2((number, is_negative),
|
---|
562 | double number,
|
---|
563 | int *is_negative)
|
---|
564 | {
|
---|
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 |
|
---|
583 | #else
|
---|
584 | # if defined(TRIO_COMPILER_DECC)
|
---|
585 | /*
|
---|
586 | * DECC has an fp_class() function.
|
---|
587 | */
|
---|
588 | # 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)
|
---|
601 | /*
|
---|
602 | * Microsoft Visual C++ and Borland C++ Builder have an _fpclass()
|
---|
603 | * function.
|
---|
604 | */
|
---|
605 | # 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)
|
---|
618 | /*
|
---|
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.
|
---|
621 | *
|
---|
622 | * AIX has class() for C, and _class() for C++, which returns the
|
---|
623 | * same values as the HP-UX fpclassify() function.
|
---|
624 | */
|
---|
625 | # 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 */
|
---|
651 | return TRIO_FP_NAN;
|
---|
652 | case TRIO_POSITIVE_INFINITY:
|
---|
653 | *is_negative = TRIO_FALSE;
|
---|
654 | return TRIO_FP_INFINITE;
|
---|
655 | case TRIO_NEGATIVE_INFINITY:
|
---|
656 | *is_negative = TRIO_TRUE;
|
---|
657 | return TRIO_FP_INFINITE;
|
---|
658 | case TRIO_POSITIVE_SUBNORMAL:
|
---|
659 | *is_negative = TRIO_FALSE;
|
---|
660 | return TRIO_FP_SUBNORMAL;
|
---|
661 | case TRIO_NEGATIVE_SUBNORMAL:
|
---|
662 | *is_negative = TRIO_TRUE;
|
---|
663 | return TRIO_FP_SUBNORMAL;
|
---|
664 | case TRIO_POSITIVE_ZERO:
|
---|
665 | *is_negative = TRIO_FALSE;
|
---|
666 | return TRIO_FP_ZERO;
|
---|
667 | case TRIO_NEGATIVE_ZERO:
|
---|
668 | *is_negative = TRIO_TRUE;
|
---|
669 | return TRIO_FP_ZERO;
|
---|
670 | case TRIO_POSITIVE_NORMAL:
|
---|
671 | *is_negative = TRIO_FALSE;
|
---|
672 | return TRIO_FP_NORMAL;
|
---|
673 | case TRIO_NEGATIVE_NORMAL:
|
---|
674 | *is_negative = TRIO_TRUE;
|
---|
675 | return TRIO_FP_NORMAL;
|
---|
676 | default:
|
---|
677 | /* Just in case... */
|
---|
678 | *is_negative = (number < 0.0);
|
---|
679 | return TRIO_FP_NORMAL;
|
---|
680 | }
|
---|
681 |
|
---|
682 | # else
|
---|
683 | /*
|
---|
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 | */
|
---|
694 | # if defined(USE_IEEE_754)
|
---|
695 | *is_negative = trio_is_negative(number);
|
---|
696 | # else
|
---|
697 | *is_negative = TRIO_FALSE; /* FIXME */
|
---|
698 | # endif
|
---|
699 | 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 |
|
---|
720 | # endif
|
---|
721 | #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 | */
|
---|
731 | TRIO_PUBLIC int
|
---|
732 | trio_signbit
|
---|
733 | TRIO_ARGS1((number),
|
---|
734 | double number)
|
---|
735 | {
|
---|
736 | int is_negative;
|
---|
737 |
|
---|
738 | (void)trio_fpclassify_and_signbit(number, &is_negative);
|
---|
739 | return is_negative;
|
---|
740 | }
|
---|
741 |
|
---|
742 | #if 0
|
---|
743 | /* Temporary fix - this routine is not used in libxml */
|
---|
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 | */
|
---|
750 | TRIO_PUBLIC int
|
---|
751 | trio_fpclassify
|
---|
752 | TRIO_ARGS1((number),
|
---|
753 | double number)
|
---|
754 | {
|
---|
755 | int dummy;
|
---|
756 |
|
---|
757 | return trio_fpclassify_and_signbit(number, &dummy);
|
---|
758 | }
|
---|
759 |
|
---|
760 | #endif
|
---|
761 |
|
---|
762 | /** @} SpecialQuantities */
|
---|
763 |
|
---|
764 | /*************************************************************************
|
---|
765 | * For test purposes.
|
---|
766 | *
|
---|
767 | * Add the following compiler option to include this test code.
|
---|
768 | *
|
---|
769 | * Unix : -DSTANDALONE
|
---|
770 | * VMS : /DEFINE=(STANDALONE)
|
---|
771 | */
|
---|
772 | #if defined(STANDALONE)
|
---|
773 | # include <stdio.h>
|
---|
774 |
|
---|
775 | static TRIO_CONST char *
|
---|
776 | getClassification
|
---|
777 | TRIO_ARGS1((type),
|
---|
778 | int type)
|
---|
779 | {
|
---|
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 |
|
---|
796 | static void
|
---|
797 | print_class
|
---|
798 | TRIO_ARGS2((prefix, number),
|
---|
799 | TRIO_CONST char *prefix,
|
---|
800 | double number)
|
---|
801 | {
|
---|
802 | printf("%-6s: %s %-15s %g\n",
|
---|
803 | prefix,
|
---|
804 | trio_signbit(number) ? "-" : "+",
|
---|
805 | getClassification(TRIO_FPCLASSIFY(number)),
|
---|
806 | number);
|
---|
807 | }
|
---|
808 |
|
---|
809 | int main(TRIO_NOARGS)
|
---|
810 | {
|
---|
811 | double my_nan;
|
---|
812 | double my_pinf;
|
---|
813 | double my_ninf;
|
---|
814 | # if defined(TRIO_PLATFORM_UNIX)
|
---|
815 | void (*signal_handler) TRIO_PROTO((int));
|
---|
816 | # endif
|
---|
817 |
|
---|
818 | my_nan = trio_nan();
|
---|
819 | my_pinf = trio_pinf();
|
---|
820 | my_ninf = trio_ninf();
|
---|
821 |
|
---|
822 | 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 |
|
---|
832 | 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));
|
---|
911 |
|
---|
912 | return 0;
|
---|
913 | }
|
---|
914 | #endif
|
---|