SYS/SYS_Math.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      Side Effects Software Inc
 *      477 Richmond Street West
 *      Toronto, Ontario
 *      Canada   M5V 3E7
 *      416-504-9876
 *
 * NAME:        SYS_Math.h ( SYS Library, C++)
 *
 * COMMENTS:    Houdini math interface. 
 *              statements like: a = sin(x); b = sinf(y);
 *              we use:          a = SYSsin(x); b = SYSsin(y);
 *
 * Functions defined here:
 *      SYSmin(a,b)                     - Returns the minimum of a or b
 *      SYSmin(a,b,c)                   - Returns the minimum of a, b or c
 *      SYSmin(a,b,c,d)                 - Returns minimum of a, b, c or d
 *      SYSargmin(a,b)                  - Returns index of minimum number
 *      SYSargmin(a,b,c)                - Returns index of minimum number
 *      SYSargmin(a,b,c,d)              - Returns index of minimum number
 *      SYSmax(a,b)                     - Returns the maximum of a or b
 *      SYSmax(a,b,c)                   - Returns the maximum of a, b or c
 *      SYSmax(a,b,c,d)                 - Returns maximum of a, b, c or d
 *      SYSsort(a, b)                   - Sorts so a < b
 *      SYSsort(a, b, c)                - Sorts so a < b < c
 *      SYSargmax(a,b)                  - Returns index of largest number
 *      SYSargmax(a,b,c)                - Returns index of largest number
 *      SYSargmax(a,b,c,d)              - Returns index of largest number
 *      SYSabs(x)                       - Works for all types (including fpreal)
 *      SYSavg(a,b,c)                   - Return average value of 3 values
 *      SYSavg(a,b,c,d)                 - Return average value of 4 values
 *
 *      SYSminmax(a,b,c,&min.&max)      - Find min and max of 3 values
 *      SYSminmax(a,b,c,d,&min.&max)    - Find min and max of 4 values
 *
 *  Comparison:
 *      SYSequalZero(a, [tol])          - Is a equal to zero
 *      SYSisEqual(a, b, tol)           - 
 *      SYSalmostEqual(a, b, ulps)      - Almost equal given units in last place
 *
 *  Interpolation/Range:
 *      SYSclamp(a, min, max)           - Clamp value between min and max
 *      SYSsmooth(min, max, a)          - Ease-in ease-out curve
 *      SYSsmooth(min, max, a, roll)    - Ease-in ease-out curve (with rolloff)
 *      SYSlerp(v1, v2, bias)           - Linear interpolation
 *      SYSfit(v, omin,omax, nmin,nmax) - Map v in (omin,omax) to (nmin,nmax)
 *
 *   The bilerp function expects:
 *        u0v1--u1v1
 *         |     |
 *        u0v0--u1v0
 *   Where u is horizontal and v is vertical in the diagram above.
 *      SYSbilerp(u0v0, u1v0, u0v1, u1v1, u, v)
 *
 *  Standard math (single/double precision signatures):
 *      SYSsin(x)
 *      SYScos(x)
 *      SYStan(x)
 *      SYSsqrt(x)
 *      SYSlog(x)
 *      SYSfabs(x)
 *      SYStrunc(x)
 *      SYSfloor(x)
 *      SYSceil(x)
 *      SYScbrt(x)
 *      SYSlog10(x)
 *      SYSfmod(x,y)
 *      SYSpow(x,y)
 *      SYSatan(x)
 *      SYSatan(y,x)    -- Note, SYSatan(y,x) is equivalent to SYSatan2(y,x)
 *      SYSatan2(y,x)
 *      SYShypot(x,y)
 *      SYSsafediv(x,y)         - Divide only if y != 0
 *      SYSsafesqrt(x)          - Compute sqrt only if x >= 0
 *
 *  Random:
 *      SYSpointerHash          - Generate a hash key for a pointer
 *      SYSwang_inthash         - 32-bit Integer hashing function
 *      SYSwang_inthash64       - 64-bit Integer hashing function
 *      SYSreal_hash            - Generate a hash for a real number
 *      SYSvector_hash          - Generate a hash for a real vector
 *      SYSfastRandom           - Really fast random number generator (0,1)
 *      SYSrandom               - Fast random number (fewer artifacts) (0,1)
 *      SYSfastRandomZero       - Really fast random number generator (-.5, .5)
 *      SYSrandomZero           - Fast random number (fewer artifacts) (-.5, .5)
 *
 *  Rounding:
 *      SYSroundDownToMultipleOf        -
 *      SYSroundUpToMultipleOf  -
 *      SYSroundAngle           -
 *
 *  Conversion:
 *      SYSatoi         - Convert to native integer format
 *      SYSatof         - Convert to native floating point format
 *      SYSatoi32       - Convert to signed 32 bit integer
 *      SYSatou32       - Convert to unsigned 32 bit integer
 *      SYSatoi64       - Convert to signed 64 bit integer
 *      SYSatou64       - Convert to unsigned 64 bit integer
 *      SYSatof32       - Convert to 32 bit floating point
 *      SYSatof64       - Convert to 64 bit floating point
 *
 *      In the strtoX conversions, the initial characters of the string are
 *      used to determine the base.  This uses the standard rules.  A leading
 *      0x or 0X indicates base 16.  A leading 0 indicates octal, and otherwise
 *      the string is assumed to be decimal.
 *      
 *      SYSstrtoi       - Convert to native integer format
 *      SYSstrtou       - Convert to native unsigned format
 *      SYSstrtoi32     - Convert to 32 bit integer
 *      SYSstrtou32     - Convert to 32 bit unsigned
 *      SYSstrtoi64     - Convert to 64 bit integer
 *      SYSstrtou64     - Convert to 64 bit unsigned
 *
 *  Misc:
 *      SYSisNan                - Is this a valid floating-point number?
 *      SYSisInt                - Is the string an integer
 *      SYSisFloat              - Is the string a float
 *      SYSisPrime              - Is it a prime?
 *      SYSsameSign(a,b)        - Do the non-zero numbers have the same sign?
 *                              (Note: does not work if either number is zero)
 *      SYSmakePrime            - Make it a prime
 *      SYSnextPrime            - Next prime
 *      SYSisPow2               - Is it a power of two?
 *      SYSmakePow2             - Make it a power of two
 *      SYSgetSinCosFromSlope   - Compute sin/cos given a slope
 *      SYShexCharToInt         - map '0'-'9,'a'-'z' to 0-15, or -1
 *      SYSintToHexChar         - map 0-15 to '0'-'9,'a'-'z'
 *
 *  For information on almostEqual, please see:
 *      http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
 *
 */

