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