root/branches/newmole/source/cddefines.h

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef _CDDEFINES_H_
#define _CDDEFINES_H_

#ifdef _MSC_VER
/* disable warning that conditional expression is constant, true or false in if */
#       pragma warning( disable : 4127 )
/* we are not using MS foundation class */
#       ifndef WIN32_LEAN_AND_MEAN
#               define WIN32_LEAN_AND_MEAN
#       endif
#endif
/* these headers are needed by all files */
/*lint -e129 these resolve several issues pclint has with my system headers */
/*lint -e78 */
/*lint -e830 */
/*lint -e38 */
/*lint -e148 */
/*lint -e114 */
/*lint -e18 */
/*lint -e49 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <float.h>
#include <limits>
#include <limits.h>
#include <time.h>
#include <signal.h>
#include <string>
#include <sstream>
#include <vector>
#include <valarray>
#include <complex>
#include <map>
#include <memory>
#include <stdexcept>
#include <algorithm>
#include <fstream>
#include <bitset>
#ifdef DMALLOC
#include <dmalloc.h>
#endif
/*lint +e18 */
/*lint +e49 */
/*lint +e38 */
/*lint +e148 */
/*lint +e830 */
/*lint +e78 */
/*lint -e129 */

using namespace std;

/** define EXTERN as nothing before including cddefines.h to allocate global data */
#ifndef EXTERN
#define EXTERN extern 
#endif

#undef STATIC
#ifdef USE_GPROF
#define STATIC
#else
#define STATIC static
#endif

#ifdef FLT_IS_DBL
typedef double realnum;
#else
typedef float realnum;
#endif

//Compile-time assertion after Alexandrescu
template<bool> struct StaticAssertFailed;
template<> struct StaticAssertFailed<true> {};
#define STATIC_ASSERT(x) ((void)StaticAssertFailed< (x) == true >())

typedef float sys_float;
// prevent explicit float's from creeping back into the code
#define float PLEASE_USE_REALNUM_NOT_FLOAT

/* make sure this is globally visible as well! */
#include "cpu.h"

//*************************************************************************
//
//! Singleton: template to construct classes where a single instance can
//! be accessed across the entire program. Construct your class as follows:
//!
//! class t_myclass : public Singleton<t_myclass>
//! {
//!     friend class Singleton<t_myclass>;
//! protected:
//!     t_myclass(); // make sure the contructor is protected!
//! public:
//!     long myfunc() const { return 43; }
//! };
//!
//! and use as follows throughout the code:
//!
//! long l = t_myclass::Inst().myfunc();
//!
//!   NB NB - This implementation is not threadsafe !!
//
// This implementation has been obtained from Wikipedia
//
//*************************************************************************

template<typename T> class Singleton
{
public:
        static T& Inst()
        {
                static T instance;  // assumes T has a protected default constructor
                return instance;
        }
};

/**************************************************************************
 *
 * these are variables and pointers for output from the code, used everywhere
 * declared extern here, and definition is in cddefines.cpp
 *
 **************************************************************************/

/** ioQQQ is the file hand to the output file itself,
 * ioQQQ is set to stdout by default in cdInit, 
 * and is reset to anything else by calling cdOutp */
EXTERN FILE *ioQQQ;

EXTERN FILE *ioStdin;

extern FILE *ioMAP;

/** we shall write errors to this file, it is set
 * to stderr in cdInit */
EXTERN FILE* ioPrnErr;

/** this set true when abort sequence is initiated - serious meltdown is happening */
EXTERN bool lgAbort;

/** flag lgTestIt turned on if routine TestCode ever called, only generates
 * comment that test code is in place */
EXTERN bool lgTestCodeCalled; 

/** flag lgTestOn set true with SET TEST command
 * for some test code to be run somewhere */
EXTERN bool lgTestCodeEnabled;

/** this is flag saying whether to print errors to 
 * standard error output */
EXTERN bool lgPrnErr;

/** nzone is zone counter, incremented in routine cloudy 
 * is zero during search phase, 1 for first zone at illuminated face */
EXTERN long int nzone;

/** this is nzone + conv.nPres2Ioniz/100 in ConvBase */
EXTERN double fnzone;

