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