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