HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SYS_Math.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: SYS_Math.h (SYS Library, C++)
7  *
8  * COMMENTS: Houdini math interface.
9  * statements like: a = sin(x); b = sinf(y);
10  * we use: a = SYSsin(x); b = SYSsin(y);
11  *
12  * Functions defined here:
13  * SYSmin(a,b) - Returns the minimum of a or b
14  * SYSmin(a,b,c) - Returns the minimum of a, b or c
15  * SYSmin(a,b,c,d) - Returns minimum of a, b, c or d
16  * SYSargmin(a,b) - Returns index of minimum number
17  * SYSargmin(a,b,c) - Returns index of minimum number
18  * SYSargmin(a,b,c,d) - Returns index of minimum number
19  * SYSmax(a,b) - Returns the maximum of a or b
20  * SYSmax(a,b,c) - Returns the maximum of a, b or c
21  * SYSmax(a,b,c,d) - Returns maximum of a, b, c or d
22  * SYSsort(a, b) - Sorts so a < b
23  * SYSsort(a, b, c) - Sorts so a < b < c
24  * SYSargmax(a,b) - Returns index of largest number
25  * SYSargmax(a,b,c) - Returns index of largest number
26  * SYSargmax(a,b,c,d) - Returns index of largest number
27  * SYSabs(x) - Works for all types (including fpreal)
28  * SYSavg(a,b,c) - Return average value of 3 values
29  * SYSavg(a,b,c,d) - Return average value of 4 values
30  *
31  * SYSminmax(a,b,&min.&max) - Find min and max of 2 values
32  * SYSminmax(a,b,c,&min.&max) - Find min and max of 3 values
33  * SYSminmax(a,b,c,d,&min.&max) - Find min and max of 4 values
34  *
35  * SYSwrapmod(a, b) - Returns [0..b) interval for +ve b.
36  *
37  * SYSsignum(a) - Numerically robust Sign function:
38  * -1, if a is negative
39  * 1, if a is positive
40  * 0, otherwise:
41  *
42  * Comparison:
43  * SYSequalZero(a, [tol]) - Is a equal to zero
44  * SYSisEqual(a, b, tol) -
45  * SYSalmostEqual(a, b, ulps) - Almost equal given units in last place
46  *
47  * Interpolation/Range:
48  * SYSclamp(a, min, max) - Clamp value between min and max
49  * SYSsmooth(min, max, a) - Ease-in ease-out curve
50  * SYSsmooth(min, max, a, roll) - Ease-in ease-out curve (with rolloff)
51  * SYSlerp(v1, v2, bias) - Linear interpolation
52  * SYSfit(v, omin,omax, nmin,nmax) - Map v in (omin,omax) to (nmin,nmax)
53  *
54  * The bilerp function expects:
55  * u0v1--u1v1
56  * | |
57  * u0v0--u1v0
58  * Where u is horizontal and v is vertical in the diagram above.
59  * SYSbilerp(u0v0, u1v0, u0v1, u1v1, u, v)
60  *
61  * Standard math (single/double precision signatures):
62  * SYSsin(x)
63  * SYScos(x)
64  * SYStan(x)
65  * SYSsqrt(x)
66  * SYSlog(x)
67  * SYSfabs(x)
68  * SYStrunc(x)
69  * SYSfloor(x)
70  * SYSceil(x)
71  * SYScbrt(x) -- Use #include <SYS/SYS_MathCbrt.h>
72  * SYSlog10(x)
73  * SYSfmod(x,y)
74  * SYSpow(x,y)
75  * SYSsafepow(x, y)
76  * SYSatan(x)
77  * SYSatan(y,x) -- Note, SYSatan(y,x) is equivalent to SYSatan2(y,x)
78  * SYSatan2(y,x)
79  * SYShypot(x,y)
80  * SYSsafediv(x,y) - Divide only if y != 0
81  * SYSsafesqrt(x) - Compute sqrt only if x >= 0
82  * SYSsafefmod(x, y) - Mod only if y != 0
83  *
84  * Random:
85  * SYSpointerHash - Generate a hash key for a pointer
86  * SYSwang_inthash - 32-bit Integer hashing function
87  * SYSwang_inthash64 - 64-bit Integer hashing function
88  * SYSreal_hash - Generate a hash for a real number
89  * SYSreal_hashseed - Generate a hash for a real number
90  * SYSvector_hash - Generate a hash for a real vector
91  * SYSfastRandom - Really fast random number generator (0,1)
92  * SYSrandom - Fast random number (fewer artifacts) (0,1)
93  * SYSfastRandomZero - Really fast random number generator (-.5, .5)
94  * SYSrandomZero - Fast random number (fewer artifacts) (-.5, .5)
95  *
96  * Rounding:
97  * SYSroundDownToMultipleOf -
98  * SYSroundUpToMultipleOf -
99  * SYSroundAngle -
100  *
101  * Misc:
102  * SYSisNan - Is this a valid floating-point number?
103  * SYSisFinite - Is not a NAN or Infinity
104  * SYSisInt - Is the string an integer
105  * SYSisFloat - Is the string a float
106  * SYSisPrime - Is it a prime?
107  * SYSsameSign(a,b) - Do the non-zero numbers have the same sign?
108  * (Note: does not work if either number is zero)
109  * SYSmakePrime - Make it a prime
110  * SYSnextPrime - Next prime
111  * SYSgetSinCosFromSlope - Compute sin/cos given a slope
112  * SYShexCharToInt - map '0'-'9,'a'-'z' to 0-15, or -1
113  * SYSintToHexChar - map 0-15 to '0'-'9,'a'-'z'
114  * SYSdivMod(n,d,q,r) - q = n/d; r = n%d; The compiler should make
115  * just one DIV or IDIV instruction.
116  *
117  * For information on almostEqual, please see:
118  * http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
119  *
120  */
121 
122 #ifndef __SYS_Math__
123 #define __SYS_Math__
124 
125 #include "SYS_API.h"
126 #include "SYS_Types.h"
127 
128 /*
129  * System dependent includes
130  */
131 #if defined(WIN32)
132 # include "SYS_NTMath.h"
133 #endif
134 
135 #if defined(LINUX)
136 # include "SYS_FastMath.h"
137 #endif
138 
139 #include <float.h>
140 #include <limits>
141 #include <math.h>
142 #include <stdio.h>
143 #include <stdlib.h>
144 
145 // We do not want to use the built-ins for srand48 or drand48 as they
146 // are not threadsafe. This version uses thread local storage to
147 // avoid that.
148 SYS_API void SYSsrand48(long seed);
149 SYS_API double SYSdrand48();
150 
151 #if !defined(SESI_ALLOW_DRAND48)
152 
153 // These macros are to aggressively prevent anyone from accidetnally
154 // using thread unsafe functions
155 #undef srand48
156 #undef drand48
157 #define srand48(X) static_assert(0, "Use SYSsrand48() instead")
158 #define drand48(X) 0; { static_assert(0, "Use SYSdrand48() instead"); }
159 
160 #endif
161 
162 // It would be useful to disable rand() and srand() as well, but
163 // they are used by std::rand :<
164 
165 template< typename F >
166 inline int SYSsignum( const F a )
167 {
168  return ( int( a > 0 ) - int( a < 0 ) );
169 }
170 
171 // Determine whether f is not-a-number (NAN) for any floating point type f
172 template< typename F > inline bool
173 SYSisNan(const F f)
174 {
175  return(f != f);
176 }
177 
178 template <> inline bool
180 {
181  return f.isNan();
182 }
183 
184 /// @{
185 /// Check whether a number is finite. That is, not Nan and not infinity.
186 /// @see SYSisNan()
187 #if defined(LINUX)
188 #include <cmath>
189 inline bool SYSisFinite(fpreal64 f) { return std::isfinite(f) != 0; }
190 inline bool SYSisFinite(fpreal32 f) { return std::isfinite(f) != 0; }
191 #elif defined(MBSD)
192 inline bool SYSisFinite(fpreal64 f) { return isfinite(f) != 0; }
193 inline bool SYSisFinite(fpreal32 f) { return isfinite(f) != 0; }
194 #elif defined(WIN32)
195 inline bool SYSisFinite(fpreal64 f) { return ::_finite(f) != 0; }
196 inline bool SYSisFinite(fpreal32 f) { return ::_finite(f) != 0; }
197 #endif
198 inline bool SYSisFinite(fpreal16 f) { return f.isFinite(); }
199 /// @}
200 
201 /// SYSisNan() checks whether the string represents a non-finite number
202 SYS_API bool SYSisNan(const char *number);
203 SYS_API bool SYSisInt(const char *str);
204 SYS_API bool SYSisFloat(const char *str);
209 SYS_API bool SYSisPrime(uint num);
213 
214 #if defined(__cplusplus)
215 
216 // NOTE:
217 // These have been carefully written so that in the case of equality
218 // we always return the first parameter. This is so that NANs in
219 // in the second parameter are suppressed.
220 #define h_min(a, b) (((a) > (b)) ? (b) : (a))
221 #define h_argmin(a, b) (((a) > (b)) ? 1 : 0)
222 #define h_max(a, b) (((a) < (b)) ? (b) : (a))
223 #define h_argmax(a, b) (((a) < (b)) ? 1 : 0)
224 // DO NOT CHANGE THE ABOVE WITHOUT READING THE COMMENT
225 #define h_abs(a) (((a) > 0) ? (a) : -(a))
226 #define h_sgn(a) (((a) > 0) ? 1 : (((a) < 0) ? -1 : 0))
227 
228 
229 static constexpr inline int16 SYSmin(int16 a, int16 b) { return h_min(a,b); }
230 static constexpr inline int16 SYSmax(int16 a, int16 b) { return h_max(a,b); }
231 static constexpr inline int16 SYSabs(int16 a) { return h_abs(a); }
232 static constexpr inline int16 SYSsgn(int16 a) { return h_sgn(a); }
233 static constexpr inline int SYSargmin(int16 a, int16 b) { return h_argmin(a,b);}
234 static constexpr inline int SYSargmax(int16 a, int16 b) { return h_argmax(a,b);}
235 static constexpr inline int32 SYSmin(int32 a, int32 b) { return h_min(a,b); }
236 static constexpr inline int32 SYSmax(int32 a, int32 b) { return h_max(a,b); }
237 static constexpr inline int32 SYSabs(int32 a) { return h_abs(a); }
238 static constexpr inline int32 SYSsgn(int32 a) { return h_sgn(a); }
239 static constexpr inline int SYSargmin(int32 a, int32 b) { return h_argmin(a,b);}
240 static constexpr inline int SYSargmax(int32 a, int32 b) { return h_argmax(a,b);}
241 static constexpr inline int64 SYSmin(int64 a, int64 b) { return h_min(a,b); }
242 static constexpr inline int64 SYSmax(int64 a, int64 b) { return h_max(a,b); }
243 static constexpr inline int64 SYSmin(int32 a, int64 b) { return h_min(a,b); }
244 static constexpr inline int64 SYSmax(int32 a, int64 b) { return h_max(a,b); }
245 static constexpr inline int64 SYSmin(int64 a, int32 b) { return h_min(a,b); }
246 static constexpr inline int64 SYSmax(int64 a, int32 b) { return h_max(a,b); }
247 static constexpr inline int64 SYSabs(int64 a) { return h_abs(a); }
248 static constexpr inline int64 SYSsgn(int64 a) { return h_sgn(a); }
249 static constexpr inline int SYSargmin(int64 a, int64 b) { return h_argmin(a,b);}
250 static constexpr inline int SYSargmax(int64 a, int64 b) { return h_argmax(a,b);}
251 static constexpr inline uint16 SYSmin(uint16 a, uint16 b) { return h_min(a,b); }
252 static constexpr inline uint16 SYSmax(uint16 a, uint16 b) { return h_max(a,b); }
253 static constexpr inline int SYSargmin(uint16 a, uint16 b) { return h_argmin(a,b);}
254 static constexpr inline int SYSargmax(uint16 a, uint16 b) { return h_argmax(a,b);}
255 static constexpr inline uint32 SYSmin(uint32 a, uint32 b) { return h_min(a,b); }
256 static constexpr inline uint32 SYSmax(uint32 a, uint32 b) { return h_max(a,b); }
257 static constexpr inline int SYSargmin(uint32 a, uint32 b) { return h_argmin(a,b);}
258 static constexpr inline int SYSargmax(uint32 a, uint32 b) { return h_argmax(a,b);}
259 static constexpr inline uint64 SYSmin(uint64 a, uint64 b) { return h_min(a,b); }
260 static constexpr inline uint64 SYSmax(uint64 a, uint64 b) { return h_max(a,b); }
261 static constexpr inline int SYSargmin(uint64 a, uint64 b) { return h_argmin(a,b);}
262 static constexpr inline int SYSargmax(uint64 a, uint64 b) { return h_argmax(a,b);}
263 static inline fpreal16 SYSmin(fpreal16 a, fpreal16 b) { return h_min(a,b); }
264 static inline fpreal16 SYSmax(fpreal16 a, fpreal16 b) { return h_max(a,b); }
265 static inline fpreal16 SYSsgn(fpreal16 a) { return h_sgn(a); }
266 static inline int SYSargmin(fpreal16 a, fpreal16 b){ return h_argmin(a,b);}
267 static inline int SYSargmax(fpreal16 a, fpreal16 b){ return h_argmax(a,b);}
268 static constexpr inline fpreal32 SYSmin(fpreal32 a, fpreal32 b) { return h_min(a,b); }
269 static constexpr inline fpreal32 SYSmax(fpreal32 a, fpreal32 b) { return h_max(a,b); }
270 static constexpr inline fpreal32 SYSsgn(fpreal32 a) { return h_sgn(a); }
271 static constexpr inline int SYSargmin(fpreal32 a, fpreal32 b){ return h_argmin(a,b);}
272 static constexpr inline int SYSargmax(fpreal32 a, fpreal32 b){ return h_argmax(a,b);}
273 static constexpr inline fpreal64 SYSmin(fpreal64 a, fpreal64 b) { return h_min(a,b); }
274 static constexpr inline fpreal64 SYSmax(fpreal64 a, fpreal64 b) { return h_max(a,b); }
275 static constexpr inline fpreal64 SYSsgn(fpreal64 a) { return h_sgn(a); }
276 static constexpr inline int SYSargmin(fpreal64 a, fpreal64 b){ return h_argmin(a,b);}
277 static constexpr inline int SYSargmax(fpreal64 a, fpreal64 b){ return h_argmax(a,b);}
278 
279 // Some systems have size_t as a seperate type from uint. Some don't.
280 #if (defined(LINUX) && defined(IA64)) || defined(MBSD)
281 static constexpr inline size_t SYSmin(size_t a, size_t b) { return h_min(a,b); }
282 static constexpr inline size_t SYSmax(size_t a, size_t b) { return h_max(a,b); }
283 static constexpr inline int SYSargmin(size_t a, size_t b) { return h_argmin(a,b);}
284 static constexpr inline int SYSargmax(size_t a, size_t b) { return h_argmax(a,b);}
285 #endif
286 
287 #undef h_min
288 #undef h_max
289 #undef h_abs
290 #undef h_sgn
291 
292 #define h_clamp(val, min, max, tol) \
293  ((val <= min+tol) ? min : ((val >= max-tol) ? max : val))
294 
295  static constexpr inline int
296  SYSclamp(int v, int min, int max)
297  { return h_clamp(v, min, max, 0); }
298 
299  static constexpr inline uint
300  SYSclamp(uint v, uint min, uint max)
301  { return h_clamp(v, min, max, 0); }
302 
303  static constexpr inline int64
304  SYSclamp(int64 v, int64 min, int64 max)
305  { return h_clamp(v, min, max, CONST_INT64(0)); }
306 
307  static constexpr inline uint64
308  SYSclamp(uint64 v, uint64 min, uint64 max)
309  { return h_clamp(v, min, max, CONST_UINT64(0)); }
310 
311  static constexpr inline fpreal32
312  SYSclamp(fpreal32 v, fpreal32 min, fpreal32 max, fpreal32 tol=(fpreal32)0)
313  { return h_clamp(v, min, max, tol); }
314 
315  static constexpr inline fpreal64
316  SYSclamp(fpreal64 v, fpreal64 min, fpreal64 max, fpreal64 tol=(fpreal64)0)
317  { return h_clamp(v, min, max, tol); }
318 
319 #undef h_clamp
320 
321 /// This converts from one integer type to another by clamping
322 /// the range, instead of the usual wrapping.
323 template<typename OUTTYPE,typename INTYPE>
324 OUTTYPE SYSclampInt(INTYPE value)
325 {
326  value = SYSclamp(value, (INTYPE)std::numeric_limits<OUTTYPE>::min(),
328  return OUTTYPE(value);
329 }
330 
331 // Wrapper for libm math function calls
332 #if defined(LINUX)
333 # define SYS_MF(X) SYS_FastMath::X
334 #else
335 # define SYS_MF(X) ::X
336 #endif
337 
338 #define SYS_UNARY(func) \
339  static inline fpreal64 SYS##func(fpreal64 arg) \
340  { return SYS_MF(func)(arg); } \
341  static inline fpreal32 SYS##func(fpreal32 arg) \
342  { return SYS_MF(func##f)(arg); } \
343  static inline fpreal64 SYS##func(int64 arg) \
344  { return SYS_MF(func)((fpreal64)arg); } \
345  static inline fpreal64 SYS##func(int32 arg) \
346  { return SYS_MF(func)((fpreal64)arg); }
347 #define SYS_BINARY(func) \
348  static inline fpreal64 SYS##func(fpreal64 a, fpreal64 b) \
349  { return SYS_MF(func)(a,b); } \
350  static inline fpreal32 SYS##func(fpreal32 a, fpreal32 b) \
351  { return SYS_MF(func##f)(a, b); }
352 
353 #if defined(WIN32)
354 #define hypotf(x,y) hypot((x),(y))
355 #endif
356 
357  SYS_UNARY(sin)
358  SYS_UNARY(cos)
359  SYS_UNARY(tan)
360  SYS_UNARY(sinh)
361  SYS_UNARY(cosh)
362  SYS_UNARY(tanh)
363  SYS_UNARY(sqrt)
364  SYS_UNARY(trunc)
365  SYS_UNARY(exp)
366  SYS_BINARY(fmod)
367  SYS_BINARY(pow)
368  SYS_BINARY(atan2)
369  SYS_BINARY(hypot)
370 
371  static inline fpreal32 SYSsafediv(fpreal32 x, fpreal32 y)
372  { return x/(y != 0 ? y : SYS_Types<fpreal32>::infinity()); }
373  static inline fpreal64 SYSsafediv(fpreal64 x, fpreal64 y)
374  { return x/(y != 0 ? y : SYS_Types<fpreal64>::infinity()); }
375  static inline fpreal32 SYSsafesqrt(fpreal32 x)
376  { return x > 0 ? SYSsqrt(x) : 0; }
377  static inline fpreal64 SYSsafesqrt(fpreal64 x)
378  { return x > 0 ? SYSsqrt(x) : 0; }
379  static inline fpreal32 SYSsafefmod(fpreal32 x, fpreal32 y)
380  { return y != 0 ? SYSfmod(x, y) : 0; }
381  static inline fpreal64 SYSsafefmod(fpreal64 x, fpreal64 y)
382  { return y != 0 ? SYSfmod(x, y) : 0; }
383 
384 #if 0
385 #include <xmmintrin.h>
386  /// Computes 1/sqrt(x) to about 11.4 bits of accuracy
387  static inline fpreal32 SYSrsqrt11(fpreal32 x)
388  {
389 #if defined(WIN32) || defined(LINUX) || defined(MBSD)
390  union {
391  fpreal32 val;
392  __m128 vec;
393  };
394  val = x;
395  vec = _mm_rsqrt_ss(vec);
396  return val;
397 #else
398  return 1.0f/SYSsqrt(x);
399 #endif
400  }
401 
402  /// Computes 1/sqrt(x) to about 22.2 bits of accuracy
403  static inline fpreal32 SYSrsqrt22(fpreal32 x)
404  {
405 #if defined(WIN32) || defined(LINUX) || defined(MBSD)
406  fpreal32 approximate = SYSrsqrt11(x);
407  // -(((1+e)/sqrt(x))^2 * x - 3) * ((1+e)/sqrt(x)) / 2
408  // -(-1 + e + (e^2)/2) * ((1+e)/sqrt(x))
409  // -(-1 + e + (e^2)/2 - e + e^2 + (e^3)/2) / sqrt(x)
410  // (1 - 3(e^2)/2 - (e^3)/2) / sqrt(x);
411  fpreal32 result = (approximate*approximate*x - 3)*approximate*-0.5f;
412 #else
413  fpreal32 result = 1.0f/SYSsqrt(x);
414 #endif
415  return result;
416  }
417 #endif
418 
419  static inline fpreal32 SYSlog(fpreal32 v) { return v <= 0 ? 0 : SYS_MF(logf)(v); }
420  static inline fpreal64 SYSlog(fpreal64 v) { return v <= 0 ? 0 : SYS_MF(log)(v); }
421  static inline fpreal32 SYSlog10(fpreal32 v) { return v <= 0 ? 0 : SYS_MF(log10f)(v); }
422  static inline fpreal64 SYSlog10(fpreal64 v) { return v <= 0 ? 0 : SYS_MF(log10)(v); }
423 
424 #if defined(WIN32)
425  static inline fpreal32 SYSexpm1(fpreal32 x) { return SYSexp(x) - 1; }
426  static inline fpreal64 SYSexpm1(fpreal64 x) { return SYSexp(x) - 1; }
427  static inline fpreal32 SYSlog1p(fpreal32 x) { return SYSlog(x+1); }
428  static inline fpreal64 SYSlog1p(fpreal64 x) { return SYSlog(x+1); }
429 #else
430  SYS_UNARY(expm1)
431  SYS_UNARY(log1p)
432 #endif
433 
434 #undef SYS_UNARY
435 #undef SYS_BINARY
436 #undef hypotf
437 
438 static inline fpreal32 SYSabs(fpreal32 a) { return SYS_MF(fabsf)(a); }
439 static inline fpreal64 SYSabs(fpreal64 a) { return SYS_MF(fabs)(a); }
440 static inline fpreal32 SYSfabs(fpreal32 a) { return SYS_MF(fabsf)(a); }
441 static inline fpreal64 SYSfabs(fpreal64 a) { return SYS_MF(fabs)(a); }
442 
443 #include "SYS_Floor.h"
444 
445 static inline fpreal32 SYSasin(fpreal32 a)
446  { return SYS_MF(asinf)(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
447 static inline fpreal32 SYSacos(fpreal32 a)
448  { return SYS_MF(acosf)(SYSclamp(a, (fpreal32)-1, (fpreal32)1)); }
449 static inline fpreal32 SYSatan(fpreal32 a)
450  { return SYS_MF(atanf)(a); }
451 static inline fpreal32 SYSatan(fpreal32 y, fpreal32 x)
452  { return SYS_MF(atan2f)(y, x); }
453 static inline fpreal64 SYSasin(fpreal64 a)
454  { return SYS_MF(asin) (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
455 static inline fpreal64 SYSacos(fpreal64 a)
456  { return SYS_MF(acos) (SYSclamp(a, (fpreal64)-1, (fpreal64)1)); }
457 static inline fpreal64 SYSatan(fpreal64 a)
458  { return SYS_MF(atan) (a); }
459 static inline fpreal64 SYSatan(fpreal64 y, fpreal64 x)
460  { return SYS_MF(atan2)(y, x); }
461 
462 static inline fpreal32 SYSsafepow(fpreal32 x, fpreal32 y)
463  { return SYSpow(x, x < 0 ? SYSrint(y) : y); }
464 static inline fpreal64 SYSsafepow(fpreal64 x, fpreal64 y)
465  { return SYSpow(x, x < 0 ? SYSrint(y) : y); }
466 
467 static inline fpreal32 SYSdegToRad(fpreal32 a) { return a*(fpreal32)(M_PI/180.0); }
468 static inline fpreal64 SYSdegToRad(fpreal64 a) { return a*(fpreal64)(M_PI/180.0); }
469 static inline fpreal64 SYSdegToRad(int32 a) { return a*(fpreal64)(M_PI/180.0); }
470 static inline fpreal64 SYSdegToRad(int64 a) { return a*(fpreal64)(M_PI/180.0); }
471 
472 static inline fpreal32 SYSradToDeg(fpreal32 a) { return a*(fpreal32)(180.0/M_PI); }
473 static inline fpreal64 SYSradToDeg(fpreal64 a) { return a*(fpreal64)(180.0/M_PI); }
474 static inline fpreal64 SYSradToDeg(int32 a) { return a*(fpreal32)(180.0/M_PI); }
475 static inline fpreal64 SYSradToDeg(int64 a) { return a*(fpreal64)(180.0/M_PI); }
476 
477 #define h_compare(func, code) \
478  static inline bool func(fpreal32 a, fpreal32 b, \
479  fpreal32 tol=SYS_FTOLERANCE) \
480  { return code; } \
481  static inline bool func(fpreal64 a, fpreal64 b, \
482  fpreal64 tol=SYS_FTOLERANCE_D) \
483  { return code; }
484 
485  static inline bool
486  SYSequalZero(fpreal32 a, fpreal32 tol=SYS_FTOLERANCE)
487  { return a >= -tol && a <= tol; }
488 
489  static inline bool
490  SYSequalZero(fpreal64 a, fpreal64 tol=SYS_FTOLERANCE_D)
491  { return a >= -tol && a <= tol; }
492 
493  static inline bool
494  SYSequalZero(int32 a, fpreal64 tol=SYS_FTOLERANCE_D)
495  { return a >= -tol && a <= tol; }
496 
497  static inline bool
498  SYSequalZero(int64 a, fpreal64 tol=SYS_FTOLERANCE_D)
499  { return a >= -tol && a <= tol; }
500 
501 
502  h_compare(SYSisEqual, SYSabs(a-b)<=tol)
503  h_compare(SYSisGreater, (a-b) > tol)
504  h_compare(SYSisGreaterOrEqual, (a-b) >= -tol)
505  h_compare(SYSisLess, (a-b) < -tol)
506  h_compare(SYSisLessOrEqual, (a-b) <= tol)
507 
508 inline bool SYSisEqual(int32 a, int32 b) { return a == b; }
509 inline bool SYSisEqual(int64 a, int64 b) { return a == b; }
510 inline bool SYSisGreater(int32 a, int32 b) { return a > b; }
511 inline bool SYSisGreater(int64 a, int64 b) { return a > b; }
512 inline bool SYSisGreaterOrEqual(int32 a, int32 b) { return a >= b; }
513 inline bool SYSisGreaterOrEqual(int64 a, int64 b) { return a >= b; }
514 inline bool SYSisLess(int32 a, int32 b) { return a < b; }
515 inline bool SYSisLess(int64 a, int64 b) { return a < b; }
516 inline bool SYSisLessOrEqual(int32 a, int32 b) { return a <= b; }
517 inline bool SYSisLessOrEqual(int64 a, int64 b) { return a <= b; }
518 
519 #undef h_compare
520 
521 static inline bool
522 SYSalmostEqual(fpreal32 a, fpreal32 b, int ulps=50)
523 {
524  SYS_FPRealUnionF ai, bi;
525  ai.fval = a;
526  bi.fval = b;
527  if (ai.ival < 0)
528  ai.ival = 0x80000000 - ai.ival;
529  if (bi.ival < 0)
530  bi.ival = 0x80000000 - bi.ival;
531  // We're very careful to avoid signed integer overflow here.
532  if ((ai.uval > bi.uval ? ai.uval - bi.uval : bi.uval - ai.uval) <= ulps)
533  return true;
534  // If both single-precision floats are very close to zero, we consider them
535  // equal.
536  return SYSabs(a) < 1e-6 && SYSabs(b) < 1e-6;
537 }
538 
539 static inline bool
540 SYSalmostEqual(fpreal64 a, fpreal64 b, int64 ulps)
541 {
542  SYS_FPRealUnionD ai, bi;
543  ai.fval = a;
544  bi.fval = b;
545  if (ai.ival < 0)
546  ai.ival = CONST_INT64(0x8000000000000000) - ai.ival;
547  if (bi.ival < 0)
548  bi.ival = CONST_INT64(0x8000000000000000) - bi.ival;
549  // We're very careful to avoid signed integer overflow here.
550  if ((ai.uval > bi.uval ? ai.uval - bi.uval : bi.uval - ai.uval) <= ulps)
551  return true;
552  // If both double-precision floats are very close to zero, we consider them
553  // equal.
554  return SYSabs(a) < 1e-15 && SYSabs(b) < 1e-15;
555 }
556 
557 static inline bool
558 SYSalmostEqual(fpreal64 a, fpreal64 b, int32 ulps=50)
559 {
560  return SYSalmostEqual(a, b, (int64)ulps);
561 }
562 
563 #define h_max3(type) \
564  static inline type \
565  SYSmax(type v0, type v1, type v2) { \
566  return SYSmax(v2, SYSmax(v0, v1)); \
567  }
568 #define h_max4(type) \
569  static inline type \
570  SYSmax(type v0, type v1, type v2, type v3) { \
571  return SYSmax(SYSmax(v0, v1), SYSmax(v2, v3)); \
572  }
573 #define h_argmax3(type) \
574  static inline int \
575  SYSargmax(type v0, type v1, type v2) { \
576  return v2 > SYSmax(v0, v1) ? 2 : SYSargmax(v0, v1); \
577  }
578 #define h_argmax4(type) \
579  static inline int \
580  SYSargmax(type v0, type v1, type v2, type v3) { \
581  return SYSmax(v0, v1) < SYSmax(v2, v3) ? \
582  (SYSargmax(v2, v3) + 2) : SYSargmax(v0, v1); \
583  }
584 #define h_min3(type) \
585  static inline type \
586  SYSmin(type v0, type v1, type v2) { \
587  return SYSmin(v2, SYSmin(v0, v1)); \
588  }
589 #define h_min4(type) \
590  static inline type \
591  SYSmin(type v0, type v1, type v2, type v3) { \
592  return SYSmin(SYSmin(v0, v1), SYSmin(v2, v3)); \
593  }
594 #define h_argmin3(type) \
595  static inline int \
596  SYSargmin(type v0, type v1, type v2) { \
597  return v2 < SYSmin(v0, v1) ? 2 : SYSargmin(v0, v1); \
598  }
599 #define h_argmin4(type) \
600  static inline int \
601  SYSargmin(type v0, type v1, type v2, type v3) { \
602  return SYSmin(v0, v1) > SYSmin(v2, v3) ? \
603  (SYSargmin(v2, v3) + 2) : SYSargmin(v0, v1); \
604  }
605 
606 
607 #define h_max(type) h_min3(type) h_min4(type) h_max3(type) h_max4(type) \
608  h_argmin3(type) h_argmin4(type) \
609  h_argmax3(type) h_argmax4(type)
610 
611 #define h_avg3(type) \
612  static inline type \
613  SYSavg(type v0, type v1, type v2) { \
614  return (v0+v1+v2) * ((type)(1.0/3.0)); \
615  }
616 #define h_avg4(type) \
617  static inline type \
618  SYSavg(type v0, type v1, type v2, type v3) { \
619  return (v0+v1+v2+v3) * ((type)0.25); \
620  }
621 
622 #define h_avg(type) h_avg3(type) h_avg4(type)
623 
624 h_max(int8)
625 h_max(uint8)
626 h_max(int16)
627 h_max(uint16)
628 h_max(int32)
629 h_max(uint32)
630 h_max(int64)
631 h_max(uint64)
632 h_max(fpreal32) h_avg(fpreal32)
633 h_max(fpreal64) h_avg(fpreal64)
634 
635 static inline int32
636 SYSavg(int32 a, int32 b, int32 c)
637 { return (a + b + c + 1) / 3; }
638 static inline int32
639 SYSavg(int32 a, int32 b, int32 c, int32 d)
640 { return (a + b + c + d + 2) / 4; }
641 
642 static inline int64
643 SYSavg(int64 a, int64 b, int64 c)
644 { return (a + b + c + 1) / 3; }
645 static inline int64
646 SYSavg(int64 a, int64 b, int64 c, int64 d)
647 { return (a + b + c + d + 2) / 4; }
648 
649 // Some systems have size_t as a seperate type from uint. Some don't.
650 #if (defined(LINUX) && defined(IA64)) || defined(MBSD)
651 h_max(size_t)
652 #endif
653 
654 #undef h_min3
655 #undef h_min4
656 #undef h_max3
657 #undef h_max4
658 #undef h_argmin3
659 #undef h_argmin4
660 #undef h_argmax3
661 #undef h_argmax4
662 #undef h_avg3
663 #undef h_avg4
664 #undef h_max
665 #undef h_avg
666 
667 /// @{
668 /// Linear interpolation between v0 and v1.
669 static inline fpreal32
671 {
672  return v1 + (v2 - v1)*t;
673 }
674 
675 static inline fpreal64
677 {
678  return v1 + (v2 - v1)*t;
679 }
680 /// @}
681 
682 /// @{
683 /// Bilinear interpolation over a quadrilateral: @code
684 /// (u=0,v=1) u0v1 u1v1 (u=1,v=1)
685 /// +--------+
686 /// | |
687 /// | |
688 /// +--------+
689 /// (u=0,v=0) u0v0 u1v0 (u=1,v=0)
690 /// @endcode
691 static inline fpreal32
692 SYSbilerp(fpreal32 u0v0, fpreal32 u1v0, fpreal32 u0v1, fpreal32 u1v1,
693  fpreal32 u, fpreal32 v)
694 {
695  return SYSlerp(SYSlerp(u0v0, u0v1, v), SYSlerp(u1v0, u1v1, v), u);
696 }
697 
698 static inline fpreal64
699 SYSbilerp(fpreal64 u0v0, fpreal64 u1v0, fpreal64 u0v1, fpreal64 u1v1,
700  fpreal64 u, fpreal64 v)
701 {
702  return SYSlerp(SYSlerp(u0v0, u0v1, v), SYSlerp(u1v0, u1v1, v), u);
703 }
704 /// @}
705 
706 /// @{
707 /// Barycentric coordinates do interpolation over the triangle specified by the
708 /// three vertices (v0, v1, v2): @code
709 /// | (u=0,v=1) |
710 /// | v2 |
711 /// | / \ |
712 /// | / \ |
713 /// | / \ |
714 /// | (u=0,v=0) v0------v1 (u=1,v=0) |
715 /// @endcode
716 static inline fpreal32
717 SYSbarycentric(fpreal32 v0, fpreal32 v1, fpreal32 v2, fpreal32 u, fpreal32 v)
718 {
719  return v0*(1-u-v) + v1*u + v2*v;
720 }
721 
722 static inline fpreal64
723 SYSbarycentric(fpreal64 v0, fpreal64 v1, fpreal64 v2, fpreal64 u, fpreal64 v)
724 {
725  return v0*(1-u-v) + v1*u + v2*v;
726 }
727 /// @}
728 
729 /// @{
730 /// Wrap mod is like mod, but returns [0..b) interval thereby
731 /// usually giving a more useful result.
732 /// This is always defined as safe!
733 /// @}
734 static inline fpreal32
735 SYSwrapmod(fpreal32 a, fpreal32 b)
736 {
737  fpreal32 r = SYSsafefmod(a, b);
738  return ((a < 0) ^ (b < 0) && r) ? r+b : r;
739 }
740 
741 static inline fpreal64
742 SYSwrapmod(fpreal64 a, fpreal64 b)
743 {
744  fpreal64 r = SYSsafefmod(a, b);
745  return ((a < 0) ^ (b < 0) && r) ? r+b : r;
746 }
747 
748 static inline int32
749 SYSwrapmod(int32 a, int32 b)
750 {
751  if (b == 0)
752  return 0;
753  int32 r = a % b;
754  return ((a < 0) ^ (b < 0) && r) ? r+b : r;
755 }
756 
757 static inline int64
758 SYSwrapmod(int64 a, int64 b)
759 {
760  if (b == 0)
761  return 0;
762  int64 r = a % b;
763  return ((a < 0) ^ (b < 0) && r) ? r+b : r;
764 }
765 
766 static inline fpreal32
767 SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 val)
768 {
769  fpreal32 t;
770  if (val <= min) return 0;
771  if (val >= max) return 1;
772  t = max - min;
773  if (SYSequalZero(t, (fpreal32)1e-8)) return (fpreal32).5;
774  t = (val - min) / t;
775  return t*t*((fpreal32)3.0 - (fpreal32)2.0*t);
776 }
777 
778 static inline fpreal64
779 SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 val)
780 {
781  fpreal64 t;
782  if (val <= min) return 0;
783  if (val >= max) return 1;
784  t = max - min;
785  if (SYSequalZero(t, (fpreal64)1e-18)) return (fpreal64).5;
786  t = (val - min) / t;
787  return t*t*((fpreal64)3.0 - (fpreal64)2.0*t);
788 }
789 
790 static inline fpreal32
791 SYSsmooth(fpreal32 min, fpreal32 max, fpreal32 value, fpreal32 roll)
792 {
793  fpreal32 f;
794  if (roll > 0)
795  {
796  f = SYSsmooth(min, max, value);
797  return roll < (fpreal32)1 ? (fpreal32)1-SYSpow((fpreal32)1-f,
798  (fpreal32)1/roll) : SYSpow(f, roll);
799  }
800  return (fpreal32)1;
801 }
802 
803 static inline fpreal64
804 SYSsmooth(fpreal64 min, fpreal64 max, fpreal64 value, fpreal64 roll)
805 {
806  fpreal64 f;
807  if (roll > 0)
808  {
809  f = SYSsmooth(min, max, value);
810  return roll < (fpreal64)1 ? (fpreal64)1-SYSpow((fpreal64)1-f,
811  (fpreal64)1/roll) : SYSpow(f, roll);
812  }
813  return (fpreal64)1;
814 }
815 
816 static inline fpreal32
817 SYSfit(fpreal32 val, fpreal32 omin, fpreal32 omax, fpreal32 nmin, fpreal32 nmax)
818 {
819  fpreal32 d = omax - omin;
820  if (SYSequalZero(d, (fpreal32)1e-8))
821  return (nmin+nmax)*(fpreal32).5;
822  if (omin < omax)
823  {
824  if (val < omin) return nmin;
825  if (val > omax) return nmax;
826  }
827  else
828  {
829  if (val < omax) return nmax;
830  if (val > omin) return nmin;
831  }
832  return nmin + (nmax-nmin)*(val-omin)/d;
833 }
834 
835 static inline fpreal64
836 SYSfit(fpreal64 val, fpreal64 omin, fpreal64 omax, fpreal64 nmin, fpreal64 nmax)
837 {
838  fpreal64 d = omax - omin;
839  if (SYSequalZero(d, (fpreal64)1e-18))
840  return (nmin+nmax)*(fpreal64).5;
841  if (omin < omax)
842  {
843  if (val < omin) return nmin;
844  if (val > omax) return nmax;
845  }
846  else
847  {
848  if (val < omax) return nmax;
849  if (val > omin) return nmin;
850  }
851  return nmin + (nmax-nmin)*(val-omin)/d;
852 }
853 
854 /// SYSefit() will not clamp the values to the range
855 static inline fpreal32
856 SYSefit(fpreal32 v, fpreal32 omin, fpreal32 omax, fpreal32 nmin, fpreal32 nmax)
857 {
858  fpreal32 d = omax - omin;
859  if (SYSequalZero(d, (fpreal32)1e-8))
860  return (nmin+nmax)*0.5f;
861  return nmin + (nmax-nmin)*(v-omin)/d;
862 }
863 
864 static inline fpreal64
865 SYSefit(fpreal64 v, fpreal64 omin, fpreal64 omax, fpreal64 nmin, fpreal64 nmax)
866 {
867  fpreal64 d = omax - omin;
868  if (SYSequalZero(d, (fpreal64)1e-18))
869  return (nmin+nmax)*0.5f;
870  return nmin + (nmax-nmin)*(v-omin)/d;
871 }
872 
873 static inline fpreal32
874 SYSfit01(fpreal32 val, fpreal32 nmin, fpreal32 nmax)
875 {
876  if (val < 0) return nmin;
877  if (val > 1) return nmax;
878  return SYSlerp(nmin, nmax, val);
879 }
880 
881 static inline fpreal64
882 SYSfit01(fpreal64 val, fpreal64 nmin, fpreal64 nmax)
883 {
884  if (val < 0) return nmin;
885  if (val > 1) return nmax;
886  return SYSlerp(nmin, nmax, val);
887 }
888 
889 // Linear hat function over kernel width dx.
890 static inline fpreal32
891 SYShat(fpreal32 x, fpreal32 dx)
892 {
893  const fpreal32 ax = SYSabs(x);
894  return (ax > dx ? 0 : 1 - ax / dx);
895 }
896 
897 static inline fpreal64
898 SYShat(fpreal64 x, fpreal64 dx)
899 {
900  const fpreal64 ax = SYSabs(x);
901  return (ax > dx ? 0 : 1 - ax / dx);
902 }
903 
904 // Derivative of the linear hat function over kernel width dx.
905 static inline fpreal32
906 SYSdhat(fpreal32 x, fpreal32 dx)
907 {
908  return (x < -dx || x > dx ? 0 : -SYSsgn(x) / dx);
909 }
910 
911 static inline fpreal64
912 SYSdhat(fpreal64 x, fpreal64 dx)
913 {
914  return (x < -dx || x > dx ? 0 : -SYSsgn(x) / dx);
915 }
916 
917 
918 // For integer types.
919 template<typename T>
920 inline T
921 SYSroundDownToMultipleOf(T val, T multiple)
922 {
923  // Only handle multiples of 2 and up.
924  if (multiple <= 1)
925  return val;
926 
927  int rem = val % multiple;
928  if (rem == 0)
929  return val;
930 
931  if (val < 0)
932  return val - multiple - rem;
933  else
934  return val - rem;
935 }
936 
937 template<>
938 inline fpreal32
939 SYSroundDownToMultipleOf<fpreal32>(fpreal32 val, fpreal32 multiple)
940 {
941  fpreal32 modulus;
942  fpreal32 retval;
943  if (SYSequalZero(multiple))
944  return val;
945  modulus = SYSfmod(val, multiple);
946  if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
947  return val;
948  retval = val - modulus;
949  if( val < (fpreal32)0 && modulus!=(fpreal32)0 )
950  retval -= multiple;
951  return retval;
952 }
953 
954 template<>
955 inline fpreal64
956 SYSroundDownToMultipleOf<fpreal64>(fpreal64 val, fpreal64 multiple)
957 {
958  fpreal64 modulus;
959  fpreal64 retval;
960  if (SYSequalZero(multiple))
961  return val;
962  modulus = SYSfmod(val, multiple);
963  if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
964  return val;
965  retval = val - modulus;
966  if( val < (fpreal64)0 && modulus!=(fpreal64)0 )
967  retval -= multiple;
968  return retval;
969 }
970 
971 template<typename T>
972 inline T
973 SYSroundUpToMultipleOf(T val, T multiple)
974 {
975  // Only handle multiples of 2 and up.
976  if (multiple <= 1)
977  return val;
978 
979  int rem = val % multiple;
980  if (rem == 0)
981  return val;
982 
983  if (val < 0)
984  return val - rem;
985  else
986  return val + multiple - rem;
987 }
988 
989 template<>
990 inline fpreal32
991 SYSroundUpToMultipleOf<fpreal32>(fpreal32 val, fpreal32 multiple)
992 {
993  fpreal32 modulus;
994  if (SYSequalZero(multiple))
995  return val;
996  modulus = SYSfmod(val, multiple);
997  if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
998  return val;
999  if (val > (fpreal32)0)
1000  val += multiple;
1001  return val - modulus;
1002 }
1003 
1004 template<>
1005 inline fpreal64
1006 SYSroundUpToMultipleOf<fpreal64>(fpreal64 val, fpreal64 multiple)
1007 {
1008  fpreal64 modulus;
1009  if (SYSequalZero(multiple))
1010  return val;
1011  modulus = SYSfmod(val, multiple);
1012  if (SYSequalZero(modulus) || SYSequalZero(modulus-multiple))
1013  return val;
1014  if (val > (fpreal64)0)
1015  val += multiple;
1016  return val - modulus;
1017 }
1018 
1019 inline uint32
1020 SYSroundUpPow2(uint32 val)
1021 {
1022  // From http://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
1023  uint32 v = val;
1024 
1025  v--;
1026  v |= v >> 1;
1027  v |= v >> 2;
1028  v |= v >> 4;
1029  v |= v >> 8;
1030  v |= v >> 16;
1031  v++;
1032 
1033  return v;
1034 }
1035 
1036 static inline uint32
1037 SYSwang_inthash(uint32 key)
1038 {
1039  // From http://www.concentric.net/~Ttwang/tech/inthash.htm
1040  key += ~(key << 16);
1041  key ^= (key >> 5);
1042  key += (key << 3);
1043  key ^= (key >> 13);
1044  key += ~(key << 9);
1045  key ^= (key >> 17);
1046  return key;
1047 }
1048 
1049 static inline uint32
1050 SYSmultiplicative_inthash(uint32 key)
1051 {
1052  // Multiply by the golden mean of 2^32 (Knuth's multiplicative method) to
1053  // get a uniformly distributed hash value. Do a google search for
1054  // 2654435761 for more information.
1055  return key * 2654435761u;
1056 }
1057 
1058 static inline uint64
1059 SYSwang_inthash64(uint64 key)
1060 {
1061  // From http://www.concentric.net/~Ttwang/tech/inthash.htm
1062  key += ~(key << 32ULL);
1063  key ^= (key >> 22);
1064  key += ~(key << 13);
1065  key ^= (key >> 8);
1066  key += (key << 3);
1067  key ^= (key >> 15);
1068  key += ~(key << 27);
1069  key ^= (key >> 31);
1070  return key;
1071 }
1072 
1073 static inline uint
1074 SYSreal_hash(fpreal16 a, int lowmask=0x3)
1075 {
1076  return SYSwang_inthash(a.bits() & (~lowmask));
1077 }
1078 
1079 static inline uint
1080 SYSreal_hashseed(fpreal16 a, uint seed, int lowmask=0x3)
1081 {
1082  return SYSwang_inthash(seed + (a.bits() & (~lowmask)));
1083 }
1084 
1085 static inline uint
1086 SYSreal_hash(fpreal32 a, int lowmask=0xf)
1087 {
1088  SYS_FPRealUnionF ai;
1089  ai.fval = a;
1090  return SYSwang_inthash(ai.uval & (~lowmask));
1091 }
1092 
1093 static inline uint
1094 SYSreal_hashseed(fpreal32 a, uint seed, int lowmask=0xf)
1095 {
1096  SYS_FPRealUnionF ai;
1097  ai.fval = a;
1098  return SYSwang_inthash(seed + (ai.uval & (~lowmask)));
1099 }
1100 
1101 static inline uint
1102 SYSreal_hash(fpreal64 a, int lowmask=0xf)
1103 {
1104  SYS_FPRealUnionD ai;
1105  ai.fval = a;
1106  return SYSwang_inthash64(ai.uval & (~lowmask));
1107 }
1108 
1109 static inline uint
1110 SYSreal_hashseed(fpreal64 a, uint seed, int64 lowmask=0xf)
1111 {
1112  SYS_FPRealUnionD ai;
1113  ai.fval = a;
1114  return SYSwang_inthash64(seed + (ai.uval & (~lowmask)));
1115 }
1116 
1117 static inline uint
1118 SYSvector_hash(const int32 *vector, int size)
1119 {
1120  uint hash = 0;
1121  for (int i = 0; i < size; ++i)
1122  hash = SYSwang_inthash(hash + vector[i]);
1123  return hash;
1124 }
1125 
1126 static inline uint
1127 SYSvector_hash(const int64 *vector, int size)
1128 {
1129  uint hash = 0;
1130  for (int i = 0; i < size; ++i)
1131  hash = (uint)SYSwang_inthash64(hash + vector[i]);
1132  return hash;
1133 }
1134 
1135 static inline uint
1136 SYSvector_hash(const fpreal16 *vector, int size)
1137 {
1138  uint hash = 0;
1139  for (int i = 0; i < size; ++i)
1140  hash = SYSreal_hashseed(vector[i], hash);
1141  return hash;
1142 }
1143 
1144 static inline uint
1145 SYSvector_hash(const fpreal32 *vector, int size)
1146 {
1147  uint hash = 0;
1148  for (int i = 0; i < size; ++i)
1149  hash = SYSreal_hashseed(vector[i], hash);
1150  return hash;
1151 }
1152 
1153 static inline uint
1154 SYSvector_hash(const fpreal64 *vector, int size)
1155 {
1156  uint hash = 0;
1157  for (int i = 0; i < size; ++i)
1158  hash = SYSreal_hashseed(vector[i], hash);
1159  return hash;
1160 }
1161 
1162 // Spatial hashing function for 3-vectors
1163 // Reference:
1164 // "Optimized Spatial Hashing for Collision Detection of Deformable Objects"
1165 //
1166 template<class P>
1167 static inline size_t
1168 SYSvector3_hash(const P &vector)
1169 {
1170  static const size_t p1 = 73856093;
1171  static const size_t p2 = 19349663;
1172  static const size_t p3 = 83492791;
1173 
1174  return size_t(vector.x()*p1) ^ size_t(vector.y()*p2) ^ size_t(vector.z()*p3);
1175 }
1176 
1177 static inline uint
1178 SYSpointerHash(const void *ptr)
1179 {
1180 #if SIZEOF_VOID_P == 8
1181  // 64 bit pointers
1182  return (uint)((int64)ptr & 0xffffffff);
1183 #else
1184  // 32 bit pointers
1185  return (uint)ptr;
1186 #endif
1187 }
1188 
1189 // Convert a uniform random bit pattern to a random float in the range [0, 1)
1190 inline static fpreal32
1191 SYShashToFloat01(uint hash)
1192 {
1193  SYS_FPRealUnionF tmp;
1194  tmp.uval = 0x3f800000 | (0x007fffff & hash);
1195  return tmp.fval-1.0F;
1196 }
1197 
1198 // Generate a random number in range [0, 1)
1199 inline static fpreal32
1200 SYSfastRandom(uint &seed)
1201 {
1202  seed = seed*1664525 + 1013904223;
1203  return SYShashToFloat01(seed);
1204 }
1205 
1206 inline static fpreal32
1207 SYSrandom(uint &seed)
1208 {
1209  seed = seed*1664525 + 1013904223;
1210  return SYShashToFloat01(SYSwang_inthash(seed));
1211 }
1212 
1213 inline static uint
1214 SYSfastRandomInt(uint &seed)
1215 {
1216  seed = seed*1664525 + 1013904223;
1217  return SYSwang_inthash(seed);
1218 }
1219 
1220 inline static fpreal32
1221 SYSfastRandomZero(uint &seed)
1222 {
1223  return SYSfastRandom(seed) - 0.5F;
1224 }
1225 
1226 inline static fpreal32
1227 SYSrandomZero(uint &seed)
1228 {
1229  return SYSrandom(seed) - 0.5F;
1230 }
1231 
1232 template <typename T>
1233 inline static void
1234 SYSminmax(T v0, T v1, T v2, T v3, T &min, T &max)
1235 {
1236  min = SYSmin(v0, v1, v2, v3);
1237  max = SYSmax(v0, v1, v2, v3);
1238 }
1239 
1240 template <typename T>
1241 inline static void
1242 SYSminmax(T v0, T v1, T v2, T &min, T &max)
1243 {
1244  min = SYSmin(v0, v1, v2);
1245  max = SYSmax(v0, v1, v2);
1246 }
1247 
1248 template <typename T>
1249 inline static void
1250 SYSminmax(T v0, T v1, T &min, T &max)
1251 {
1252  min = SYSmin(v0, v1);
1253  max = SYSmax(v0, v1);
1254 }
1255 
1256 inline static void
1257 SYSgetSinCosFromSlope(fpreal32 slope, fpreal32 &sintheta, fpreal32 &costheta)
1258 {
1259  fpreal32 one_over_m;
1260  sintheta = slope / SYSsqrt(slope*slope + (fpreal32)1);
1261  if ((slope = SYSabs(slope)) > (fpreal32)1)
1262  {
1263  one_over_m = (fpreal32)1 / slope;
1264  costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
1265  }
1266  else
1267  costheta = SYSsqrt((fpreal32)1 - sintheta*sintheta);
1268 }
1269 
1270 inline static void
1271 SYSgetSinCosFromSlope(fpreal64 slope, fpreal64 &sintheta, fpreal64 &costheta)
1272 {
1273  fpreal64 one_over_m;
1274  sintheta = slope / SYSsqrt(slope*slope + (fpreal64)1);
1275  if ((slope = SYSabs(slope)) > (fpreal64)1)
1276  {
1277  one_over_m = (fpreal64)1 / slope;
1278  costheta = one_over_m / SYSsqrt(one_over_m*one_over_m + 1);
1279  }
1280  else
1281  costheta = SYSsqrt((fpreal64)1 - sintheta*sintheta);
1282 }
1283 
1284 inline static bool
1285 SYSsameSign( fpreal32 v0, fpreal32 v1 )
1286 {
1287  return (v0*v1)>0;
1288 }
1289 
1290 inline static bool
1291 SYSsameSign( fpreal64 v0, fpreal64 v1 )
1292 {
1293  return (v0*v1)>0;
1294 }
1295 
1296 inline static bool
1297 SYSsameSign( int32 v0, int32 v1 )
1298 {
1299  return (v0 ^ v1) >= 0;
1300 }
1301 
1302 inline static bool
1303 SYSsameSign( int64 v0, int64 v1 )
1304 {
1305  return (v0 ^ v1) >= 0;
1306 }
1307 
1308 static inline uint
1309 SYSnextPrime(uint num)
1310 {
1311  return SYSmakePrime(num+1);
1312 }
1313 
1314 static inline int
1315 SYShexCharToInt(char c)
1316 {
1317  // Given a hexidecimal character, return it's corresponding value between
1318  // 0 and 15, or -1 if it's not a hexidecimal character.
1319  if (c >= '0' && c <= '9')
1320  return c - '0';
1321  if (c >= 'a' && c <= 'f')
1322  return c - 'a' + 10;
1323  if (c >= 'A' && c <= 'F')
1324  return c - 'A' + 10;
1325  return -1;
1326 }
1327 
1328 static inline char
1329 SYSintToHexChar(int value)
1330 {
1331  // The value must be from 0-15 for this function to return a valid result.
1332  return value < 10 ? '0' + value : 'a' + value - 10;
1333 }
1334 
1335 SYS_API void SYSsort(int &a, int &b);
1336 SYS_API void SYSsort(int &a, int &b, int &c);
1337 SYS_API void SYSsort(int64 &a, int64 &b);
1338 SYS_API void SYSsort(int64 &a, int64 &b, int64 &c);
1339 SYS_API void SYSsort(float &a, float &b);
1340 SYS_API void SYSsort(float &a, float &b, float &c);
1341 SYS_API void SYSsort(double &a, double &b);
1342 SYS_API void SYSsort(double &a, double &b, double &c);
1343 
1344 // Compute both integer division and integer modulus
1345 // They are compiled into a single instruction
1346 static inline void
1347 SYSdivMod(int numerator, int denominator, int &quotient, int &remainder)
1348 {
1349  quotient = numerator / denominator;
1350  remainder = numerator % denominator;
1351 }
1352 
1353 static inline void
1354 SYSdivMod(int64 numerator, int64 denominator, int64 &quotient, int64 &remainder)
1355 {
1356  quotient = numerator / denominator;
1357  remainder = numerator % denominator;
1358 }
1359 
1360 // Include permutations of fpreal32/fpreal64
1361 #include "SYS_MathPermute.h"
1362 #include "SYS_MathRestrictive.h"
1363 
1364 #else /* cplusplus */
1365 #define SYSmax(a,b) ((a) > (b) ? (a) : (b))
1366 #define SYSmin(a,b) ((a) < (b) ? (a) : (b))
1367 #define SYSabs(a) ((a) < 0 ? (a) : -(a))
1368 #endif
1369 
1370 #endif
SYS_API float acosf(float x)
#define SYSmax(a, b)
Definition: SYS_Math.h:1365
SYS_API double cos(double x)
SYS_API double fmod(double x, double y)
SYS_API double atan2(double y, double x)
SYS_API double expm1(double x)
const hboost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:128
int trunc(T x)
Definition: ImathFun.h:165
const GLdouble * v
Definition: glcorearb.h:836
SYS_API void SYSsrand48(long seed)
png_voidp ptr
Definition: png.h:2145
#define CONST_INT64(x)
Definition: SYS_Types.h:312
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
bool SYSisNan(const F f)
Definition: SYS_Math.h:173
#define SYSabs(a)
Definition: SYS_Math.h:1367
SYS_API float atan2f(float y, float x)
UT_Matrix2T< T > SYSbilerp(const UT_Matrix2T< T > &u0v0, const UT_Matrix2T< T > &u1v0, const UT_Matrix2T< T > &u0v1, const UT_Matrix2T< T > &u1v1, S u, S v)
Definition: UT_Matrix2.h:70
GLint y
Definition: glcorearb.h:102
#define SYS_FTOLERANCE_D
Definition: SYS_Types.h:203
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:817
UT_Matrix2T< T > SYSlerp(const UT_Matrix2T< T > &v1, const UT_Matrix2T< T > &v2, S t)
Definition: UT_Matrix2.h:595
GLfloat GLfloat GLfloat GLfloat v3
Definition: glcorearb.h:818
SYS_API double log10(double x)
SYS_API double pow(double x, double y)
png_uint_32 i
Definition: png.h:2877
#define M_PI
Definition: ImathPlatform.h:51
GLsizeiptr size
Definition: glcorearb.h:663
SYS_API float log10f(float x)
const hboost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
signed char int8
Definition: SYS_Types.h:30
SYS_API double log(double x)
long long int64
Definition: SYS_Types.h:106
GLfloat f
Definition: glcorearb.h:1925
SYS_API double asin(double x)
unsigned long long uint64
Definition: SYS_Types.h:107
int SYSsignum(const F a)
Definition: SYS_Math.h:166
SYS_API double sinh(double x)
SYS_API float atanf(float x)
double fpreal64
Definition: SYS_Types.h:191
SYS_API fpreal32 SYSroundAngle(fpreal32 base, fpreal32 source)
GLsizei GLsizei GLchar * source
Definition: glcorearb.h:802
SYS_API double cosh(double x)
SYS_API bool SYSisPrime(uint num)
SYS_API double exp(double x)
bool isFinite() const
Definition: fpreal16.h:640
SYS_API double SYSdrand48()
SYS_API float logf(float x)
SYS_API double acos(double x)
SYS_API bool SYSisInt(const char *str)
int int32
Definition: SYS_Types.h:34
SYS_API double hypot(double x, double y)
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
fpreal32 SYSrint(fpreal32 val)
Definition: SYS_Floor.h:163
unsigned int uint
Definition: SYS_Types.h:39
SYS_API double log1p(double x)
GLfloat v0
Definition: glcorearb.h:815
SYS_API double tanh(double x)
SYS_API uint SYSmakePrime(uint num)
SYS_API double tan(double x)
unsigned short uint16
Definition: SYS_Types.h:33
GLsizei const GLfloat * value
Definition: glcorearb.h:823
unsigned char uint8
Definition: SYS_Types.h:31
SYS_API double atan(double x)
GLint GLenum GLint x
Definition: glcorearb.h:408
GLfloat GLfloat v1
Definition: glcorearb.h:816
GLuint GLfloat * val
Definition: glcorearb.h:1607
bool isNan() const
Definition: fpreal16.h:672
SYS_API float asinf(float x)
png_infop png_uint_32 int num
Definition: png.h:2158
#define SYS_FTOLERANCE
Definition: SYS_Types.h:202
#define SYS_API
Definition: SYS_API.h:11
#define CONST_UINT64(x)
Definition: SYS_Types.h:313
GLboolean r
Definition: glcorearb.h:1221
unsigned short bits() const
Definition: fpreal16.h:733
short int16
Definition: SYS_Types.h:32
#define SYSmin(a, b)
Definition: SYS_Math.h:1366
SYS_API bool SYSisFloat(const char *str)
float fpreal32
Definition: SYS_Types.h:190
SYS_API double sin(double x)
unsigned int uint32
Definition: SYS_Types.h:35
bool SYSisFinite(fpreal16 f)
Definition: SYS_Math.h:198