/** the iteration counter, set and incremented in routine cloudy,
 * ==1 during first iteration, 2 during second, etc */
EXTERN long int iteration;

/** number of simulations done with this coreload - 1 for first sim */
EXTERN long int nSimThisCoreload;

/**
 * this is the number zero, used to trick clever compilers when 
 * dividing by it to crash program
 * there is a routine called zero - this name cannot overlap
 * definition is in cddefines.cpp
 */
extern const double ZeroNum;

/**************************************************************************
 *
 * these are constants used to dimension several vectors and index arrays
 *
 **************************************************************************/

/** FILENAME_PATH_LENGTH is the size of the string that holds the path.  The longest
 * string that can be held is one less than this, due to the end
 * end of string sentinel in C.  Increase this is a larger string
 * is needed to hold the path on your system */
const int FILENAME_PATH_LENGTH = 200; 

/** twice the above, so that we can add file name to end of full path */
const int FILENAME_PATH_LENGTH_2  = FILENAME_PATH_LENGTH*2; 

/** this is limit to longest line of information that is scanned in,
 * end of line char is actually at this +1, dim of vector is this + 1 
 * all routines that scan in information should use this to dim vars */ 
const int INPUT_LINE_LENGTH = 200;

/** This is the number of elements included in the code,
 * is used to set lengths of many vectors */
const int LIMELM = 30;

/** the number of iso sequences now in the code */
const int NISO = 2;

/** following is real limit to how many levels can be computed for
 * model hydrogen atom - this has the one extra included, so at most
 * iso.numLevels_max[ipH_LIKE] can be NHYDRO_MAX_LEVEL-1 and vectors should be dim NHYDRO_MAX_LEVEL*/
const int NHYDRO_MAX_LEVEL = 401;

/** this is the maximum particle density allowed in cm^-3 */
const double MAX_DENSITY = 1.e21;

#if 0
/** The default value of the stopping temperature */
const double TEMP_STOP_DEFAULT = 4000.;
/** highest temperature to ever allow */
const double TEMP_LIMIT_HIGH = 1e10;
/** lowest temperature to ever allow */
const double TEMP_LIMIT_LOW = 2.8;
#endif

/** this is used to add to depth to prevent div or log of zero */
const double DEPTH_OFFSET = 1.e-30;

 /* these are flags for various colliders that are used across the code */
enum collider {
        ipELECTRON,
        ipPROTON,
        ipHE_PLUS,
        ipALPHA,
        //ipATOMH,
        //ipHE_2PLUS,
        //ipH2_ORTHO,
        //ipH2_PARA,
        ipNCOLLIDER
};

/* indices within recombination coefficient array */
/* ipRecEsc is state specific escape probability*/
const int ipRecEsc = 2;
/* the net escaping, including destruction by background and optical deepth*/
const int ipRecNetEsc = 1;
/* ipRecRad is state specific radiative recombination rate*/
const int ipRecRad = 0;
/** with the above, the total radiative rec per ion is
 * iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
   iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc]*dense.eden; */

/* these specify the form of the line redistribution function */
/* partial redistribution with wings */
const int ipPRD = 1;
/* complete redistribution, core only, no wings, Hummer's K2 function */
const int ipCRD = -1;
/* complete redistribution with wings */
const int ipCRDW = 2;
/* redistribution function for Lya, calls Hummer routine for H-like series only */
const int ipLY_A = -2;
/* core function for K2 destruction */
const int ipDEST_K2 = 1;
/* core function for complete redist destruction */
const int ipDEST_INCOM = 2;
/* core function for simple destruction */
const int ipDEST_SIMPL = 3;

/** some levels for hydrogenic species */

const int ipH1s = 0;
const int ipH2s = 1;
const int ipH2p = 2;
const int ipH3s = 3;
const int ipH3p = 4;
const int ipH3d = 5;
const int ipH4s = 6;
const int ipH4p = 7;
const int ipH4d = 8;
const int ipH4f = 9;
const int ipH5s = 10;
const int ipH5p = 11;
const int ipH5d = 12;
const int ipH5f = 13;
const int ipH5g = 14;
const int ipH6s = 15;

/** some levels for he-like species */

/* level 1 */
const int ipHe1s1S = 0;