#ifndef __SYS_Math__
#define __SYS_Math__

#include "SYS_API.h"
#include <math.h>
#include <limits.h>
#include <float.h>
#include <stdlib.h>
#include "SYS_Types.h"

/*
 * System dependent includes
 */
#if defined(SOLARIS)
    #include <sunmath.h>
#endif

#if defined(WIN32) || defined(GAMEOS)
    #include "SYS_NTMath.h"
#endif

/* Some common constants */

#define FP32_TOLERANCE  ((fpreal32)0.00001)
#define FP64_TOLERANCE  ((fpreal64)0.00001)
#define FP32_REAL_MIN   ((fpreal32)2e-45)
#define FP64_REAL_MIN   ((fpreal64)2e-324)
#define FP32_REAL_MAX   ((fpreal32)FLT_MAX)
#define FP64_REAL_MAX   ((fpreal64)DBL_MAX)

// These are the precision you should use for writing out
// streams to guarantee that you fully save a double or a float.
// Note that limits.h defines DBL_DIG and FLT_DIG - they are
// defined incorrectly!
#define SYS_DBL_DIG     17
#define SYS_FLT_DIG     9

// Determine whether f is not-a-number (NAN) for any floating point type f
template< typename F >
inline bool SYSisNan(const F f)
{
    return(f != f);
}

SYS_API bool SYSisNan(const char *number);

SYS_API bool            SYSisInt(const char *str);
SYS_API bool            SYSisFloat(const char *str);
SYS_API fpreal32        SYSroundAngle(fpreal32 base, fpreal32 source);
SYS_API fpreal64        SYSroundAngle(fpreal64 base, fpreal64 source);
SYS_API fpreal64        SYSroundAngle(int64 base, int64 source);
SYS_API bool            SYSisPrime(uint num);
SYS_API bool            SYSisPrime(uint64 num);
SYS_API uint            SYSmakePrime(uint num);
SYS_API uint64          SYSmakePrime(uint64 num);
SYS_API uint            SYSmakePow2(uint num);
SYS_API uint64          SYSmakePow2(uint64 num);
#if defined(MBSD)
    SYS_API size_t      SYSmakePow2(size_t num);
#endif

static inline bool      SYSisPow2(int n)        { return n>0 && !(n&(n-1)); }
static inline bool      SYSisPow2(int64 n)      { return n>0 && !(n&(n-1)); }
static inline bool      SYSisPow2(uint n)       { return !(n&(n-1)); }
static inline bool      SYSisPow2(uint64 n)     { return !(n&(n-1)); }
#if defined(MBSD)
    static inline bool  SYSisPow2(size_t n)     { return !(n&(n-1)); }
#endif

static inline int32 SYSatoi32(const char *buf, int base=10)
{
    return base == 10 ? atoi(buf) : strtol(buf, 0, base);
}

static inline uint32 SYSatou32(const char *buf, int base=10)
{
    return strtoul(buf, 0, base);
}

static inline int64 SYSatoi64(const char *buf, int base=10)
{
#if defined(WIN32)
    return base == 10 ? _atoi64(buf) : _strtoi64(buf, 0, base);
#else
    return base == 10 ? atoll(buf) : strtoll(buf, 0, base);
#endif
}

static inline uint64 SYSatou64(const char *buf, int base=10)
{
#if defined(WIN32)
    return _strtoui64(buf, 0, base);
#else
    return strtoull(buf, 0, base);
#endif
}

static inline int32  SYSstrtoi32(const char *buf) { return strtol(buf, 0, 0); }
static inline uint32 SYSstrtou32(const char *buf) { return strtoul(buf, 0, 0); }
static inline int64  SYSstrtoi64(const char *buf)
{
#if defined(WIN32)
    return _strtoi64(buf, 0, 0);
#else
    return strtoll(buf, 0, 0);
#endif
}
static inline uint64 SYSstrtou64(const char *buf)
{
#if defined(WIN32)
    return _strtoui64(buf, 0, 0);
#else
    return strtoull(buf, 0, 0);
#endif
}

static inline fpreal32 SYSatof32(const char *buf)
{
    return atof(buf);
}

static inline fpreal64 SYSatof64(const char *buf)
{
    return atof(buf);
}

static inline  int SYSstrtoi(const char *buf) { return SYSstrtoi32(buf); }
static inline uint SYSstrtou(const char *buf) { return SYSstrtou32(buf); }

static inline int SYSatoi(const char *buf, int base=10)
            { return SYSatoi32(buf, base); }
