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