/* level 2 */
const int ipHe2s3S = 1;
const int ipHe2s1S = 2;
const int ipHe2p3P0 = 3;
const int ipHe2p3P1 = 4;
const int ipHe2p3P2 = 5;
const int ipHe2p1P = 6;

/* level 3 */
const int ipHe3s3S = 7;
const int ipHe3s1S = 8;
const int ipHe3p3P = 9;
const int ipHe3d3D = 10;
const int ipHe3d1D = 11;
const int ipHe3p1P = 12;

/* level 4 */
const int ipHe4s3S = 13;
const int ipHe4s1S = 14;
const int ipHe4p3P = 15;
const int ipHe4d3D = 16;
const int ipHe4d1D = 17;
const int ipHe4f3F = 18;
const int ipHe4f1F = 19;
const int ipHe4p1P = 20;

/** these are array indices for isoelectronic sequences,
 * same as element but used for array addressing to make
 * context totally clear */
const int ipH_LIKE = 0;
const int ipHE_LIKE = 1;
const int ipLI_LIKE = 2;
const int ipBE_LIKE = 3;
const int ipB_LIKE = 4;
const int ipC_LIKE = 5;
const int ipN_LIKE = 6;
const int ipO_LIKE = 7;
const int ipF_LIKE = 8;
const int ipNE_LIKE = 9;
const int ipNA_LIKE = 10;
const int ipMG_LIKE = 11;
const int ipAL_LIKE = 12;
const int ipSI_LIKE = 13;
const int ipP_LIKE = 14;
const int ipS_LIKE = 15;
const int ipCL_LIKE = 16;
const int ipAR_LIKE = 17;

/** these are indices for some elements, on the C scale */
const int ipHYDROGEN = 0;
const int ipHELIUM = 1;
const int ipLITHIUM = 2;
const int ipBERYLLIUM = 3;
const int ipBORON = 4;
const int ipCARBON = 5;
const int ipNITROGEN = 6;
const int ipOXYGEN = 7;
const int ipFLUORINE = 8;
const int ipNEON = 9;
const int ipSODIUM = 10;
const int ipMAGNESIUM = 11;
const int ipALUMINIUM = 12;
const int ipSILICON = 13;
const int ipPHOSPHORUS = 14;
const int ipSULPHUR = 15;
const int ipCHLORINE = 16;
const int ipARGON = 17;
const int ipPOTASSIUM = 18;
const int ipCALCIUM = 19;
const int ipSCANDIUM = 20;
const int ipTITANIUM = 21;
const int ipVANADIUM = 22;
const int ipCHROMIUM = 23;
const int ipMANGANESE = 24;
const int ipIRON = 25;
const int ipCOBALT = 26;
const int ipNICKEL = 27;
const int ipCOPPER = 28;
const int ipZINC = 29;

/***************************************************************************
 * the following are prototypes for some routines that are part of the
 * debugging process - they come and go in any particular sub.  
 * it is not necessary to declare them when used since they are defined here
 **************************************************************************/

/**
 fudge enter fudge factors, or some arbitrary number, with fudge command
 return value is the fudge factor 
 fudge(-1) queries the routine for the number of fudge parameters that 
 were entered, zero returned if none
 \param ipnt integer saying which of the possible numbers on the fudge
 command to use - 0 would be the first
 */ 
double fudge(long int ipnt);

/**
 broken set flag saying that the code is broken 
 */ 
void broken(void);

/**fixit set flag saying that this code needs attention, but is not broken,
 * code is in service.cpp */
void fixit(void);

/**CodeReview - placed next to code that needs to be checked */ 
void CodeReview(void);

/**TestCode set flag saying that test code is in place */
void TestCode(void);

/**MyMalloc wrapper for malloc().  Returns a good pointer or dies.
\param size use same type as library function malloc
\param file
\param line
*/
void *MyMalloc(size_t size, const char *file, int line);

/**MyCalloc wrapper for calloc().  Returns a good pointer or dies. 
\param num use same type as library function CALLOC
\param size
*/
void *MyCalloc(size_t num, size_t size);

/**MyRealloc wrapper for realloc().  Returns a good pointer or dies. 
\param num use same type as library function REALLOC
\param size
*/
void *MyRealloc(void *p, size_t size);