static inline uint SYSatou(const char *buf, int base=10)
            { return SYSatou32(buf, base); }
static inline fpreal SYSatof(const char *buf)
            { return SYSatof32(buf); }

#if defined(__cplusplus)

// NOTE:
// These have been carefully written so that in the case of equality
// we always return the first parameter.  This is so that NANs in
// in the second parameter are suppressed.
#define h_min(a, b)     (((a) > (b)) ? (b) : (a))
#define h_argmin(a, b)  (((a) > (b)) ? 1 : 0)
#define h_max(a, b)     (((a) < (b)) ? (b) : (a))
#define h_argmax(a, b)  (((a) < (b)) ? 1 : 0)
// DO NOT CHANGE THE ABOVE WITHOUT READING THE COMMENT
#define h_abs(a)        (((a) > 0) ? (a) : -(a))

static inline  int16 SYSmin(int16 a, int16 b)           { return h_min(a,b); }
static inline  int16 SYSmax(int16 a, int16 b)           { return h_max(a,b); }
static inline  int16 SYSabs(int16 a)                    { return h_abs(a); }
static inline    int SYSargmin(int16 a, int16 b)        { return h_argmin(a,b);}
static inline    int SYSargmax(int16 a, int16 b)        { return h_argmax(a,b);}
static inline  int32 SYSmin(int32 a, int32 b)           { return h_min(a,b); }
static inline  int32 SYSmax(int32 a, int32 b)           { return h_max(a,b); }
static inline  int32 SYSabs(int32 a)                    { return h_abs(a); }
static inline    int SYSargmin(int32 a, int32 b)        { return h_argmin(a,b);}
static inline    int SYSargmax(int32 a, int32 b)        { return h_argmax(a,b);}
static inline  int64 SYSmin(int64 a, int64 b)           { return h_min(a,b); }
static inline  int64 SYSmax(int64 a, int64 b)           { return h_max(a,b); }
static inline  int64 SYSmin(int32 a, int64 b)           { return h_min(a,b); }
static inline  int64 SYSmax(int32 a, int64 b)           { return h_max(a,b); }
static inline  int64 SYSmin(int64 a, int32 b)           { return h_min(a,b); }
static inline  int64 SYSmax(int64 a, int32 b)           { return h_max(a,b); }
static inline  int64 SYSabs(int64 a)                    { return h_abs(a); }
static inline    int SYSargmin(int64 a, int64 b)        { return h_argmin(a,b);}
static inline    int SYSargmax(int64 a, int64 b)        { return h_argmax(a,b);}
static inline uint16 SYSmin(uint16 a, uint16 b)         { return h_min(a,b); }
static inline uint16 SYSmax(uint16 a, uint16 b)         { return h_max(a,b); }
static inline    int SYSargmin(uint16 a, uint16 b)      { return h_argmin(a,b);}
static inline    int SYSargmax(uint16 a, uint16 b)      { return h_argmax(a,b);}
static inline uint32 SYSmin(uint32 a, uint32 b)         { return h_min(a,b); }
static inline uint32 SYSmax(uint32 a, uint32 b)         { return h_max(a,b); }
static inline    int SYSargmin(uint32 a, uint32 b)      { return h_argmin(a,b);}
static inline    int SYSargmax(uint32 a, uint32 b)      { return h_argmax(a,b);}
static inline uint64 SYSmin(uint64 a, uint64 b)         { return h_min(a,b); }
static inline uint64 SYSmax(uint64 a, uint64 b)         { return h_max(a,b); }
static inline    int SYSargmin(uint64 a, uint64 b)      { return h_argmin(a,b);}
static inline    int SYSargmax(uint64 a, uint64 b)      { return h_argmax(a,b);}
static inline fpreal16 SYSmin(fpreal16 a, fpreal16 b)   { return h_min(a,b); }
static inline fpreal16 SYSmax(fpreal16 a, fpreal16 b)   { return h_max(a,b); }
static inline      int SYSargmin(fpreal16 a, fpreal16 b){ return h_argmin(a,b);}
static inline      int SYSargmax(fpreal16 a, fpreal16 b){ return h_argmax(a,b);}
static inline fpreal32 SYSmin(fpreal32 a, fpreal32 b)   { return h_min(a,b); }
static inline fpreal32 SYSmax(fpreal32 a, fpreal32 b)   { return h_max(a,b); }
static inline      int SYSargmin(fpreal32 a, fpreal32 b){ return h_argmin(a,b);}
static inline      int SYSargmax(fpreal32 a, fpreal32 b){ return h_argmax(a,b);}
static inline fpreal64 SYSmin(fpreal64 a, fpreal64 b)   { return h_min(a,b); }
static inline fpreal64 SYSmax(fpreal64 a, fpreal64 b)   { return h_max(a,b); }
static inline      int SYSargmin(fpreal64 a, fpreal64 b){ return h_argmin(a,b);}
static inline      int SYSargmax(fpreal64 a, fpreal64 b){ return h_argmax(a,b);}

// Some systems have size_t as a seperate type from uint.  Some don't.
#if (defined(LINUX) && defined(IA64)) || defined(MBSD)
static inline size_t SYSmin(size_t a, size_t b)         { return h_min(a,b); }
static inline size_t SYSmax(size_t a, size_t b)         { return h_max(a,b); }
static inline    int SYSargmin(size_t a, size_t b)      { return h_argmin(a,b);}
static inline    int SYSargmax(size_t a, size_t b)      { return h_argmax(a,b);}
#endif

#undef h_min
#undef h_max
#undef h_abs

