1 #ifndef SimTK_SIMMATRIX_SCALAR_H_ 2 #define SimTK_SIMMATRIX_SCALAR_H_ 270 inline bool signBit(
unsigned char u) {
return false;}
271 inline bool signBit(
unsigned short u) {
return false;}
272 inline bool signBit(
unsigned int u) {
return false;}
273 inline bool signBit(
unsigned long u) {
return false;}
274 inline bool signBit(
unsigned long long u) {
return false;}
284 inline bool signBit(
signed char i) {
return ((
unsigned char)i & 0x80U) != 0;}
285 inline bool signBit(
short i) {
return ((
unsigned short)i & 0x8000U) != 0;}
286 inline bool signBit(
int i) {
return ((
unsigned int)i & 0x80000000U) != 0;}
287 inline bool signBit(
long long i) {
return ((
unsigned long long)i & 0x8000000000000000ULL) != 0;}
289 inline bool signBit(
long i) {
return ((
unsigned long)i
290 & (
unsigned long)LONG_MIN) != 0;}
292 inline bool signBit(
const float& f) {
return std::signbit(f);}
293 inline bool signBit(
const double& d) {
return std::signbit(d);}
311 inline unsigned int sign(
unsigned char u) {
return u==0 ? 0 : 1;}
312 inline unsigned int sign(
unsigned short u) {
return u==0 ? 0 : 1;}
313 inline unsigned int sign(
unsigned int u) {
return u==0 ? 0 : 1;}
314 inline unsigned int sign(
unsigned long u) {
return u==0 ? 0 : 1;}
315 inline unsigned int sign(
unsigned long long u) {
return u==0 ? 0 : 1;}
320 inline int sign(
signed char i) {
return i>0 ? 1 : (i<0 ? -1 : 0);}
321 inline int sign(
short i) {
return i>0 ? 1 : (i<0 ? -1 : 0);}
322 inline int sign(
int i) {
return i>0 ? 1 : (i<0 ? -1 : 0);}
323 inline int sign(
long i) {
return i>0 ? 1 : (i<0 ? -1 : 0);}
324 inline int sign(
long long i) {
return i>0 ? 1 : (i<0 ? -1 : 0);}
326 inline int sign(
const float& x) {
return x>0 ? 1 : (x<0 ? -1 : 0);}
327 inline int sign(
const double& x) {
return x>0 ? 1 : (x<0 ? -1 : 0);}
328 inline int sign(
const long double& x) {
return x>0 ? 1 : (x<0 ? -1 : 0);}
351 inline unsigned char square(
unsigned char u) {
return u*u;}
352 inline unsigned short square(
unsigned short u) {
return u*u;}
353 inline unsigned int square(
unsigned int u) {
return u*u;}
354 inline unsigned long square(
unsigned long u) {
return u*u;}
355 inline unsigned long long square(
unsigned long long u) {
return u*u;}
359 inline signed char square(
signed char i) {
return i*i;}
360 inline short square(
short i) {
return i*i;}
363 inline long long square(
long long i) {
return i*i;}
365 inline float square(
const float& x) {
return x*x;}
366 inline double square(
const double& x) {
return x*x;}
367 inline long double square(
const long double& x) {
return x*x;}
378 template <
class P>
inline 379 std::complex<P>
square(
const std::complex<P>& x) {
380 const P re=x.real(), im=x.imag();
381 return std::complex<P>(re*re-im*im, 2*re*im);
387 template <
class P>
inline 389 const P re=x.real(), nim=x.negImag();
390 return std::complex<P>(re*re-nim*nim, -2*re*nim);
393 template <
class P>
inline 400 template <
class P>
inline 424 inline unsigned char cube(
unsigned char u) {
return u*u*u;}
425 inline unsigned short cube(
unsigned short u) {
return u*u*u;}
426 inline unsigned int cube(
unsigned int u) {
return u*u*u;}
427 inline unsigned long cube(
unsigned long u) {
return u*u*u;}
428 inline unsigned long long cube(
unsigned long long u) {
return u*u*u;}
430 inline char cube(
char c) {
return c*c*c;}
432 inline signed char cube(
signed char i) {
return i*i*i;}
433 inline short cube(
short i) {
return i*i*i;}
434 inline int cube(
int i) {
return i*i*i;}
435 inline long cube(
long i) {
return i*i*i;}
436 inline long long cube(
long long i) {
return i*i*i;}
438 inline float cube(
const float& x) {
return x*x*x;}
439 inline double cube(
const double& x) {
return x*x*x;}
440 inline long double cube(
const long double& x) {
return x*x*x;}
458 template <
class P>
inline 459 std::complex<P>
cube(
const std::complex<P>& x) {
460 const P re=x.real(), im=x.imag(), rr=re*re, ii=im*im;
461 return std::complex<P>(re*(rr-3*ii), im*(3*rr-ii));
468 template <
class P>
inline 471 const P nre=(-x).
real(), nim=(-x).
imag(), rr=nre*nre, ii=nim*nim;
472 return std::complex<P>(nre*(3*ii-rr), nim*(ii-3*rr));
477 template <
class P>
inline 479 const P re=x.real(), nim=x.negImag(), rr=re*re, ii=nim*nim;
480 return std::complex<P>(re*(rr-3*ii), nim*(ii-3*rr));
486 template <
class P>
inline 489 const P nre=(-x).
real(), im=(-x).negImag(), rr=nre*nre, ii=im*im;
490 return std::complex<P>(nre*(3*ii-rr), im*(3*rr-ii));
542 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
545 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
547 inline long double&
clampInPlace(
long double low,
long double& v,
long double high)
548 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
564 {
return clampInPlace((
long double)low,v,(
long double)high); }
576 inline long double&
clampInPlace(
int low,
long double& v,
long double high)
589 inline long double&
clampInPlace(
long double low,
long double& v,
int high)
593 inline unsigned char&
clampInPlace(
unsigned char low,
unsigned char& v,
unsigned char high)
594 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
596 inline unsigned short&
clampInPlace(
unsigned short low,
unsigned short& v,
unsigned short high)
597 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
599 inline unsigned int&
clampInPlace(
unsigned int low,
unsigned int& v,
unsigned int high)
600 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
602 inline unsigned long&
clampInPlace(
unsigned long low,
unsigned long& v,
unsigned long high)
603 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
605 inline unsigned long long&
clampInPlace(
unsigned long long low,
unsigned long long& v,
unsigned long long high)
606 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
610 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
612 inline signed char&
clampInPlace(
signed char low,
signed char& v,
signed char high)
613 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
617 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
620 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
623 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
625 inline long long&
clampInPlace(
long long low,
long long& v,
long long high)
626 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
630 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
633 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
636 { assert(low<=high);
if (v<low) v=low;
else if (v>high) v=high;
return v; }
660 inline double clamp(
double low,
double v,
double high)
663 inline float clamp(
float low,
float v,
float high)
666 inline long double clamp(
long double low,
long double v,
long double high)
671 inline double clamp(
int low,
double v,
int high)
675 inline float clamp(
int low,
float v,
int high)
679 inline long double clamp(
int low,
long double v,
int high)
684 inline double clamp(
int low,
double v,
double high)
688 inline float clamp(
int low,
float v,
float high)
692 inline long double clamp(
int low,
long double v,
long double high)
697 inline double clamp(
double low,
double v,
int high)
701 inline float clamp(
float low,
float v,
int high)
705 inline long double clamp(
long double low,
long double v,
int high)
709 inline unsigned char clamp(
unsigned char low,
unsigned char v,
unsigned char high)
712 inline unsigned short clamp(
unsigned short low,
unsigned short v,
unsigned short high)
715 inline unsigned int clamp(
unsigned int low,
unsigned int v,
unsigned int high)
718 inline unsigned long clamp(
unsigned long low,
unsigned long v,
unsigned long high)
721 inline unsigned long long clamp(
unsigned long long low,
unsigned long long v,
unsigned long long high)
725 inline char clamp(
char low,
char v,
char high)
728 inline signed char clamp(
signed char low,
signed char v,
signed char high)
732 inline short clamp(
short low,
short v,
short high)
735 inline int clamp(
int low,
int v,
int high)
738 inline long clamp(
long low,
long v,
long high)
741 inline long long clamp(
long long low,
long long v,
long long high)
756 {
return clamp(low,(
float)v,high); }
762 {
return clamp(low,(
double)v,high); }
768 {
return clamp(low,(
long double)v,high); }
825 { assert(0 <= x && x <= 1);
826 return x*x*x*(10+x*(6*x-15)); }
920 inline double stepAny(
double y0,
double yRange,
921 double x0,
double oneOverXRange,
923 {
double xadj = (x-x0)*oneOverXRange;
927 return y0 + yRange*
stepUp(xadj); }
933 assert(0 <= x && x <= 1);
934 const double xxm1=x*(x-1);
947 double x0,
double oneOverXRange,
949 {
double xadj = (x-x0)*oneOverXRange;
953 return yRange*oneOverXRange*
dstepUp(xadj); }
959 assert(0 <= x && x <= 1);
960 return 60*x*(1+x*(2*x-3));
972 double x0,
double oneOverXRange,
974 {
double xadj = (x-x0)*oneOverXRange;
984 assert(0 <= x && x <= 1);
985 return 60+360*x*(x-1);
997 double x0,
double oneOverXRange,
999 {
double xadj = (x-x0)*oneOverXRange;
1009 { assert(0 <= x && x <= 1);
1010 return x*x*x*(10+x*(6*x-15)); }
1015 float x0,
float oneOverXRange,
1017 {
float xadj = (x-x0)*oneOverXRange;
1021 return y0 + yRange*
stepUp(xadj); }
1025 { assert(0 <= x && x <= 1);
1026 const float xxm1=x*(x-1);
1027 return 30*xxm1*xxm1; }
1032 float x0,
float oneOverXRange,
1034 {
float xadj = (x-x0)*oneOverXRange;
1038 return yRange*oneOverXRange*
dstepUp(xadj); }
1042 { assert(0 <= x && x <= 1);
1043 return 60*x*(1+x*(2*x-3)); }
1048 float x0,
float oneOverXRange,
1050 {
float xadj = (x-x0)*oneOverXRange;
1058 { assert(0 <= x && x <= 1);
1059 return 60+360*x*(x-1); }
1064 float x0,
float oneOverXRange,
1066 {
float xadj = (x-x0)*oneOverXRange;
1076 { assert(0 <= x && x <= 1);
1077 return x*x*x*(10+x*(6*x-15)); }
1081 inline long double stepAny(
long double y0,
long double yRange,
1082 long double x0,
long double oneOverXRange,
1084 {
long double xadj = (x-x0)*oneOverXRange;
1088 return y0 + yRange*
stepUp(xadj); }
1093 { assert(0 <= x && x <= 1);
1094 const long double xxm1=x*(x-1);
1095 return 30*xxm1*xxm1; }
1100 long double x0,
long double oneOverXRange,
1102 {
long double xadj = (x-x0)*oneOverXRange;
1106 return yRange*oneOverXRange*
dstepUp(xadj); }
1110 { assert(0 <= x && x <= 1);
1111 return 60*x*(1+x*(2*x-3)); }
1116 long double x0,
long double oneOverXRange,
1118 {
long double xadj = (x-x0)*oneOverXRange;
1127 { assert(0 <= x && x <= 1);
1128 return 60+360*x*(x-1); }
1133 long double x0,
long double oneOverXRange,
1135 {
long double xadj = (x-x0)*oneOverXRange;
1211 static inline std::pair<T,T> approxCompleteEllipticIntegralsKE_T(T m) {
1212 static const T a[] =
1213 { (T)1.38629436112, (T)0.09666344259, (T)0.03590092383,
1214 (T)0.03742563713, (T)0.01451196212 };
1215 static const T b[] =
1216 { (T)0.5, (T)0.12498593597, (T)0.06880248576,
1217 (T)0.03328355346, (T)0.00441787012 };
1218 static const T c[] =
1219 { (T)1, (T)0.44325141463, (T)0.06260601220,
1220 (T)0.04757383546, (T)0.01736506451 };
1221 static const T d[] =
1222 { (T)0, (T)0.24998368310, (T)0.09200180037,
1223 (T)0.04069697526, (T)0.00526449639 };
1230 "approxCompleteEllipticIntegralsKE()",
1231 "Argument m (%g) is outside the legal range [0,1].", (
double)m);
1232 if (m >= 1)
return std::make_pair(Infinity, (T)1);
1233 if (m <= 0)
return std::make_pair(PiOver2, PiOver2);
1235 const T m1=1-m, m2=m1*m1, m3=m1*m2, m4=m2*m2;
1236 const T lnm2 = std::log(m1);
1239 const T K = (a[0] + a[1]*m1 + a[2]*m2 + a[3]*m3 + a[4]*m4)
1240 - lnm2*(b[0] + b[1]*m1 + b[2]*m2 + b[3]*m3 + b[4]*m4);
1241 const T E = (c[0] + c[1]*m1 + c[2]*m2 + c[3]*m3 + c[4]*m4)
1242 - lnm2*( d[1]*m1 + d[2]*m2 + d[3]*m3 + d[4]*m4);
1244 return std::make_pair(K,E);
1284 inline std::pair<double,double>
1286 {
return approxCompleteEllipticIntegralsKE_T<double>(m); }
1292 inline std::pair<float,float>
1294 {
return approxCompleteEllipticIntegralsKE_T<float>(m); }
1301 inline std::pair<double,double>
1303 {
return approxCompleteEllipticIntegralsKE_T<double>((double)m); }
1308 static inline std::pair<T,T> completeEllipticIntegralsKE_T(T m) {
1318 "completeEllipticIntegralsKE()",
1319 "Argument m (%g) is outside the legal range [0,1].", (
double)m);
1320 if (m >= 1)
return std::make_pair(Infinity, (T)1);
1321 if (m <= 0)
return std::make_pair(PiOver2, PiOver2);
1323 const T k = std::sqrt(1-m);
1324 T v1=1, w1=k, c1=1, d1=k*k;
1327 T w2 = std::sqrt(v1*w1);
1329 T d2 = (w1*c1+v1*d1)/(v1+w1);
1330 v1=v2; w1=w2; c1=c2; d1=d2;
1331 }
while(
std::abs(v1-w1) >= TenEps);
1333 const T K = PiOver2/v1;
1336 return std::make_pair(K,E);
1365 {
return completeEllipticIntegralsKE_T<double>(m); }
1374 {
return completeEllipticIntegralsKE_T<float>(m); }
1382 {
return completeEllipticIntegralsKE_T<double>((double)m); }
1388 #endif //SimTK_SIMMATRIX_SCALAR_H_ const Real MostNegativeReal
This is the smallest finite negative real number that can be expressed in values of type Real...
const Real Ln2
Real(ln(2)) (natural log of 2)
#define SimTK_SimTKCOMMON_EXPORT
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:218
const Real Sqrt3
Real(sqrt(3))
double d2stepDown(double x)
Second derivative of stepDown(): d^2/dx^2 stepDown(x).
Definition: Scalar.h:966
const Real LeastNegativeReal
This is the largest negative real number (that is, closest to zero) that can be expressed in values o...
double d2stepAny(double yRange, double x0, double oneOverXRange, double x)
Second derivative of stepAny(): d^2/dx^2 stepAny(x).
Definition: Scalar.h:971
bool signBit(unsigned char u)
Definition: Scalar.h:270
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
SimTK::conjugate<R> should be instantiated only for float, double, long double.
Definition: String.h:45
This file defines the negator<N> template which is an adaptor for the numeric types N (Real...
This file defines the conjugate<R> template class, where R is one of the three built-in real types...
double & clampInPlace(double low, double &v, double high)
Check that low <= v <= high and modify v in place if necessary to bring it into that range...
Definition: Scalar.h:541
std::complex< Real > Complex
This is the default complex type for SimTK, with precision for the real and imaginary parts set to th...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:596
conjugate< Real > Conjugate
Definition: Scalar.h:57
const Real OneFifth
Real(1)/5.
const Real Sqrt2
Real(sqrt(2))
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
SimTK_Real Real
This is the default compiled-in floating point type for SimTK, either float or double.
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:593
bool exactlyOneBitIsSet(unsigned char v)
Definition: Scalar.h:230
std::pair< double, double > completeEllipticIntegralsKE(double m)
Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m)...
Definition: Scalar.h:1364
This file contains classes and typedefs needed to provide uniform handling of floating point numeric ...
const Real LeastPositiveReal
This is the smallest positive real number that can be expressed in the type Real; it is ~1e-308 when ...
const Real NaN
This is the IEEE "not a number" constant for this implementation of the default-precision Real type; ...
unsigned char square(unsigned char u)
Definition: Scalar.h:351
This file defines the Array_<T,X> class and related support classes including base classes ArrayViewC...
This file contains macros which are convenient to use for sprinkling error checking around liberally ...
const Real CubeRoot2
Real(2^(1/3)) (cube root of 2)
bool atMostOneBitIsSet(unsigned char v)
Definition: Scalar.h:203
const Complex I
We only need one complex constant, i = sqrt(-1). For the rest just multiply the real constant by i...
const Real OneSixth
Real(1)/6.
#define SimTK_ERRCHK1_ALWAYS(cond, whereChecked, fmt, a1)
Definition: ExceptionMacros.h:285
double dstepAny(double yRange, double x0, double oneOverXRange, double x)
First derivative of stepAny(): d/dx stepAny(x).
Definition: Scalar.h:946
High precision mathematical and physical constants.
The purpose of the CNT<T> class is to hide the differences between built-in numerical types and compo...
const Real OneEighth
Real(1)/8.
double d2stepUp(double x)
Second derivative of stepUp(): d^2/dx^2 stepUp(x).
Definition: Scalar.h:958
double dstepDown(double x)
First derivative of stepDown(): d/dx stepDown(x).
Definition: Scalar.h:941
const Real SqrtEps
This is the square root of Eps, ~1e-8 if Real is double, ~3e-4 if Real is float.
const Real Log10E
Real(log10(e)) (log base 10)
const int NumDigitsReal
This is the number of decimal digits that can be reliably stored and retrieved in the default Real pr...
const Real TinyReal
TinyReal is a floating point value smaller than the floating point precision; it is defined as Eps^(5...
const Real OneFourth
Real(1)/4.
double stepUp(double x)
Interpolate smoothly from 0 up to 1 as the input argument goes from 0 to 1, with first and second der...
Definition: Scalar.h:824
const float & real(const conjugate< float > &c)
Definition: conjugate.h:771
const Real MinusOne
Real(-1)
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
const Real Infinity
This is the IEEE positive infinity constant for this implementation of the default-precision Real typ...
double d3stepAny(double yRange, double x0, double oneOverXRange, double x)
Third derivative of stepAny(): d^3/dx^3 stepAny(x).
Definition: Scalar.h:996
const Real Eps
Epsilon is the size of roundoff noise; it is the smallest positive number of default-precision type R...
double clamp(double low, double v, double high)
If v is in range low <= v <= high then return v, otherwise return the nearest bound; this function do...
Definition: Scalar.h:660
const Real SignificantReal
SignificantReal is the smallest value that we consider to be clearly distinct from roundoff error whe...
const Real E
e = Real(exp(1))
Mandatory first inclusion for any Simbody source or header file.
static const negator< N > & recast(const N &val)
Definition: negator.h:236
const Real Ln10
Real(ln(10)) (natural log of 10)
const Real Log2E
Real(log2(e)) (log base 2)
const Real OneHalf
Real(1)/2.
const Real OneSeventh
Real(1)/7.
unsigned int sign(unsigned char u)
Definition: Scalar.h:311
double stepAny(double y0, double yRange, double x0, double oneOverXRange, double x)
Interpolate smoothly from y0 to y1 as the input argument goes from x0 to x1, with first and second de...
Definition: Scalar.h:920
const Real OneOverPi
1/Real(pi)
const Real OneOverSqrt2
1/sqrt(2)==sqrt(2)/2 as Real
std::pair< double, double > approxCompleteEllipticIntegralsKE(double m)
Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m)...
Definition: Scalar.h:1285
const Real MostPositiveReal
This is the largest finite positive real number that can be expressed in the Real type; ~1e+308 when ...
const negator< float > & imag(const conjugate< float > &c)
Definition: conjugate.h:772
double d3stepUp(double x)
Third derivative of stepUp(): d^3/dx^3 stepUp(x).
Definition: Scalar.h:983
const int LosslessNumDigitsReal
This is the smallest number of decimal digits you should store in a text file if you want to be able ...
const Real OneOverSqrt3
Real(1/sqrt(3))
const Real OneNinth
Real(1)/9.
double stepDown(double x)
Interpolate smoothly from 1 down to 0 as the input argument goes from 0 to 1, with first and second d...
Definition: Scalar.h:843
const Real OneThird
Real(1)/3.
double dstepUp(double x)
First derivative of stepUp(): d/dx stepUp(x).
Definition: Scalar.h:932
const Real CubeRoot3
Real(3^(1/3)) (cube root of 3)
double d3stepDown(double x)
Third derivative of stepDown(): d^3/dx^3 stepDown(x).
Definition: Scalar.h:991
unsigned char cube(unsigned char u)
Definition: Scalar.h:424