/** MyAssert a version of assert that fails gracefully
\param *file
\param line
*/ 
void MyAssert(const char *file, int line);

/** the Cloudy exit handler 
\param iexit, EXIT_FAILURE for failure, EXIT_SUCCESS for success 
*/
NORETURN void cdExit(int iexit);

class cloudy_exit
{
        const char* p_routine;
        const char* p_file;
        long p_line;
        int p_exit;
public:
        cloudy_exit(const char* routine, const char* file, long line, int exit_code)
        {
                p_routine = routine;
                p_file = file;
                p_line = line;
                p_exit = exit_code;
        }
        virtual ~cloudy_exit() throw()
        {
                p_routine = NULL;
                p_file = NULL;
        }
        const char* routine() const throw()
        {
                return p_routine;
        }
        const char* file() const throw()
        {
                return p_file;
        }
        long line() const
        {
                return p_line;
        }
        int exit_status() const
        {
                return p_exit;
        }
};

// workarounds for __func__ are defined in cpu.h
#define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL )

// calls like puts( "[Stop in MyRoutine]" ) have been integrated in cdEXIT above
#define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed

/** print comment asking to show output to me */
void ShowMe(void);

/**TotalInsanity general error handler for something that cannot happen, exits */
NORETURN void TotalInsanity(void);

/**BadRead tried to read internal data and failed */
NORETURN void BadRead(void);

/**BadOpen tried to open data files and failed */
NORETURN void BadOpen(void);

/**
   NoNumb general error handler for no numbers on input line 
   \param *chCard
 */ 
NORETURN void NoNumb(char *chCard);

/** dbg_printf is a debug print routine that was provided by Peter Teuben,
 * as a component from his NEMO package.  It offers run-time specification
 * of the level of debugging */ 
int dbg_printf(int debug, const char *fmt, ...);

/** dprintf -- version of fprintf which prepends DEBUG */
int dprintf(FILE *fp, const char *format, ...);

/**read_whole_line safe version of fgets - read a line, 
 * return null if cannot read line or if input line is too long 
 \param char *chLine - previously allocated string where the line
 image will be stored
 \param int nChar size of chLine, we will return NULL if input line is
 longer than this
 \param FILE *ioIN a previously opened file handle, will read from from here
*/
char *read_whole_line( char *chLine , int nChar , FILE *ioIN );

/**************************************************************************
 *
 * various macros used by the code
 *
 **************************************************************************/

/** to avoid errors introduced by C's infamous double-negative logic,
 * this uses NDEBUG (the ANSI std macro used to tell assert that
 * we are not debugging) to define DEBUG */
#ifndef NDEBUG
#       define DEBUG
#else
#       undef DEBUG
#endif

/** use special version of malloc - it tests result and dies if cannot allocate space */
#if defined(malloc)
/* ...but if malloc is a macro, assume it is instrumented by a memory debugging tool
 * (e.g. dmalloc) */
#       define MALLOC(exp) (malloc(exp))
#else
/* Otherwise instrument and protect it ourselves */
#       define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
#endif

/** now special version of calloc - it dies if cannot allocate space. */
#if defined(calloc)
/* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */
#       define CALLOC calloc
#else
/* Otherwise instrument and protect it ourselves */
#       define CALLOC MyCalloc 
#endif

/** now special version of calloc - it dies if cannot allocate space. */
#if defined(realloc)
/* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */
#       define REALLOC realloc
#else
/* Otherwise instrument and protect it ourselves */
#       define REALLOC MyRealloc 
#endif

class bad_signal
{
        int p_sig;
public:
        explicit bad_signal(int sig)
        {
                p_sig = sig;
        }
        virtual ~bad_signal() throw() {}
        int sig() const throw()
        {
                return p_sig;
        }
};

class bad_assert
{
        const char* p_file;
        long p_line;
public:
        bad_assert(const char* file, long line)
        {
                p_file = file;
                p_line = line;
        }
        virtual ~bad_assert() throw()
        {
                p_file = NULL;
        }
        const char* file() const throw()
        {
                return p_file;
        }
        long line() const throw()
        {
                return p_line;
        }
};