#define h_clamp(val, min, max, tol)     \
            ((val <= min+tol) ? min : ((val >= max-tol) ? max : val))

    static inline int
    SYSclamp(int v, int min, int max)
        { return h_clamp(v, min, max, 0); }

    static inline uint
    SYSclamp(uint v, uint min, uint max)
        { return h_clamp(v, min, max, 0); }

    static inline int64
    SYSclamp(int64 v, int64 min, int64 max)
        { return h_clamp(v, min, max, CONST_INT64(0)); }

    static inline uint64
    SYSclamp(uint64 v, uint64 min, uint64 max)
        { return h_clamp(v, min, max, CONST_UINT64(0)); }

    static inline fpreal32
    SYSclamp(fpreal32 v, fpreal32 min, fpreal32 max, fpreal32 tol=(fpreal32)0)
        { return h_clamp(v, min, max, tol); }

    static inline fpreal64
    SYSclamp(fpreal64 v, fpreal64 min, fpreal64 max, fpreal64 tol=(fpreal64)0)
        { return h_clamp(v, min, max, tol); }

#undef h_clamp

#define SYS_UNARY(func)         \
    static inline fpreal64 SYS##func(fpreal64 arg) { return func(arg); } \
    static inline fpreal32 SYS##func(fpreal32 arg) { return func##f(arg); } \
    static inline fpreal64 SYS##func(int64 arg) { return func((fpreal64)arg); } \
    static inline fpreal64 SYS##func(int32 arg) { return func((fpreal64)arg); }

// No platform (to my knowledge) has cbrtf() defined, therefore:
#define cbrtf(x)        cbrt(x)
#define hypotf(x,y)     hypot((x),(y))

#define SYS_BINARY(func)                \
    static inline fpreal64 SYS##func(fpreal64 a, fpreal64 b) { return func(a,b); } \
    static inline fpreal32 SYS##func(fpreal32 a, fpreal32 b) \
                                    { return func##f(a, b); }

    SYS_UNARY(sin)
    SYS_UNARY(cos)
    SYS_UNARY(tan)
    SYS_UNARY(sinh)
    SYS_UNARY(cosh)
    SYS_UNARY(tanh)
    SYS_UNARY(sqrt)
    SYS_UNARY(trunc)
    SYS_UNARY(cbrt)
    SYS_UNARY(exp)
    SYS_BINARY(fmod)
    SYS_BINARY(pow)
    SYS_BINARY(atan2)
    SYS_BINARY(hypot)

    static inline fpreal32 SYSsafediv(fpreal32 x, fpreal32 y)
        { return y != 0 ? x/y : 0; }
    static inline fpreal64 SYSsafediv(fpreal64 x, fpreal64 y)
        { return y != 0 ? x/y : 0; }
    static inline fpreal32 SYSsafesqrt(fpreal32 x)
        { return x > 0 ? SYSsqrt(x) : 0; }
    static inline fpreal64 SYSsafesqrt(fpreal64 x)
        { return x > 0 ? SYSsqrt(x) : 0; }

    static inline fpreal32 SYSlog(fpreal32 v)   { return v <= 0 ? 0:logf(v); }
    static inline fpreal32 SYSlog(fpreal64 v)   { return v <= 0 ? 0:log(v); }
    static inline fpreal32 SYSlog10(fpreal32 v) { return v <= 0 ? 0:log10f(v); }
    static inline fpreal32 SYSlog10(fpreal64 v) { return v <= 0 ? 0:log10(v); }

#if defined(WIN32)
    static inline fpreal32 SYSexpm1(fpreal32 x) { return SYSexp(x) - 1; }
    static inline fpreal64 SYSexpm1(fpreal64 x) { return SYSexp(x) - 1; }
    static inline fpreal32 SYSlog1p(fpreal32 x) { return SYSlog(x+1); }
    static inline fpreal64 SYSlog1p(fpreal64 x) { return SYSexp(x+1); }
#else
    SYS_UNARY(expm1)
    SYS_UNARY(log1p)
#endif

#undef SYS_UNARY
#undef SYS_BINARY
#undef cbrtf
#undef hypotf

static inline fpreal32 SYSabs(fpreal32 a)       { return fabsf(a); }
static inline fpreal32 SYSabs(fpreal64 a)       { return  fabs(a); }
static inline fpreal32 SYSfabs(fpreal32 a)      { return fabsf(a); }
static inline fpreal32 SYSfabs(fpreal64 a)      { return  fabs(a); }

#include "SYS_Floor.h"

static inline fpreal32 SYSasin(fpreal32 a)
    { return asinf(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
static inline fpreal32 SYSacos(fpreal32 a)
    { return acosf(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
static inline fpreal32 SYSatan(fpreal32 a)
    { return atanf(a); }
static inline fpreal32 SYSatan(fpreal32 y, fpreal32 x)
    { return atan2f(y, x); }
static inline fpreal64 SYSasin(fpreal64 a)
    { return asin (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
static inline fpreal64 SYSacos(fpreal64 a)
    { return acos (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
static inline fpreal64 SYSatan(fpreal64 a)
    { return atan (a); }
static inline fpreal64 SYSatan(fpreal64 y, fpreal64 x)
    { return atan2(y, x); }

static inline fpreal32 SYSdegToRad(fpreal32 a) { return a*(fpreal32)(M_PI/180.0); }
static inline fpreal64 SYSdegToRad(fpreal64 a) { return a*(fpreal64)(M_PI/180.0); }
static inline fpreal64 SYSdegToRad(int64 a) { return a*(fpreal64)(M_PI/180.0); }

static inline fpreal32 SYSradToDeg(fpreal32 a) { return a*(fpreal32)(180.0/M_PI); }
static inline fpreal64 SYSradToDeg(fpreal64 a) { return a*(fpreal64)(180.0/M_PI); }
static inline fpreal64 SYSradToDeg(int64 a) { return a*(fpreal64)(180.0/M_PI); }

#define h_compare(func, code)   \
        static inline bool func(fpreal32 a, fpreal32 b, \
                                fpreal32 tol=FP32_TOLERANCE) \
            { return code; } \
        static inline bool func(fpreal64 a, fpreal64 b, \
                                fpreal64 tol=FP64_TOLERANCE) \
            { return code; }

    static inline bool
    SYSequalZero(fpreal32 a, fpreal32 tol=FP32_TOLERANCE)
        { return a >= -tol && a <= tol; }

    static inline bool
    SYSequalZero(fpreal64 a, fpreal64 tol=FP64_TOLERANCE)
        { return a >= -tol && a <= tol; }

    static inline bool
    SYSequalZero(int64 a, fpreal64 tol=FP64_TOLERANCE)
        { return a >= -tol && a <= tol; }


    h_compare(SYSisEqual, SYSabs(a-b)<=tol)
    h_compare(SYSisGreater, (a-b) > tol)
    h_compare(SYSisGreaterOrEqual, (a-b) >= -tol)
    h_compare(SYSisLess, (a-b) < -tol)
    h_compare(SYSisLessOrEqual, (a-b) <= tol)

inline bool SYSisEqual(int64 a, int64 b) { return a == b; }
inline bool SYSisGreater(int64 a, int64 b) { return a > b; }
inline bool SYSisGreaterOrEqual(int64 a, int64 b) { return a >= b; }
inline bool SYSisLess(int64 a, int64 b) { return a < b; }
inline bool SYSisLessOrEqual(int64 a, int64 b) { return a <= b; }

#undef h_compare

static inline bool
SYSalmostEqual(fpreal32 a, fpreal32 b, int ulps=50)
{
    SYS_FPReal32Union   ai, bi;
    ai.fval = a;
    bi.fval = b;
    if (ai.ival < 0)
        ai.ival = 0x80000000 - ai.ival;
    if (bi.ival < 0)
        bi.ival = 0x80000000 - bi.ival;
    return SYSabs(ai.ival-bi.ival) <= ulps;
}

static inline bool
SYSalmostEqual(fpreal64 a, fpreal64 b, int64 ulps)
{
    SYS_FPReal64Union   ai, bi;
    ai.fval = a;
    bi.fval = b;
    if (ai.ival < 0)
        ai.ival = CONST_INT64(0x8000000000000000) - ai.ival;
    if (bi.ival < 0)
        bi.ival = CONST_INT64(0x8000000000000000) - bi.ival;
    return SYSabs(ai.ival-bi.ival) <= ulps;
}

static inline bool
SYSalmostEqual(fpreal64 a, fpreal64 b, int32 ulps=50)
{
    return SYSalmostEqual(a, b, (int64)ulps);
}

#define h_max3(type)    \
            static inline type \
            SYSmax(type v0, type v1, type v2) { \
                if (v1 > v0) v0 = v1; \
                return SYSmax(v2, v0); \
            }
#define h_max4(type)    \
            static inline type \
            SYSmax(type v0, type v1, type v2, type v3) { \
                if (v1 > v0) v0 = v1; \
                if (v2 > v0) v0 = v2; \
                return SYSmax(v3, v0); \
            }
#define h_argmax3(type) \
            static inline int \
            SYSargmax(type v0, type v1, type v2) { \
                if (v1 > v0) return SYSargmax(v1, v2) + 1; \
                if (v2 > v0) return 2; \
                return 0; \
            }
#define h_argmax4(type) \
            static inline int \
            SYSargmax(type v0, type v1, type v2, type v3) { \
                if (v1 > v0) return SYSargmax(v1, v2, v3) + 1; \
                if (v2 > v0) return SYSargmax(v2, v3) + 2; \
                if (v3 > v0) return 3; \
                return 0; \
            }
#define h_min3(type)    \
            static inline type \
            SYSmin(type v0, type v1, type v2) { \
                if (v1 < v0) v0 = v1; \
                return SYSmin(v2, v0); \
            }
#define h_min4(type)    \
            static inline type \
            SYSmin(type v0, type v1, type v2, type v3) { \
                if (v1 < v0) v0 = v1; \
                if (v2 < v0) v0 = v2; \
                return SYSmin(v3, v0); \
            }
#define h_argmin3(type) \
            static inline int \
            SYSargmin(type v0, type v1, type v2) { \
                if (v1 < v0) return SYSargmin(v1, v2) + 1; \
                if (v2 < v0) return 2; \
                return 0; \
            }
#define h_argmin4(type) \
            static inline int \
            SYSargmin(type v0, type v1, type v2, type v3) { \
                if (v1 < v0) return SYSargmin(v1, v2, v3) + 1; \
                if (v2 < v0) return SYSargmin(v2, v3) + 2; \
                if (v3 < v0) return 3; \
                return 0; \
            }


#define h_max(type)     h_min3(type) h_min4(type) h_max3(type) h_max4(type) \
                        h_argmin3(type) h_argmin4(type)  \
                        h_argmax3(type) h_argmax4(type)

#define h_avg3(type)    \
            static inline type \
            SYSavg(type v0, type v1, type v2) { \
                return (v0+v1+v2) * ((type)(1.0/3.0)); \
            }
#define h_avg4(type)    \
            static inline type \
            SYSavg(type v0, type v1, type v2, type v3) { \
                return (v0+v1+v2+v3) * ((type)0.25); \
            }

#define h_avg(type)     h_avg3(type) h_avg4(type)

h_max(int8)
h_max(uint8)
h_max(int16)
h_max(uint16)
h_max(int32)
h_max(uint32)
h_max(int64)
h_max(uint64)
h_max(fpreal32) h_avg(fpreal32)
h_max(fpreal64) h_avg(fpreal64)

static inline int64
SYSavg(int64 a, int64 b, int64 c)
{ return (a + b + c + 1) / 3; }
static inline int64
SYSavg(int64 a, int64 b, int64 c, int64 d)
{ return (a + b + c + d + 2) / 4; }

// Some systems have size_t as a seperate type from uint.  Some don't.
#if (defined(LINUX) && defined(IA64)) || defined(MBSD)
h_max(size_t)
#endif

#undef h_min3
#undef h_min4
#undef h_max3
#undef h_max4
#undef h_argmin3
#undef h_argmin4
#undef h_argmax3
#undef h_argmax4
#undef h_avg3
#undef h_avg4
#undef h_max
#undef h_avg

/// @{
/// Linear interpolation between v0 and v1.
static inline fpreal32
SYSlerp(fpreal32 v1, fpreal32 v2, fpreal32 t)
{
    return v1 + (v2 - v1)*t;
}

static inline fpreal64
SYSlerp(fpreal64 v1, fpreal64 v2, fpreal64 t)
{
    return v1 + (v2 - v1)*t;
}
/// @}

/// @{
/// Bilinear interpolation over a quadrilateral: @code
///  (u=0,v=1) u0v1     u1v1 (u=1,v=1)
///              +--------+
///              |        |
///              |        |
///              +--------+
///  (u=0,v=0) u0v0     u1v0 (u=1,v=0)
/// @endcode
static inline fpreal32
SYSbilerp(fpreal32 u0v0, fpreal32 u1v0, fpreal32 u0v1, fpreal32 u1v1,
                        fpreal32 u, fpreal32 v)
{
    return SYSlerp(SYSlerp(u0v0, u0v1, v), SYSlerp(u1v0, u1v1, v), u);
}

static inline fpreal64
SYSbilerp(fpreal64 u0v0, fpreal64 u1v0, fpreal64 u0v1, fpreal64 u1v1,
                        fpreal64 u, fpreal64 v)
{
    return SYSlerp(SYSlerp(u0v0, u0v1, v), SYSlerp(u1v0, u1v1, v), u);
}
/// @}

/// @{
/// Barycentric coordinates do interpolation over the triangle specified by the
/// three vertices (v0, v1, v2): @code
/// |              (u=0,v=1)               |
/// |                  v2                  |
/// |                 / \                  |
/// |                /   \                 |
/// |               /     \                |
/// |    (u=0,v=0) v0------v1 (u=1,v=0)    |
/// @endcode
static inline fpreal32
SYSbarycentric(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 u, fpreal32 v)
{
    return v0*(1-u-v) + v1*u + v2*v;
}

static inline fpreal64
SYSbarycentric(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 u, fpreal64 v)
{
    return v0*(1-u-v) + v1*u + v2*v;
}
/// @}

static inline fpreal32
SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 val)
{
    fpreal32    t;
    if (val <= min) return 0;
    if (val >= max) return 1;
    t = max - min;
    if (SYSequalZero(t, (fpreal32)1e-8)) return (fpreal32).5;
    t = (val - min) / t;
    return t*t*((fpreal32)3.0 - (fpreal32)2.0*t);
}

static inline fpreal64
SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 val)
{
    fpreal64    t;
    if (val <= min) return 0;
    if (val >= max) return 1;
    t = max - min;
    if (SYSequalZero(t, (fpreal64)1e-18)) return (fpreal64).5;
    t = (val - min) / t;
    return t*t*((fpreal64)3.0 - (fpreal64)2.0*t);
}

static inline fpreal32
SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 value, fpreal32 roll)
{
    fpreal32    f;
    if (roll > 0)
    {
        f = SYSsmooth(min, max, value);
        return roll < (fpreal32)1 ? (fpreal32)1-SYSpow((fpreal32)1-f,
                                    (fpreal32)1/roll) : SYSpow(f, roll);
    }
    return (fpreal32)1;
}

static inline fpreal64
SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 value, fpreal64 roll)
{
    fpreal64    f;
    if (roll > 0)
    {
        f = SYSsmooth(min, max, value);
        return roll < (fpreal64)1 ? (fpreal64)1-SYSpow((fpreal64)1-f,
                                    (fpreal64)1/roll) : SYSpow(f, roll);
    }
    return (fpreal64)1;
}

static inline fpreal32
SYSfit(fpreal32 val, fpreal32 omin, fpreal32 omax, fpreal32 nmin, fpreal32 nmax)
{
    fpreal32    d = omax - omin;
    if (SYSequalZero(d, (fpreal32)1e-8))
        return (nmin+nmax)*(fpreal32).5;
    if (omin < omax)
    {
        if (val < omin) return nmin;
        if (val > omax) return nmax;
    }
    else
    {
        if (val < omax) return nmax;
        if (val > omin) return nmin;
    }
    return nmin + (nmax-nmin)*(val-omin)/d;
}

static inline fpreal64
SYSfit(fpreal64 val, fpreal64 omin, fpreal64 omax, fpreal64 nmin, fpreal64 nmax)
{
    fpreal64    d = omax - omin;
    if (SYSequalZero(d, (fpreal64)1e-18))
        return (nmin+nmax)*(fpreal64).5;
    if (omin < omax)
    {
        if (val < omin) return nmin;
        if (val > omax) return nmax;
    }
    else
    {
        if (val < omax) return nmax;
        if (val > omin) return nmin;
    }
    return nmin + (nmax-nmin)*(val-omin)/d;
}

static inline fpreal32
SYSfit01(fpreal32 val, fpreal32 nmin, fpreal32 nmax)
{
    if (val < 0) return nmin;
    if (val > 1) return nmax;
    return SYSlerp(nmin, nmax, val);
}

static inline fpreal64
SYSfit01(fpreal64 val, fpreal64 nmin, fpreal64 nmax)
{
    if (val < 0) return nmin;
    if (val > 1) return nmax;
    return SYSlerp(nmin, nmax, val);
}

static inline fpreal32
SYSroundDownToMultipleOf(fpreal32 val, fpreal32 multiple)
{
    fpreal32    modulus;
    fpreal32    retval;
    if (SYSequalZero(multiple))
        return val;
    modulus = SYSfmod(val, multiple);
    if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
        return val;
    retval = val - modulus;
    if( val < (fpreal32)0 && modulus!=(fpreal32)0 )
        retval -= multiple;
    return retval;
}

static inline fpreal64
SYSroundDownToMultipleOf(fpreal64 val, fpreal64 multiple)
{
    fpreal64    modulus;
    fpreal64    retval;
    if (SYSequalZero(multiple))
        return val;
    modulus = SYSfmod(val, multiple);
    if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
        return val;
    retval = val - modulus;
    if( val < (fpreal64)0 && modulus!=(fpreal64)0 )
        retval -= multiple;
    return retval;
}

static inline fpreal32
SYSroundUpToMultipleOf(fpreal32 val, fpreal32 multiple)
{
    fpreal32    modulus;
    if (SYSequalZero(multiple))
        return val;
    modulus = SYSfmod(val, multiple);
    if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
        return val;
    if (val > (fpreal32)0)
        val += multiple;
    return val - modulus;
}

static inline fpreal64
SYSroundUpToMultipleOf(fpreal64 val, fpreal64 multiple)
{
    fpreal64    modulus;
    if (SYSequalZero(multiple))
        return val;
    modulus = SYSfmod(val, multiple);
    if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
        return val;
    if (val > (fpreal64)0)
        val += multiple;
    return val - modulus;
}

static inline uint32
SYSwang_inthash(uint32 key)
{
    // From http://www.concentric.net/~Ttwang/tech/inthash.htm
    key += ~(key << 16);
    key ^=  (key >> 5);
    key +=  (key << 3);
    key ^=  (key >> 13);
    key += ~(key << 9);
    key ^=  (key >> 17);
    return key;
}

static inline uint32
SYSmultiplicative_inthash(uint32 key)
{
    // Multiply by the golden mean of 2^32 (Knuth's multiplicative method) to
    // get a uniformly distributed hash value.  Do a google search for
    // 2654435761 for more information.
    return key * 2654435761u;
}

static inline uint64
SYSwang_inthash64(uint64 key)
{
    // From http://www.concentric.net/~Ttwang/tech/inthash.htm
    key += ~(key << 32);
    key ^= (key >> 22);
    key += ~(key << 13);
    key ^= (key >> 8);
    key += (key << 3);
    key ^= (key >> 15);
    key += ~(key << 27);
    key ^= (key >> 31);
    return key;
}

static inline uint
SYSreal_hash(fpreal32 a, int lowmask=0xf)
{
    SYS_FPReal32Union   ai;
    ai.fval = a;
    return SYSwang_inthash(ai.uval & (~lowmask));
}

static inline uint
SYSreal_hash(fpreal64 a, int lowmask=0xf)
{
    SYS_FPReal64Union   ai;
    int32       lo, hi;
    ai.fval = a;
    lo = ai.uval & 0xfffffff;
    hi = ai.uval >> 32;
    return SYSwang_inthash(lo & (~lowmask)) ^ SYSwang_inthash(hi);
}

static inline uint
SYSvector_hash(const int32 *vector, int size)
{
    uint        hash = size;
    for (int i = 0; i < size; ++i)
        hash = hash*37 + SYSwang_inthash(vector[i]);
    return hash;
}

static inline uint
SYSvector_hash(const int64 *vector, int size)
{
    uint        hash = size;
    for (int i = 0; i < size; ++i)
        hash = hash*37 + SYSwang_inthash64(vector[i]);
    return hash;
}


static inline uint
SYSvector_hash(const fpreal32 *vector, int size)
{
    uint        hash = size;
    for (int i = 0; i < size; ++i)
        hash = hash*37 + SYSreal_hash(vector[i]);
    return hash;
}

static inline uint
SYSvector_hash(const fpreal64 *vector, int size)
{
    uint        hash = size;
    for (int i = 0; i < size; ++i)
        hash = hash*37 + SYSreal_hash(vector[i]);
    return hash;
}

static inline uint
SYSpointerHash(const void *ptr)
{
#if SIZEOF_VOID_P == 8
    // 64 bit pointers
    return (uint)((int64)ptr & 0xffffffff);
#else
    // 32 bit pointers
    return (uint)ptr;
#endif
}

// Generate a random number in range [0, 1)
inline static fpreal32
SYSfastRandom(uint &seed)
{
    SYS_FPReal32Union   tmp;
    seed = seed*1664525 + 1013904223;
    tmp.uval = 0x3f800000 | (0x007fffff & seed);
    return tmp.fval-1.0F;
}

inline static fpreal32
SYSrandom(uint &seed)
{
    SYS_FPReal32Union   tmp;
    seed = seed*1664525 + 1013904223;
    tmp.uval = SYSwang_inthash(seed);
    tmp.uval = 0x3f800000 | (0x007fffff & tmp.uval);
    return tmp.fval - 1.0F;
}

inline static uint
SYSfastRandomInt(uint &seed)
{
    seed = seed*1664525 + 1013904223;
    return SYSwang_inthash(seed);
}

inline static fpreal32
SYSfastRandomZero(uint &seed)
{
    SYS_FPReal32Union   tmp;
    seed = seed*1664525 + 1013904223;
    tmp.uval = 0x3f800000 | (0x007fffff & seed);
    return tmp.fval - 1.5F;
}

inline static fpreal32
SYSrandomZero(uint &seed)
{
    SYS_FPReal32Union   tmp;
    seed = seed*1664525 + 1013904223;
    tmp.uval = SYSwang_inthash(seed);
    tmp.uval = 0x3f800000 | (0x007fffff & tmp.uval);
    return tmp.fval - 1.5F;
}

// Find the minimum/maximum of four floating point values using the minimum
// number of comparisons.
inline static void
SYSminmax(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 v3,
         fpreal32 &min, fpreal32 &max)
{
    if (v0 < v1) { min = v0; max = v1; }
    else         { min = v1; max = v0; }
    if (v2 < v3)
    {
        if (v2 < min) min = v2;
        if (v3 > max) max = v3;
    }
    else
    {
        if (v2 > max) max = v2;
        if (v3 < min) min = v3;
    }
}

inline static void
SYSminmax(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 v3,
         fpreal64 &min, fpreal64 &max)
{
    if (v0 < v1) { min = v0; max = v1; }
    else         { min = v1; max = v0; }
    if (v2 < v3)
    {
        if (v2 < min) min = v2;
        if (v3 > max) max = v3;
    }
    else
    {
        if (v2 > max) max = v2;
        if (v3 < min) min = v3;
    }
}

inline static void
SYSminmax(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 &min, fpreal32 &max)
{
    if (v0 < v1) { min = v0; max = v1; }
    else         { min = v1; max = v0; }
    if (v2 < min)       { min = v2; }
    else if (v2 > max)  { max = v2; }
}

inline static void
SYSminmax(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 &min, fpreal64 &max)
{
    if (v0 < v1) { min = v0; max = v1; }
    else         { min = v1; max = v0; }
    if (v2 < min)       { min = v2; }
    else if (v2 > max)  { max = v2; }
}

inline static void
SYSgetSinCosFromSlope(fpreal32 slope, fpreal32 &sintheta, fpreal32 &costheta)
{
    fpreal32    one_over_m;
    sintheta = slope / SYSsqrt(slope*slope + (fpreal32)1);
    if ((slope = SYSabs(slope)) > (fpreal32)1)
    {
        one_over_m = (fpreal32)1 / slope;
        costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
    }
    else
        costheta = SYSsqrt((fpreal32)1 - sintheta*sintheta);
}

inline static void
SYSgetSinCosFromSlope(fpreal64 slope, fpreal64 &sintheta, fpreal64 &costheta)
{
    fpreal64    one_over_m;
    sintheta = slope / SYSsqrt(slope*slope + (fpreal64)1);
    if ((slope = SYSabs(slope)) > (fpreal64)1)
    {
        one_over_m = (fpreal64)1 / slope;
        costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
    }
    else
        costheta = SYSsqrt((fpreal64)1 - sintheta*sintheta);
}

inline static bool
SYSsameSign( fpreal32 v0, fpreal32 v1 )
{
    return (v0*v1)>0;
}

inline static bool
SYSsameSign( fpreal64 v0, fpreal64 v1 )
{
    return (v0*v1)>0;
}

inline static bool
SYSsameSign( int32 v0, int32 v1 )
{
    return (v0 ^ v1) >= 0;
}

inline static bool
SYSsameSign( int64 v0, int64 v1 )
{
    return (v0 ^ v1) >= 0;
}

static inline uint      SYSnextPrime(uint num)  { return SYSmakePrime(num+1); }

static inline int SYShexCharToInt(char c)
{
    // Given a hexidecimal character, return it's corresponding value between
    // 0 and 15, or -1 if it's not a hexidecimal character.
    if (c >= '0' && c <= '9')
        return c - '0';
    if (c >= 'a' && c <= 'f')
        return c - 'a' + 10;
    if (c >= 'A' && c <= 'F')
        return c - 'A' + 10;
    return -1;
}

static inline char SYSintToHexChar(int value)
{
    // The value must be from 0-15 for this function to return a valid result.
    return value < 10 ? '0' + value : 'a' + value - 10;
}

SYS_API void SYSsort(int &a, int &b);
SYS_API void SYSsort(int &a, int &b, int &c);
SYS_API void SYSsort(int64 &a, int64 &b);
SYS_API void SYSsort(int64 &a, int64 &b, int64 &c);
SYS_API void SYSsort(float &a, float &b);
SYS_API void SYSsort(float &a, float &b, float &c);
SYS_API void SYSsort(double &a, double &b);
SYS_API void SYSsort(double &a, double &b, double &c);

// Include permutations of fpreal32/fpreal64
#include "SYS_MathPermute.h"
#include "SYS_MathRestrictive.h"

#else   /* cplusplus */
    #define SYSmax(a,b) ((a) > (b) ? (a) : (b))
    #define SYSmin(a,b) ((a) < (b) ? (a) : (b))
    #define SYSabs(a)   ((a) < 0 ? (a) : -(a))
#endif

#endif

---