HDK: SYS/SYS_Math.h Source File

00001 /*
00002  * PROPRIETARY INFORMATION.  This software is proprietary to
00003  * Side Effects Software Inc., and is not to be reproduced,
00004  * transmitted, or disclosed in any way without written permission.
00005  *
00006  * Produced by:
00007  *      Side Effects Software Inc
00008  *      477 Richmond Street West
00009  *      Toronto, Ontario
00010  *      Canada   M5V 3E7
00011  *      416-504-9876
00012  *
00013  * NAME:        SYS_Math.h ( SYS Library, C++)
00014  *
00015  * COMMENTS:    Houdini math interface. 
00016  *              statements like: a = sin(x); b = sinf(y);
00017  *              we use:          a = SYSsin(x); b = SYSsin(y);
00018  *
00019  * Functions defined here:
00020  *      SYSmin(a,b)                     - Returns the minimum of a or b
00021  *      SYSmin(a,b,c)                   - Returns the minimum of a, b or c
00022  *      SYSmin(a,b,c,d)                 - Returns minimum of a, b, c or d
00023  *      SYSargmin(a,b)                  - Returns index of minimum number
00024  *      SYSargmin(a,b,c)                - Returns index of minimum number
00025  *      SYSargmin(a,b,c,d)              - Returns index of minimum number
00026  *      SYSmax(a,b)                     - Returns the maximum of a or b
00027  *      SYSmax(a,b,c)                   - Returns the maximum of a, b or c
00028  *      SYSmax(a,b,c,d)                 - Returns maximum of a, b, c or d
00029  *      SYSargmax(a,b)                  - Returns index of largest number
00030  *      SYSargmax(a,b,c)                - Returns index of largest number
00031  *      SYSargmax(a,b,c,d)              - Returns index of largest number
00032  *      SYSabs(x)                       - Works for all types (including fpreal)
00033  *      SYSavg(a,b,c)                   - Return average value of 3 values
00034  *      SYSavg(a,b,c,d)                 - Return average value of 4 values
00035  *
00036  *      SYSminmax(a,b,c,&min.&max)      - Find min and max of 3 values
00037  *      SYSminmax(a,b,c,d,&min.&max)    - Find min and max of 4 values
00038  *
00039  *  Comparison:
00040  *      SYSequalZero(a, [tol])          - Is a equal to zero
00041  *      SYSisEqual(a, b, tol)           - 
00042  *      SYSalmostEqual(a, b, ulps)      - Almost equal given units in last place
00043  *
00044  *  Interpolation/Range:
00045  *      SYSclamp(a, min, max)           - Clamp value between min and max
00046  *      SYSsmooth(min, max, a)          - Ease-in ease-out curve
00047  *      SYSsmooth(min, max, a, roll)    - Ease-in ease-out curve (with rolloff)
00048  *      SYSlerp(v1, v2, bias)           - Linear interpolation
00049  *      SYSfit(v, omin,omax, nmin,nmax) - Map v in (omin,omax) to (nmin,nmax)
00050  *
00051  *   The bilerp function expects:
00052  *        v01--v11
00053  *         |    |
00054  *        v00--v10
00055  *   Where u is horizontal and v is vertical in the diagram above.
00056  *      SYSbilerp(v00, v10, v01, v11, u, v)
00057  *
00058  *  Standard math (single/double precision signatures):
00059  *      SYSsin(x)
00060  *      SYScos(x)
00061  *      SYStan(x)
00062  *      SYSsqrt(x)
00063  *      SYSlog(x)
00064  *      SYSfabs(x)
00065  *      SYStrunc(x)
00066  *      SYSfloor(x)
00067  *      SYSceil(x)
00068  *      SYScbrt(x)
00069  *      SYSlog10(x)
00070  *      SYSfmod(x,y)
00071  *      SYSpow(x,y)
00072  *      SYSatan(x)
00073  *      SYSatan(y,x)    -- Note, SYSatan(y,x) is equivalent to SYSatan2(y,x)
00074  *      SYSatan2(y,x)
00075  *      SYShypot(x,y)
00076  *      SYSsafediv(x,y)         - Divide only if y != 0
00077  *      SYSsafesqrt(x)          - Compute sqrt only if x >= 0
00078  *
00079  *  Random:
00080  *      SYSpointerHash          - Generate a hash key for a pointer
00081  *      SYSwang_inthash         - 32-bit Integer hashing function
00082  *      SYSwang_inthash64       - 64-bit Integer hashing function
00083  *      SYSreal_hash            - Generate a hash for a real number
00084  *      SYSfastRandom           - Really fast random number generator (0,1)
00085  *      SYSrandom               - Fast random number (fewer artifacts) (0,1)
00086  *      SYSfastRandomZero       - Really fast random number generator (-.5, .5)
00087  *      SYSrandomZero           - Fast random number (fewer artifacts) (-.5, .5)
00088  *
00089  *  Rounding:
00090  *      SYSroundDownToMultipleOf        -
00091  *      SYSroundUpToMultipleOf  -
00092  *      SYSroundAngle           -
00093  *
00094  *  Conversion:
00095  *      SYSatoi         - Convert to native integer format
00096  *      SYSatof         - Convert to native floating point format
00097  *      SYSatoi32       - Convert to signed 32 bit integer
00098  *      SYSatou32       - Convert to unsigned 32 bit integer
00099  *      SYSatoi64       - Convert to signed 64 bit integer
00100  *      SYSatou64       - Convert to unsigned 64 bit integer
00101  *      SYSatof32       - Convert to 32 bit floating point
00102  *      SYSatof64       - Convert to 64 bit floating point
00103  *
00104  *      In the strtoX conversions, the initial characters of the string are
00105  *      used to determine the base.  This uses the standard rules.  A leading
00106  *      0x or 0X indicates base 16.  A leading 0 indicates octal, and otherwise
00107  *      the string is assumed to be decimal.
00108  *      
00109  *      SYSstrtoi       - Convert to native integer format
00110  *      SYSstrtou       - Convert to native unsigned format
00111  *      SYSstrtoi32     - Convert to 32 bit integer
00112  *      SYSstrtou32     - Convert to 32 bit unsigned
00113  *      SYSstrtoi64     - Convert to 64 bit integer
00114  *      SYSstrtou64     - Convert to 64 bit unsigned
00115  *
00116  *  Misc:
00117  *      SYSisNan                - Is this a valid floating-point number?
00118  *      SYSisInt                - Is the string an integer
00119  *      SYSisFloat              - Is the string a float
00120  *      SYSisPrime              - Is it a prime?
00121  *      SYSsameSign(a,b)        - Do the non-zero numbers have the same sign?
00122  *                              (Note: does not work if either number is zero)
00123  *      SYSmakePrime            - Make it a prime
00124  *      SYSnextPrime            - Next prime
00125  *      SYSisPow2               - Is it a power of two?
00126  *      SYSmakePow2             - Make it a power of two
00127  *      SYSgetSinCosFromSlope   - Compute sin/cos given a slope
00128  *
00129  *  For information on almostEqual, please see:
00130  *      http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
00131  *
00132  */
00133 
00134 #ifndef __SYS_Math__
00135 #define __SYS_Math__
00136 
00137 #include "SYS_API.h"
00138 #include <math.h>
00139 #include <limits.h>
00140 #include <float.h>
00141 #include <stdlib.h>
00142 #include "SYS_Types.h"
00143 
00144 /*
00145  * System dependent includes
00146  */
00147 #if defined(SOLARIS)
00148     #include <sunmath.h>
00149 #endif
00150 
00151 #if defined(WIN32) || defined(GAMEOS)
00152     #include "SYS_NTMath.h"
00153 #endif
00154 
00155 /* Some common constants */
00156 
00157 #define FP32_TOLERANCE  ((fpreal32)0.00001)
00158 #define FP64_TOLERANCE  ((fpreal64)0.00001)
00159 #define FP32_REAL_MIN   ((fpreal32)2e-45)
00160 #define FP64_REAL_MIN   ((fpreal64)2e-324)
00161 #define FP32_REAL_MAX   ((fpreal32)FLT_MAX)
00162 #define FP64_REAL_MAX   ((fpreal64)DBL_MAX)
00163 
00164 // These are the precision you should use for writing out
00165 // streams to guarantee that you fully save a double or a float.
00166 // Note that limits.h defines DBL_DIG and FLT_DIG - they are
00167 // defined incorrectly!
00168 #define SYS_DBL_DIG     17
00169 #define SYS_FLT_DIG     9
00170 
00171 // Determine whether f is not-a-number (NAN) for any floating point type f
00172 template< typename F >
00173 inline bool SYSisNan(const F f)
00174 {
00175     return(f != f);
00176 }
00177 
00178 SYS_API bool            SYSisInt(const char *str);
00179 SYS_API bool            SYSisFloat(const char *str);
00180 SYS_API fpreal32        SYSroundAngle(fpreal32 base, fpreal32 source);
00181 SYS_API fpreal64        SYSroundAngle(fpreal64 base, fpreal64 source);
00182 SYS_API bool            SYSisPrime(uint num);
00183 SYS_API bool            SYSisPrime(uint64 num);
00184 SYS_API uint            SYSmakePrime(uint num);
00185 SYS_API uint64          SYSmakePrime(uint64 num);
00186 SYS_API bool            SYSisPow2(uint num);
00187 SYS_API bool            SYSisPow2(uint64 num);
00188 SYS_API uint            SYSmakePow2(uint num);
00189 SYS_API uint64          SYSmakePow2(uint64 num);
00190 #if defined(MBSD)
00191     SYS_API bool        SYSisPow2(size_t num);
00192     SYS_API size_t      SYSmakePow2(size_t num);
00193 #endif
00194 
00195 static inline int32 SYSatoi32(const char *buf, int base=10)
00196 {
00197     return base == 10 ? atoi(buf) : strtol(buf, 0, base);
00198 }
00199 
00200 static inline uint32 SYSatou32(const char *buf, int base=10)
00201 {
00202     return strtoul(buf, 0, base);
00203 }
00204 
00205 static inline int64 SYSatoi64(const char *buf, int base=10)
00206 {
00207 #if defined(WIN32)
00208     return base == 10 ? _atoi64(buf) : _strtoi64(buf, 0, base);
00209 #else
00210     return base == 10 ? atoll(buf) : strtoll(buf, 0, base);
00211 #endif
00212 }
00213 
00214 static inline uint64 SYSatou64(const char *buf, int base=10)
00215 {
00216 #if defined(WIN32)
00217     return _strtoui64(buf, 0, base);
00218 #else
00219     return strtoull(buf, 0, base);
00220 #endif
00221 }
00222 
00223 static inline int32  SYSstrtoi32(const char *buf) { return strtol(buf, 0, 0); }
00224 static inline uint32 SYSstrtou32(const char *buf) { return strtoul(buf, 0, 0); }
00225 static inline int64  SYSstrtoi64(const char *buf)
00226 {
00227 #if defined(WIN32)
00228     return _strtoi64(buf, 0, 0);
00229 #else
00230     return strtoll(buf, 0, 0);
00231 #endif
00232 }
00233 static inline uint64 SYSstrtou64(const char *buf)
00234 {
00235 #if defined(WIN32)
00236     return _strtoui64(buf, 0, 0);
00237 #else
00238     return strtoull(buf, 0, 0);
00239 #endif
00240 }
00241 
00242 static inline fpreal32 SYSatof32(const char *buf)
00243 {
00244     return atof(buf);
00245 }
00246 
00247 static inline fpreal64 SYSatof64(const char *buf)
00248 {
00249     return atof(buf);
00250 }
00251 
00252 static inline  int SYSstrtoi(const char *buf) { return SYSstrtoi32(buf); }
00253 static inline uint SYSstrtou(const char *buf) { return SYSstrtou32(buf); }
00254 
00255 static inline int SYSatoi(const char *buf, int base=10)
00256             { return SYSatoi32(buf, base); }
00257 static inline uint SYSatou(const char *buf, int base=10)
00258             { return SYSatou32(buf, base); }
00259 static inline fpreal SYSatof(const char *buf)
00260             { return SYSatof32(buf); }
00261 
00262 #if defined(__cplusplus)
00263 #define h_min(a, b)     (((a) > (b)) ? (b) : (a))
00264 #define h_argmin(a, b)  (((a) > (b)) ? 1 : 0)
00265 #define h_max(a, b)     (((a) < (b)) ? (b) : (a))
00266 #define h_argmax(a, b)  (((a) < (b)) ? 1 : 0)
00267 #define h_abs(a)        (((a) > 0) ? (a) : -(a))
00268 
00269 static inline  int16 SYSmin(int16 a, int16 b)           { return h_min(a,b); }
00270 static inline  int16 SYSmax(int16 a, int16 b)           { return h_max(a,b); }
00271 static inline  int16 SYSabs(int16 a)                    { return h_abs(a); }
00272 static inline    int SYSargmin(int16 a, int16 b)        { return h_argmin(a,b);}
00273 static inline    int SYSargmax(int16 a, int16 b)        { return h_argmax(a,b);}
00274 static inline  int32 SYSmin(int32 a, int32 b)           { return h_min(a,b); }
00275 static inline  int32 SYSmax(int32 a, int32 b)           { return h_max(a,b); }
00276 static inline  int32 SYSabs(int32 a)                    { return h_abs(a); }
00277 static inline    int SYSargmin(int32 a, int32 b)        { return h_argmin(a,b);}
00278 static inline    int SYSargmax(int32 a, int32 b)        { return h_argmax(a,b);}
00279 static inline  int64 SYSmin(int64 a, int64 b)           { return h_min(a,b); }
00280 static inline  int64 SYSmax(int64 a, int64 b)           { return h_max(a,b); }
00281 static inline  int64 SYSabs(int64 a)                    { return h_abs(a); }
00282 static inline    int SYSargmin(int64 a, int64 b)        { return h_argmin(a,b);}
00283 static inline    int SYSargmax(int64 a, int64 b)        { return h_argmax(a,b);}
00284 static inline uint16 SYSmin(uint16 a, uint16 b)         { return h_min(a,b); }
00285 static inline uint16 SYSmax(uint16 a, uint16 b)         { return h_max(a,b); }
00286 static inline    int SYSargmin(uint16 a, uint16 b)      { return h_argmin(a,b);}
00287 static inline    int SYSargmax(uint16 a, uint16 b)      { return h_argmax(a,b);}
00288 static inline uint32 SYSmin(uint32 a, uint32 b)         { return h_min(a,b); }
00289 static inline uint32 SYSmax(uint32 a, uint32 b)         { return h_max(a,b); }
00290 static inline    int SYSargmin(uint32 a, uint32 b)      { return h_argmin(a,b);}
00291 static inline    int SYSargmax(uint32 a, uint32 b)      { return h_argmax(a,b);}
00292 static inline uint64 SYSmin(uint64 a, uint64 b)         { return h_min(a,b); }
00293 static inline uint64 SYSmax(uint64 a, uint64 b)         { return h_max(a,b); }
00294 static inline    int SYSargmin(uint64 a, uint64 b)      { return h_argmin(a,b);}
00295 static inline    int SYSargmax(uint64 a, uint64 b)      { return h_argmax(a,b);}
00296 static inline fpreal16 SYSmin(fpreal16 a, fpreal16 b)   { return h_min(a,b); }
00297 static inline fpreal16 SYSmax(fpreal16 a, fpreal16 b)   { return h_max(a,b); }
00298 static inline      int SYSargmin(fpreal16 a, fpreal16 b){ return h_argmin(a,b);}
00299 static inline      int SYSargmax(fpreal16 a, fpreal16 b){ return h_argmax(a,b);}
00300 static inline fpreal32 SYSmin(fpreal32 a, fpreal32 b)   { return h_min(a,b); }
00301 static inline fpreal32 SYSmax(fpreal32 a, fpreal32 b)   { return h_max(a,b); }
00302 static inline      int SYSargmin(fpreal32 a, fpreal32 b){ return h_argmin(a,b);}
00303 static inline      int SYSargmax(fpreal32 a, fpreal32 b){ return h_argmax(a,b);}
00304 static inline fpreal64 SYSmin(fpreal64 a, fpreal64 b)   { return h_min(a,b); }
00305 static inline fpreal64 SYSmax(fpreal64 a, fpreal64 b)   { return h_max(a,b); }
00306 static inline      int SYSargmin(fpreal64 a, fpreal64 b){ return h_argmin(a,b);}
00307 static inline      int SYSargmax(fpreal64 a, fpreal64 b){ return h_argmax(a,b);}
00308 
00309 // Some systems have size_t as a seperate type from uint.  Some don't.
00310 #if (defined(LINUX) && defined(IA64)) || defined(MBSD)
00311 static inline size_t SYSmin(size_t a, size_t b)         { return h_min(a,b); }
00312 static inline size_t SYSmax(size_t a, size_t b)         { return h_max(a,b); }
00313 static inline    int SYSargmin(size_t a, size_t b)      { return h_argmin(a,b);}
00314 static inline    int SYSargmax(size_t a, size_t b)      { return h_argmax(a,b);}
00315 #endif
00316 
00317 #undef h_min
00318 #undef h_max
00319 #undef h_abs
00320 
00321 #define h_clamp(val, min, max, tol)     \
00322             ((val <= min+tol) ? min : ((val >= max-tol) ? max : val))
00323 
00324     static inline int
00325     SYSclamp(int v, int min, int max)
00326         { return h_clamp(v, min, max, 0); }
00327 
00328     static inline uint
00329     SYSclamp(uint v, uint min, uint max)
00330         { return h_clamp(v, min, max, 0); }
00331 
00332     static inline int64
00333     SYSclamp(int64 v, int64 min, int64 max)
00334         { return h_clamp(v, min, max, CONST_INT64(0)); }
00335 
00336     static inline uint64
00337     SYSclamp(uint64 v, uint64 min, uint64 max)
00338         { return h_clamp(v, min, max, CONST_UINT64(0)); }
00339 
00340     static inline fpreal32
00341     SYSclamp(fpreal32 v, fpreal32 min, fpreal32 max, fpreal32 tol=(fpreal32)0)
00342         { return h_clamp(v, min, max, tol); }
00343 
00344     static inline fpreal64
00345     SYSclamp(fpreal64 v, fpreal64 min, fpreal64 max, fpreal64 tol=(fpreal64)0)
00346         { return h_clamp(v, min, max, tol); }
00347 
00348 #undef h_clamp
00349 
00350 #define SYS_UNARY(func)         \
00351     static inline fpreal64 SYS##func(fpreal64 arg) { return func(arg); } \
00352     static inline fpreal32 SYS##func(fpreal32 arg) { return func##f(arg); }
00353 
00354 // No platform (to my knowledge) has cbrtf() defined, therefore:
00355 #define cbrtf(x)        cbrt(x)
00356 #define hypotf(x,y)     hypot((x),(y))
00357 
00358 #define SYS_BINARY(func)                \
00359     static inline fpreal64 SYS##func(fpreal64 a, fpreal64 b) { return func(a,b); } \
00360     static inline fpreal32 SYS##func(fpreal32 a, fpreal32 b) \
00361                                     { return func##f(a, b); }
00362 
00363     SYS_UNARY(sin)
00364     SYS_UNARY(cos)
00365     SYS_UNARY(tan)
00366     SYS_UNARY(sinh)
00367     SYS_UNARY(cosh)
00368     SYS_UNARY(tanh)
00369     SYS_UNARY(sqrt)
00370     SYS_UNARY(trunc)
00371     SYS_UNARY(cbrt)
00372     SYS_UNARY(exp)
00373     SYS_BINARY(fmod)
00374     SYS_BINARY(pow)
00375     SYS_BINARY(atan2)
00376     SYS_BINARY(hypot)
00377 
00378     static inline fpreal32 SYSsafediv(fpreal32 x, fpreal32 y)
00379         { return y != 0 ? x/y : 0; }
00380     static inline fpreal64 SYSsafediv(fpreal64 x, fpreal64 y)
00381         { return y != 0 ? x/y : 0; }
00382     static inline fpreal32 SYSsafesqrt(fpreal32 x)
00383         { return x > 0 ? SYSsqrt(x) : 0; }
00384     static inline fpreal64 SYSsafesqrt(fpreal64 x)
00385         { return x > 0 ? SYSsqrt(x) : 0; }
00386 
00387     static inline fpreal32 SYSlog(fpreal32 v)   { return v <= 0 ? 0:logf(v); }
00388     static inline fpreal32 SYSlog(fpreal64 v)   { return v <= 0 ? 0:log(v); }
00389     static inline fpreal32 SYSlog10(fpreal32 v) { return v <= 0 ? 0:log10f(v); }
00390     static inline fpreal32 SYSlog10(fpreal64 v) { return v <= 0 ? 0:log10(v); }
00391 
00392 #if defined(WIN32)
00393     static inline fpreal32 SYSexpm1(fpreal32 x) { return SYSexp(x) - 1; }
00394     static inline fpreal64 SYSexpm1(fpreal64 x) { return SYSexp(x) - 1; }
00395     static inline fpreal32 SYSlog1p(fpreal32 x) { return SYSlog(x+1); }
00396     static inline fpreal64 SYSlog1p(fpreal64 x) { return SYSexp(x+1); }
00397 #else
00398     SYS_UNARY(expm1)
00399     SYS_UNARY(log1p)
00400 #endif
00401 
00402 #undef SYS_UNARY
00403 #undef SYS_BINARY
00404 #undef cbrtf
00405 #undef hypotf
00406 
00407 static inline fpreal32 SYSabs(fpreal32 a)       { return fabsf(a); }
00408 static inline fpreal32 SYSabs(fpreal64 a)       { return  fabs(a); }
00409 static inline fpreal32 SYSfabs(fpreal32 a)      { return fabsf(a); }
00410 static inline fpreal32 SYSfabs(fpreal64 a)      { return  fabs(a); }
00411 
00412 #include "SYS_Floor.h"
00413 
00414 static inline fpreal32 SYSasin(fpreal32 a)
00415     { return asinf(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
00416 static inline fpreal32 SYSacos(fpreal32 a)
00417     { return acosf(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
00418 static inline fpreal32 SYSatan(fpreal32 a)
00419     { return atanf(a); }
00420 static inline fpreal32 SYSatan(fpreal32 y, fpreal32 x)
00421     { return atan2f(y, x); }
00422 static inline fpreal64 SYSasin(fpreal64 a)
00423     { return asin (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
00424 static inline fpreal64 SYSacos(fpreal64 a)
00425     { return acos (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
00426 static inline fpreal64 SYSatan(fpreal64 a)
00427     { return atan (a); }
00428 static inline fpreal64 SYSatan(fpreal64 y, fpreal64 x)
00429     { return atan2(y, x); }
00430 
00431 static inline fpreal32 SYSdegToRad(fpreal32 a) { return a*(fpreal32)(M_PI/180.0); }
00432 static inline fpreal64 SYSdegToRad(fpreal64 a) { return a*(fpreal64)(M_PI/180.0); }
00433 
00434 static inline fpreal32 SYSradToDeg(fpreal32 a) { return a*(fpreal32)(180.0/M_PI); }
00435 static inline fpreal64 SYSradToDeg(fpreal64 a) { return a*(fpreal64)(180.0/M_PI); }
00436 
00437 #define h_compare(func, code)   \
00438         static inline bool func(fpreal32 a, fpreal32 b, \
00439                                 fpreal32 tol=FP32_TOLERANCE) \
00440             { return code; } \
00441         static inline bool func(fpreal64 a, fpreal64 b, \
00442                                 fpreal64 tol=FP64_TOLERANCE) \
00443             { return code; }
00444 
00445     static inline bool
00446     SYSequalZero(fpreal32 a, fpreal32 tol=FP32_TOLERANCE)
00447         { return a >= -tol && a <= tol; }
00448 
00449     static inline bool
00450     SYSequalZero(fpreal64 a, fpreal64 tol=FP64_TOLERANCE)
00451         { return a >= -tol && a <= tol; }
00452 
00453     h_compare(SYSisEqual, SYSabs(a-b)<=tol)
00454     h_compare(SYSisGreater, (a-b) > tol)
00455     h_compare(SYSisGreaterOrEqual, (a-b) >= -tol)
00456     h_compare(SYSisLess, (a-b) < -tol)
00457     h_compare(SYSisLessOrEqual, (a-b) <= tol)
00458 
00459 #undef h_compare
00460 
00461 static inline bool
00462 SYSalmostEqual(fpreal32 a, fpreal32 b, int ulps=50)
00463 {
00464     SYS_FPReal32Union   ai, bi;
00465     ai.fval = a;
00466     bi.fval = b;
00467     if (ai.ival < 0)
00468         ai.ival = 0x80000000 - ai.ival;
00469     if (bi.ival < 0)
00470         bi.ival = 0x80000000 - bi.ival;
00471     return SYSabs(ai.ival-bi.ival) <= ulps;
00472 }
00473 
00474 static inline bool
00475 SYSalmostEqual(fpreal64 a, fpreal64 b, int64 ulps=50)
00476 {
00477     SYS_FPReal64Union   ai, bi;
00478     ai.fval = a;
00479     bi.fval = b;
00480     if (ai.ival < 0)
00481         ai.ival = CONST_INT64(0x8000000000000000) - ai.ival;
00482     if (bi.ival < 0)
00483         bi.ival = CONST_INT64(0x8000000000000000) - bi.ival;
00484     return SYSabs(ai.ival-bi.ival) <= ulps;
00485 }
00486 
00487 static inline bool
00488 SYSalmostEqual(fpreal64 a, fpreal64 b, int32 ulps=50)
00489 {
00490     return SYSalmostEqual(a, b, (int64)ulps);
00491 }
00492 
00493 #define h_max3(type)    \
00494             static inline type \
00495             SYSmax(type v0, type v1, type v2) { \
00496                 if (v1 > v0) v0 = v1; \
00497                 return SYSmax(v2, v0); \
00498             }
00499 #define h_max4(type)    \
00500             static inline type \
00501             SYSmax(type v0, type v1, type v2, type v3) { \
00502                 if (v1 > v0) v0 = v1; \
00503                 if (v2 > v0) v0 = v2; \
00504                 return SYSmax(v3, v0); \
00505             }
00506 #define h_argmax3(type) \
00507             static inline int \
00508             SYSargmax(type v0, type v1, type v2) { \
00509                 if (v1 > v0) return SYSargmax(v1, v2) + 1; \
00510                 if (v2 > v0) return 2; \
00511                 return 0; \
00512             }
00513 #define h_argmax4(type) \
00514             static inline int \
00515             SYSargmax(type v0, type v1, type v2, type v3) { \
00516                 if (v1 > v0) return SYSargmax(v1, v2, v3) + 1; \
00517                 if (v2 > v0) return SYSargmax(v2, v3) + 2; \
00518                 if (v3 > v0) return 3; \
00519                 return 0; \
00520             }
00521 #define h_min3(type)    \
00522             static inline type \
00523             SYSmin(type v0, type v1, type v2) { \
00524                 if (v1 < v0) v0 = v1; \
00525                 return SYSmin(v2, v0); \
00526             }
00527 #define h_min4(type)    \
00528             static inline type \
00529             SYSmin(type v0, type v1, type v2, type v3) { \
00530                 if (v1 < v0) v0 = v1; \
00531                 if (v2 < v0) v0 = v2; \
00532                 return SYSmin(v3, v0); \
00533             }
00534 #define h_argmin3(type) \
00535             static inline int \
00536             SYSargmin(type v0, type v1, type v2) { \
00537                 if (v1 < v0) return SYSargmin(v1, v2) + 1; \
00538                 if (v2 < v0) return 2; \
00539                 return 0; \
00540             }
00541 #define h_argmin4(type) \
00542             static inline int \
00543             SYSargmin(type v0, type v1, type v2, type v3) { \
00544                 if (v1 < v0) return SYSargmin(v1, v2, v3) + 1; \
00545                 if (v2 < v0) return SYSargmin(v2, v3) + 2; \
00546                 if (v3 < v0) return 3; \
00547                 return 0; \
00548             }
00549 
00550 
00551 #define h_max(type)     h_min3(type) h_min4(type) h_max3(type) h_max4(type) \
00552                         h_argmin3(type) h_argmin4(type)  \
00553                         h_argmax3(type) h_argmax4(type)
00554 
00555 #define h_avg3(type)    \
00556             static inline type \
00557             SYSavg(type v0, type v1, type v2) { \
00558                 return (v0+v1+v2) * ((type)(1.0/3.0)); \
00559             }
00560 #define h_avg4(type)    \
00561             static inline type \
00562             SYSavg(type v0, type v1, type v2, type v3) { \
00563                 return (v0+v1+v2+v3) * ((type)0.25); \
00564             }
00565 
00566 #define h_avg(type)     h_avg3(type) h_avg4(type)
00567 
00568 h_max(int8)
00569 h_max(uint8)
00570 h_max(int16)
00571 h_max(uint16)
00572 h_max(int32)
00573 h_max(uint32)
00574 h_max(int64)
00575 h_max(uint64)
00576 h_max(fpreal32) h_avg(fpreal32)
00577 h_max(fpreal64) h_avg(fpreal64)
00578 
00579 // Some systems have size_t as a seperate type from uint.  Some don't.
00580 #if (defined(LINUX) && defined(IA64)) || defined(MBSD)
00581 h_max(size_t)
00582 #endif
00583 
00584 #undef h_min3
00585 #undef h_min4
00586 #undef h_max3
00587 #undef h_max4
00588 #undef h_argmin3
00589 #undef h_argmin4
00590 #undef h_argmax3
00591 #undef h_argmax4
00592 #undef h_avg3
00593 #undef h_avg4
00594 #undef h_max
00595 #undef h_avg
00596 
00597 static inline fpreal32
00598 SYSlerp(fpreal32 v1, fpreal32 v2, fpreal32 bias)
00599 {
00600     return v1 + (v2 - v1)*bias;
00601 }
00602 
00603 static inline fpreal64
00604 SYSlerp(fpreal64 v1, fpreal64 v2, fpreal64 bias)
00605 {
00606     return v1 + (v2 - v1)*bias;
00607 }
00608 
00609 static inline fpreal32
00610 SYSbilerp(fpreal32 v00, fpreal32 v10, fpreal32 v01, fpreal32 v11,
00611                         fpreal32 u, fpreal32 v)
00612 {
00613     return SYSlerp(SYSlerp(v00, v01, v), SYSlerp(v10, v11, v), u);
00614 }
00615 
00616 static inline fpreal64
00617 SYSbilerp(fpreal64 v00, fpreal64 v10, fpreal64 v01, fpreal64 v11,
00618                         fpreal64 u, fpreal64 v)
00619 {
00620     return SYSlerp(SYSlerp(v00, v01, v), SYSlerp(v10, v11, v), u);
00621 }
00622 
00623 static inline fpreal32
00624 SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 val)
00625 {
00626     fpreal32    t;
00627     if (val <= min) return 0;
00628     if (val >= max) return 1;
00629     t = max - min;
00630     if (SYSequalZero(t, (fpreal32)1e-8)) return (fpreal32).5;
00631     t = (val - min) / t;
00632     return t*t*((fpreal32)3.0 - (fpreal32)2.0*t);
00633 }
00634 
00635 static inline fpreal64
00636 SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 val)
00637 {
00638     fpreal64    t;
00639     if (val <= min) return 0;
00640     if (val >= max) return 1;
00641     t = max - min;
00642     if (SYSequalZero(t, (fpreal64)1e-18)) return (fpreal64).5;
00643     t = (val - min) / t;
00644     return t*t*((fpreal64)3.0 - (fpreal64)2.0*t);
00645 }
00646 
00647 static inline fpreal32
00648 SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 value, fpreal32 roll)
00649 {
00650     fpreal32    f;
00651     if (roll > 0)
00652     {
00653         f = SYSsmooth(min, max, value);
00654         return roll < (fpreal32)1 ? (fpreal32)1-SYSpow((fpreal32)1-f,
00655                                     (fpreal32)1/roll) : SYSpow(f, roll);
00656     }
00657     return (fpreal32)1;
00658 }
00659 
00660 static inline fpreal64
00661 SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 value, fpreal64 roll)
00662 {
00663     fpreal64    f;
00664     if (roll > 0)
00665     {
00666         f = SYSsmooth(min, max, value);
00667         return roll < (fpreal64)1 ? (fpreal64)1-SYSpow((fpreal64)1-f,
00668                                     (fpreal64)1/roll) : SYSpow(f, roll);
00669     }
00670     return (fpreal64)1;
00671 }
00672 
00673 static inline fpreal32
00674 SYSfit(fpreal32 val, fpreal32 omin, fpreal32 omax, fpreal32 nmin, fpreal32 nmax)
00675 {
00676     fpreal32    d = omax - omin;
00677     if (SYSequalZero(d, (fpreal32)1e-8))
00678         return (nmin+nmax)*(fpreal32).5;
00679     if (omin < omax)
00680     {
00681         if (val < omin) return nmin;
00682         if (val > omax) return nmax;
00683     }
00684     else
00685     {
00686         if (val < omax) return nmax;
00687         if (val > omin) return nmin;
00688     }
00689     return nmin + (nmax-nmin)*(val-omin)/d;
00690 }
00691 
00692 static inline fpreal64
00693 SYSfit(fpreal64 val, fpreal64 omin, fpreal64 omax, fpreal64 nmin, fpreal64 nmax)
00694 {
00695     fpreal64    d = omax - omin;
00696     if (SYSequalZero(d, (fpreal64)1e-18))
00697         return (nmin+nmax)*(fpreal64).5;
00698     if (omin < omax)
00699     {
00700         if (val < omin) return nmin;
00701         if (val > omax) return nmax;
00702     }
00703     else
00704     {
00705         if (val < omax) return nmax;
00706         if (val > omin) return nmin;
00707     }
00708     return nmin + (nmax-nmin)*(val-omin)/d;
00709 }
00710 
00711 static inline fpreal32
00712 SYSroundDownToMultipleOf(fpreal32 val, fpreal32 multiple)
00713 {
00714     fpreal32    modulus;
00715     fpreal32    retval;
00716     if (SYSequalZero(multiple))
00717         return val;
00718     modulus = SYSfmod(val, multiple);
00719     if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
00720         return val;
00721     retval = val - modulus;
00722     if( val < (fpreal32)0 && modulus!=(fpreal32)0 )
00723         retval -= multiple;
00724     return retval;
00725 }
00726 
00727 static inline fpreal64
00728 SYSroundDownToMultipleOf(fpreal64 val, fpreal64 multiple)
00729 {
00730     fpreal64    modulus;
00731     fpreal64    retval;
00732     if (SYSequalZero(multiple))
00733         return val;
00734     modulus = SYSfmod(val, multiple);
00735     if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
00736         return val;
00737     retval = val - modulus;
00738     if( val < (fpreal64)0 && modulus!=(fpreal64)0 )
00739         retval -= multiple;
00740     return retval;
00741 }
00742 
00743 static inline fpreal32
00744 SYSroundUpToMultipleOf(fpreal32 val, fpreal32 multiple)
00745 {
00746     fpreal32    modulus;
00747     if (SYSequalZero(multiple))
00748         return val;
00749     modulus = SYSfmod(val, multiple);
00750     if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
00751         return val;
00752     if (val > (fpreal32)0)
00753         val += multiple;
00754     return val - modulus;
00755 }
00756 
00757 static inline fpreal64
00758 SYSroundUpToMultipleOf(fpreal64 val, fpreal64 multiple)
00759 {
00760     fpreal64    modulus;
00761     if (SYSequalZero(multiple))
00762         return val;
00763     modulus = SYSfmod(val, multiple);
00764     if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
00765         return val;
00766     if (val > (fpreal64)0)
00767         val += multiple;
00768     return val - modulus;
00769 }
00770 
00771 static inline uint32
00772 SYSwang_inthash(uint32 key)
00773 {
00774     // From http://www.concentric.net/~Ttwang/tech/inthash.htm
00775     key += ~(key << 16);
00776     key ^=  (key >> 5);
00777     key +=  (key << 3);
00778     key ^=  (key >> 13);
00779     key += ~(key << 9);
00780     key ^=  (key >> 17);
00781     return key;
00782 }
00783 
00784 static inline uint32
00785 SYSmultiplicative_inthash(uint32 key)
00786 {
00787     // Multiply by the golden mean of 2^32 (Knuth's multiplicative method) to
00788     // get a uniformly distributed hash value.  Do a google search for
00789     // 2654435761 for more information.
00790     return key * 2654435761u;
00791 }
00792 
00793 static inline uint64
00794 SYSwang_inthash64(uint64 key)
00795 {
00796     // From http://www.concentric.net/~Ttwang/tech/inthash.htm
00797     key += ~(key << 32);
00798     key ^= (key >> 22);
00799     key += ~(key << 13);
00800     key ^= (key >> 8);
00801     key += (key << 3);
00802     key ^= (key >> 15);
00803     key += ~(key << 27);
00804     key ^= (key >> 31);
00805     return key;
00806 }
00807 
00808 static inline uint
00809 SYSreal_hash(fpreal32 a, int lowmask=0xf)
00810 {
00811     SYS_FPReal32Union   ai;
00812     ai.fval = a;
00813     return SYSwang_inthash(ai.uval & (~lowmask));
00814 }
00815 
00816 static inline uint
00817 SYSreal_hash(fpreal64 a, int lowmask=0xf)
00818 {
00819     SYS_FPReal64Union   ai;
00820     int32       lo, hi;
00821     ai.fval = a;
00822     lo = ai.uval & 0xfffffff;
00823     hi = ai.uval >> 32;
00824     return SYSwang_inthash(lo & (~lowmask)) ^ SYSwang_inthash(hi);
00825 }
00826 
00827 static inline uint
00828 SYSpointerHash(const void *ptr)
00829 {
00830 #if SIZEOF_VOID_P == 8
00831     // 64 bit pointers
00832     return (uint)((int64)ptr & 0xffffffff);
00833 #else
00834     // 32 bit pointers
00835     return (uint)ptr;
00836 #endif
00837 }
00838 
00839 // Generate a random number in range [0, 1)
00840 inline static fpreal32
00841 SYSfastRandom(uint &seed)
00842 {
00843     SYS_FPReal32Union   tmp;
00844     seed = seed*1664525 + 1013904223;
00845     tmp.uval = 0x3f800000 | (0x007fffff & seed);
00846     return tmp.fval-1.0F;
00847 }
00848 
00849 inline static fpreal32
00850 SYSrandom(uint &seed)
00851 {
00852     SYS_FPReal32Union   tmp;
00853     seed = seed*1664525 + 1013904223;
00854     tmp.uval = SYSwang_inthash(seed);
00855     tmp.uval = 0x3f800000 | (0x007fffff & tmp.uval);
00856     return tmp.fval - 1.0F;
00857 }
00858 
00859 inline static uint
00860 SYSfastRandomInt(uint &seed)
00861 {
00862     seed = seed*1664525 + 1013904223;
00863     return SYSwang_inthash(seed);
00864 }
00865 
00866 inline static fpreal32
00867 SYSfastRandomZero(uint &seed)
00868 {
00869     SYS_FPReal32Union   tmp;
00870     seed = seed*1664525 + 1013904223;
00871     tmp.uval = 0x3f800000 | (0x007fffff & seed);
00872     return tmp.fval - 1.5F;
00873 }
00874 
00875 inline static fpreal32
00876 SYSrandomZero(uint &seed)
00877 {
00878     SYS_FPReal32Union   tmp;
00879     seed = seed*1664525 + 1013904223;
00880     tmp.uval = SYSwang_inthash(seed);
00881     tmp.uval = 0x3f800000 | (0x007fffff & tmp.uval);
00882     return tmp.fval - 1.5F;
00883 }
00884 
00885 // Find the minimum/maximum of four floating point values using the minimum
00886 // number of comparisons.
00887 inline static void
00888 SYSminmax(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 v3,
00889          fpreal32 &min, fpreal32 &max)
00890 {
00891     if (v0 < v1) { min = v0; max = v1; }
00892     else         { min = v1; max = v0; }
00893     if (v2 < v3)
00894     {
00895         if (v2 < min) min = v2;
00896         if (v3 > max) max = v3;
00897     }
00898     else
00899     {
00900         if (v2 > max) max = v2;
00901         if (v3 < min) min = v3;
00902     }
00903 }
00904 
00905 inline static void
00906 SYSminmax(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 v3,
00907          fpreal64 &min, fpreal64 &max)
00908 {
00909     if (v0 < v1) { min = v0; max = v1; }
00910     else         { min = v1; max = v0; }
00911     if (v2 < v3)
00912     {
00913         if (v2 < min) min = v2;
00914         if (v3 > max) max = v3;
00915     }
00916     else
00917     {
00918         if (v2 > max) max = v2;
00919         if (v3 < min) min = v3;
00920     }
00921 }
00922 
00923 inline static void
00924 SYSminmax(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 &min, fpreal32 &max)
00925 {
00926     if (v0 < v1) { min = v0; max = v1; }
00927     else         { min = v1; max = v0; }
00928     if (v2 < min)       { min = v2; }
00929     else if (v2 > max)  { max = v2; }
00930 }
00931 
00932 inline static void
00933 SYSminmax(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 &min, fpreal64 &max)
00934 {
00935     if (v0 < v1) { min = v0; max = v1; }
00936     else         { min = v1; max = v0; }
00937     if (v2 < min)       { min = v2; }
00938     else if (v2 > max)  { max = v2; }
00939 }
00940 
00941 inline static void
00942 SYSgetSinCosFromSlope(fpreal32 slope, fpreal32 &sintheta, fpreal32 &costheta)
00943 {
00944     fpreal32    one_over_m;
00945     sintheta = slope / SYSsqrt(slope*slope + (fpreal32)1);
00946     if ((slope = SYSabs(slope)) > (fpreal32)1)
00947     {
00948         one_over_m = (fpreal32)1 / slope;
00949         costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
00950     }
00951     else
00952         costheta = SYSsqrt((fpreal32)1 - sintheta*sintheta);
00953 }
00954 
00955 inline static void
00956 SYSgetSinCosFromSlope(fpreal64 slope, fpreal64 &sintheta, fpreal64 &costheta)
00957 {
00958     fpreal64    one_over_m;
00959     sintheta = slope / SYSsqrt(slope*slope + (fpreal64)1);
00960     if ((slope = SYSabs(slope)) > (fpreal64)1)
00961     {
00962         one_over_m = (fpreal64)1 / slope;
00963         costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
00964     }
00965     else
00966         costheta = SYSsqrt((fpreal64)1 - sintheta*sintheta);
00967 }
00968 
00969 inline static bool
00970 SYSsameSign( fpreal32 v0, fpreal32 v1 )
00971 {
00972     return (v0*v1)>0;
00973 }
00974 
00975 inline static bool
00976 SYSsameSign( fpreal64 v0, fpreal64 v1 )
00977 {
00978     return (v0*v1)>0;
00979 }
00980 
00981 inline static bool
00982 SYSsameSign( int32 v0, int32 v1 )
00983 {
00984     return (v0 ^ v1) >= 0;
00985 }
00986 
00987 inline static bool
00988 SYSsameSign( int64 v0, int64 v1 )
00989 {
00990     return (v0 ^ v1) >= 0;
00991 }
00992 
00993 static inline uint      SYSnextPrime(uint num)  { return SYSmakePrime(num+1); }
00994 
00995 // Include permutations of fpreal32/fpreal64
00996 #include "SYS_MathPermute.h"
00997 #include "SYS_MathRestrictive.h"
00998 
00999 #else   /* cplusplus */
01000     #define SYSmax(a,b) ((a) > (b) ? (a) : (b))
01001     #define SYSmin(a,b) ((a) < (b) ? (a) : (b))
01002     #define SYSabs(a)   ((a) < 0 ? (a) : -(a))
01003 #endif
01004 
01005 #endif