class bad_mpi {
        long p_failMode;
public:
        bad_mpi(long failMode) : p_failMode(failMode) {}
        long failMode() const { return p_failMode; }
};

/* the do { ... } while ( 0 ) construct prevents bugs in code like this:
 * if( test )
 *      ASSERT( n == 10 );
 * else
 *      do something else...
 */
#undef  ASSERT
#ifndef OLD_ASSERT
/** a cloudy version of the assert macro that can be recovered from */
#       ifdef NDEBUG
#               define ASSERT(exp) ((void)0)
#       else
#               define ASSERT(exp) \
                        do { \
                                if (!(exp)) \
                                { \
                                        if( cpu.lgAssertAbort() ) \
                                                abort();          \
                                        else \
                                                throw bad_assert(__FILE__,__LINE__); \
                                } \
                        } while( 0 )
#       endif
#else
/** the old version of the assert macro that terminates safely */
#       ifdef NDEBUG
#               define ASSERT(exp) ((void)0)
#       else
#               define ASSERT(exp) \
                        do { \
                                if (!(exp)) \
                                        MyAssert(__FILE__, __LINE__); \
                        } while( 0 )
#       endif
#endif

NORETURN inline void OUT_OF_RANGE(const char* str)
{
        if( cpu.lgAssertAbort() )
                abort();
        else
                throw out_of_range( str );
}

/* Windows does not define isnan */
/* use our version on all platforms since the isnanf
 * function does not exist under Solaris 9 either */
#undef isnan
#define isnan MyIsnan

/** entry and exit of each routine will go here,
 * macros enabled if compiler-set flag DEBUG_FUN is defined */
class t_debug : public Singleton<t_debug>
{
        friend class Singleton<t_debug>;
        FILE *p_fp;
        int p_callLevel;
protected:
        t_debug()
        {
                p_fp = stderr;
                p_callLevel = 0;
        }
public:
        void enter(const char *name)
        {
                ++p_callLevel;
                fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name);
        }
        void leave(const char *name)
        {
                fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name);
                --p_callLevel;
        }               
};

template<bool lgTrace>
class debugtrace
{
        const char *p_name;
public:
        explicit debugtrace(const char *funcname)
        {
                p_name = funcname;
#               ifdef _MSC_VER
                /* disable warning that conditional expression is constant, true or false in if */
#               pragma warning( disable : 4127 )
#               endif
                if( lgTrace )
                        t_debug::Inst().enter(p_name);
        }
        ~debugtrace()
        {
#               ifdef _MSC_VER
                /* disable warning that conditional expression is constant, true or false in if */
#               pragma warning( disable : 4127 )
#               endif
                if( lgTrace )
                        t_debug::Inst().leave(p_name);
                p_name = NULL;
        }
        const char* name() const
        {
                return p_name;
        }
};

#ifdef DEBUG_FUN
#define DEBUG_ENTRY( funcname ) debugtrace<true> DEBUG_ENTRY( funcname )
#else
#ifdef __GNUC__
#define DEBUG_ENTRY( funcname ) ((void)0)
#else
#define DEBUG_ENTRY( funcname ) debugtrace<false> DEBUG_ENTRY( funcname )
#endif
#endif

/* TorF(l) returns a 'T' or 'F' depending on the 'logical' expr 'l' */
inline char TorF( bool l ) { return l ? 'T' : 'F'; }
/* */

/** checks whether argument is odd */
inline bool is_odd( int j ) { return (j&1) == 1; }
inline bool is_odd( long j ) { return (j&1L) == 1L; }
/* */

/** nint rounds to the nearest long int */
inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
/* */

/* define min for mixed arguments, the rest already exists */
inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); }

#undef MIN2
/** MIN2 takes two arguments, returns the smaller of the two */
#define MIN2 min
/* */

#undef MIN3
/** MIN3 takes 3 arguments, returns the smallest of the 3 */
#define MIN3(a,b,c) (min(min(a,b),c))
/* */

#undef MIN4
/** MIN2 takes 4 arguments, returns the smallest of the 4 */
#define MIN4(a,b,c,d) (min(min(a,b),min(c,d)))
/* */

/* define max for mixed arguments, the rest already exists */
inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); }

#undef MAX2
/** MAX2 takes two arguments, returns the larger of the two */
#define MAX2 max
/* */

#undef MAX3
/** MAX3 takes 3 arguments, returns the largest of the 3 */
#define MAX3(a,b,c) (max(max(a,b),c))
/* */

#undef MAX4
/** MAX2 takes 4 arguments, returns the largest of the 4 */
#define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
/* */

/** FP sign transfer (fortran sign function) - sign of y times abs value of x 
\param x
\param y
*/
template<class T>
inline T sign( T x, T y )
{
        return ( y < T() ) ? -abs(x) : abs(x);
}
/* */

/** sign3 returns -1 for negative arguments, +1 for positive, and 0 for zero (pascal sign function) */
template<class T>
inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
/* */

/** checks whether two FP numbers are "equal" (differ no more than n epsilon) */
inline bool fp_equal( sys_float x, sys_float y, int n=3 )
{
#ifdef _MSC_VER
        /* disable warning that conditional expression is constant, true or false in if */
#       pragma warning( disable : 4127 )
        /* we are not using MS foundation class */
#       ifndef WIN32_LEAN_AND_MEAN
#               define WIN32_LEAN_AND_MEAN
#       endif
#endif
        ASSERT( n >= 1 );
        // mimic IEEE behavior
        if( isnan(x) || isnan(y) )
                return false;
        int sx = sign3(x);
        int sy = sign3(y);
        // treat zero cases first to avoid division by zero below
        if( sx == 0 && sy == 0 )
                return true;
        // either x or y is zero (but not both), or x and y have different sign
        if( sx*sy != 1 )
                return false;
        x = abs(x);
        y = abs(y);
        return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*FLT_EPSILON );
}

inline bool fp_equal( double x, double y, int n=3 )
{
        ASSERT( n >= 1 );
        // mimic IEEE behavior
        if( isnan(x) || isnan(y) )
                return false;
        int sx = sign3(x);
        int sy = sign3(y);
        // treat zero cases first to avoid division by zero below
        if( sx == 0 && sy == 0 )
                return true;
        // either x or y is zero (but not both), or x and y have different sign
        if( sx*sy != 1 )
                return false;
        x = abs(x);
        y = abs(y);
        return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*DBL_EPSILON );
}

inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol )
{
        ASSERT( tol > 0.f );
        // mimic IEEE behavior
        if( isnan(x) || isnan(y) )
                return false;
        // make sure the tolerance is not too stringent
        ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) );
        return ( abs( x-y ) <= tol );
}

inline bool fp_equal_tol( double x, double y, double tol )
{
        ASSERT( tol > 0. );
        // mimic IEEE behavior
        if( isnan(x) || isnan(y) )
                return false;
        // make sure the tolerance is not too stringent
        ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) );
        return ( abs( x-y ) <= tol );
}

/** checks whether a number is within bounds */
inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
{
        ASSERT( n >= 1 );
        // mimic IEEE behavior
        if( isnan(x) || isnan(lo) || isnan(hi) )
                return false;
        if( fp_equal(lo,hi,n) )
                return fp_equal(0.5f*(lo+hi),x,n);
        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
                return false;
        return true;
}
inline bool fp_bound( double lo, double x, double hi, int n=3 )
{
        ASSERT( n >= 1 );
        // mimic IEEE behavior
        if( isnan(x) || isnan(lo) || isnan(hi) )
                return false;
        if( fp_equal(lo,hi,n) )
                return fp_equal(0.5*(lo+hi),x,n);
        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
                return false;
        return true;
}
inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
{
        ASSERT( tol > 0.f );
        // mimic IEEE behavior
        if( isnan(x) || isnan(lo) || isnan(hi) )
                return false;
        if( fp_equal_tol(lo,hi,tol) )
                return fp_equal_tol(0.5f*(lo+hi),x,tol);
        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
                return false;
        return true;
}
inline bool fp_bound_tol( double lo, double x, double hi, double tol )
{
        ASSERT( tol > 0. );
        // mimic IEEE behavior
        if( isnan(x) || isnan(lo) || isnan(hi) )
                return false;
        if( fp_equal_tol(lo,hi,tol) )
                return