📄 trionan.c.svn-base
字号:
/************************************************************************* * * $Id: trionan.c,v 1.14 2003/10/15 08:17:58 veillard Exp $ * * Copyright (C) 2001 Bjorn Reese <breese@users.sourceforge.net> * * Permission to use, copy, modify, and distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF * MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE AUTHORS AND * CONTRIBUTORS ACCEPT NO RESPONSIBILITY IN ANY CONCEIVABLE MANNER. * ************************************************************************ * * Functions to handle special quantities in floating-point numbers * (that is, NaNs and infinity). They provide the capability to detect * and fabricate special quantities. * * Although written to be as portable as possible, it can never be * guaranteed to work on all platforms, as not all hardware supports * special quantities. * * The approach used here (approximately) is to: * * 1. Use C99 functionality when available. * 2. Use IEEE 754 bit-patterns if possible. * 3. Use platform-specific techniques. * ************************************************************************//* * TODO: * o Put all the magic into trio_fpclassify_and_signbit(), and use this from * trio_isnan() etc. *//************************************************************************* * Include files */#include "triodef.h"#include "trionan.h"#include <math.h>#include <string.h>#include <limits.h>#include <float.h>#if defined(TRIO_PLATFORM_UNIX)# include <signal.h>#endif#if defined(TRIO_COMPILER_DECC)# if defined(__linux__)# include <cpml.h># else# include <fp_class.h># endif#endif#include <assert.h>#if defined(TRIO_DOCUMENTATION)# include "doc/doc_nan.h"#endif/** @addtogroup SpecialQuantities @{*//************************************************************************* * Definitions */#define TRIO_TRUE (1 == 1)#define TRIO_FALSE (0 == 1)/* * We must enable IEEE floating-point on Alpha */#if defined(__alpha) && !defined(_IEEE_FP)# if defined(TRIO_COMPILER_DECC)# if defined(TRIO_PLATFORM_VMS)# error "Must be compiled with option /IEEE_MODE=UNDERFLOW_TO_ZERO/FLOAT=IEEE"# else# if !defined(_CFE)# error "Must be compiled with option -ieee"# endif# endif# elif defined(TRIO_COMPILER_GCC) && (defined(__osf__) || defined(__linux__))# error "Must be compiled with option -mieee"# endif#endif /* __alpha && ! _IEEE_FP *//* * In ANSI/IEEE 754-1985 64-bits double format numbers have the * following properties (amoungst others) * * o FLT_RADIX == 2: binary encoding * o DBL_MAX_EXP == 1024: 11 bits exponent, where one bit is used * to indicate special numbers (e.g. NaN and Infinity), so the * maximum exponent is 10 bits wide (2^10 == 1024). * o DBL_MANT_DIG == 53: The mantissa is 52 bits wide, but because * numbers are normalized the initial binary 1 is represented * implicitly (the so-called "hidden bit"), which leaves us with * the ability to represent 53 bits wide mantissa. */#if (FLT_RADIX == 2) && (DBL_MAX_EXP == 1024) && (DBL_MANT_DIG == 53)# define USE_IEEE_754#endif/************************************************************************* * Constants */static TRIO_CONST char rcsid[] = "@(#)$Id: trionan.c,v 1.14 2003/10/15 08:17:58 veillard Exp $";#if defined(USE_IEEE_754)/* * Endian-agnostic indexing macro. * * The value of internalEndianMagic, when converted into a 64-bit * integer, becomes 0x0706050403020100 (we could have used a 64-bit * integer value instead of a double, but not all platforms supports * that type). The value is automatically encoded with the correct * endianess by the compiler, which means that we can support any * kind of endianess. The individual bytes are then used as an index * for the IEEE 754 bit-patterns and masks. */#define TRIO_DOUBLE_INDEX(x) (((unsigned char *)&internalEndianMagic)[7-(x)])static TRIO_CONST double internalEndianMagic = 7.949928895127363e-275;/* Mask for the exponent */static TRIO_CONST unsigned char ieee_754_exponent_mask[] = { 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};/* Mask for the mantissa */static TRIO_CONST unsigned char ieee_754_mantissa_mask[] = { 0x00, 0x0F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF};/* Mask for the sign bit */static TRIO_CONST unsigned char ieee_754_sign_mask[] = { 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};/* Bit-pattern for negative zero */static TRIO_CONST unsigned char ieee_754_negzero_array[] = { 0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};/* Bit-pattern for infinity */static TRIO_CONST unsigned char ieee_754_infinity_array[] = { 0x7F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};/* Bit-pattern for quiet NaN */static TRIO_CONST unsigned char ieee_754_qnan_array[] = { 0x7F, 0xF8, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};/************************************************************************* * Functions *//* * trio_make_double */TRIO_PRIVATE doubletrio_make_doubleTRIO_ARGS1((values), TRIO_CONST unsigned char *values){ TRIO_VOLATILE double result; int i; for (i = 0; i < (int)sizeof(double); i++) { ((TRIO_VOLATILE unsigned char *)&result)[TRIO_DOUBLE_INDEX(i)] = values[i]; } return result;}/* * trio_is_special_quantity */TRIO_PRIVATE inttrio_is_special_quantityTRIO_ARGS2((number, has_mantissa), double number, int *has_mantissa){ unsigned int i; unsigned char current; int is_special_quantity = TRIO_TRUE; *has_mantissa = 0; for (i = 0; i < (unsigned int)sizeof(double); i++) { current = ((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)]; is_special_quantity &= ((current & ieee_754_exponent_mask[i]) == ieee_754_exponent_mask[i]); *has_mantissa |= (current & ieee_754_mantissa_mask[i]); } return is_special_quantity;}/* * trio_is_negative */TRIO_PRIVATE inttrio_is_negativeTRIO_ARGS1((number), double number){ unsigned int i; int is_negative = TRIO_FALSE; for (i = 0; i < (unsigned int)sizeof(double); i++) { is_negative |= (((unsigned char *)&number)[TRIO_DOUBLE_INDEX(i)] & ieee_754_sign_mask[i]); } return is_negative;}#endif /* USE_IEEE_754 *//** Generate negative zero. @return Floating-point representation of negative zero.*/TRIO_PUBLIC doubletrio_nzero(TRIO_NOARGS){#if defined(USE_IEEE_754) return trio_make_double(ieee_754_negzero_array);#else TRIO_VOLATILE double zero = 0.0; return -zero;#endif}/** Generate positive infinity. @return Floating-point representation of positive infinity.*/TRIO_PUBLIC doubletrio_pinf(TRIO_NOARGS){ /* Cache the result */ static double result = 0.0; if (result == 0.0) { #if defined(INFINITY) && defined(__STDC_IEC_559__) result = (double)INFINITY;#elif defined(USE_IEEE_754) result = trio_make_double(ieee_754_infinity_array);#else /* * If HUGE_VAL is different from DBL_MAX, then HUGE_VAL is used * as infinity. Otherwise we have to resort to an overflow * operation to generate infinity. */# if defined(TRIO_PLATFORM_UNIX) void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);# endif result = HUGE_VAL; if (HUGE_VAL == DBL_MAX) { /* Force overflow */ result += HUGE_VAL; } # if defined(TRIO_PLATFORM_UNIX) signal(SIGFPE, signal_handler);# endif#endif } return result;}/** Generate negative infinity. @return Floating-point value of negative infinity.*/TRIO_PUBLIC doubletrio_ninf(TRIO_NOARGS){ static double result = 0.0; if (result == 0.0) { /* * Negative infinity is calculated by negating positive infinity, * which can be done because it is legal to do calculations on * infinity (for example, 1 / infinity == 0). */ result = -trio_pinf(); } return result;}/** Generate NaN. @return Floating-point representation of NaN.*/TRIO_PUBLIC doubletrio_nan(TRIO_NOARGS){ /* Cache the result */ static double result = 0.0; if (result == 0.0) { #if defined(TRIO_COMPILER_SUPPORTS_C99) result = nan("");#elif defined(NAN) && defined(__STDC_IEC_559__) result = (double)NAN; #elif defined(USE_IEEE_754) result = trio_make_double(ieee_754_qnan_array);#else /* * There are several ways to generate NaN. The one used here is * to divide infinity by infinity. I would have preferred to add * negative infinity to positive infinity, but that yields wrong * result (infinity) on FreeBSD. * * This may fail if the hardware does not support NaN, or if * the Invalid Operation floating-point exception is unmasked. */# if defined(TRIO_PLATFORM_UNIX) void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);# endif result = trio_pinf() / trio_pinf(); # if defined(TRIO_PLATFORM_UNIX) signal(SIGFPE, signal_handler);# endif #endif } return result;}/** Check for NaN. @param number An arbitrary floating-point number. @return Boolean value indicating whether or not the number is a NaN.*/TRIO_PUBLIC inttrio_isnanTRIO_ARGS1((number), double number){#if (defined(TRIO_COMPILER_SUPPORTS_C99) && defined(isnan)) \ || defined(TRIO_COMPILER_SUPPORTS_UNIX95) /* * C99 defines isnan() as a macro. UNIX95 defines isnan() as a * function. This function was already present in XPG4, but this * is a bit tricky to detect with compiler defines, so we choose * the conservative approach and only use it for UNIX95. */ return isnan(number); #elif defined(TRIO_COMPILER_MSVC) || defined(TRIO_COMPILER_BCB) /* * Microsoft Visual C++ and Borland C++ Builder have an _isnan() * function. */ return _isnan(number) ? TRIO_TRUE : TRIO_FALSE;#elif defined(USE_IEEE_754) /* * Examine IEEE 754 bit-pattern. A NaN must have a special exponent * pattern, and a non-empty mantissa. */ int has_mantissa; int is_special_quantity; is_special_quantity = trio_is_special_quantity(number, &has_mantissa); return (is_special_quantity && has_mantissa); #else /* * Fallback solution */ int status; double integral, fraction; # if defined(TRIO_PLATFORM_UNIX) void (*signal_handler)(int) = signal(SIGFPE, SIG_IGN);# endif status = (/* * NaN is the only number which does not compare to itself */ ((TRIO_VOLATILE double)number != (TRIO_VOLATILE double)number) || /* * Fallback solution if NaN compares to NaN */ ((number != 0.0) && (fraction = modf(number, &integral), integral == fraction))); # if defined(TRIO_PLATFORM_UNIX) signal(SIGFPE, signal_handler);# endif return status; #endif}/** Check for infinity. @param number An arbitrary floating-point number. @return 1 if positive infinity, -1 if negative infinity, 0 otherwise.*/TRIO_PUBLIC inttrio_isinfTRIO_ARGS1((number), double number){#if defined(TRIO_COMPILER_DECC) && !defined(__linux__) /* * DECC has an isinf() macro, but it works differently than that * of C99, so we use the fp_class() function instead. */ return ((fp_class(number) == FP_POS_INF) ? 1 : ((fp_class(number) == FP_NEG_INF) ? -1 : 0));#elif defined(isinf) /* * C99 defines isinf() as a macro. */ return isinf(number) ? ((number > 0.0) ? 1 : -1)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -