|
Simbody
|
NTraits.h
Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_NTRAITS_H_ 00002 #define SimTK_SIMMATRIX_NTRAITS_H_ 00003 00004 /* -------------------------------------------------------------------------- * 00005 * SimTK Core: SimTK Simmatrix(tm) * 00006 * -------------------------------------------------------------------------- * 00007 * This is part of the SimTK Core biosimulation toolkit originating from * 00008 * Simbios, the NIH National Center for Physics-Based Simulation of * 00009 * Biological Structures at Stanford, funded under the NIH Roadmap for * 00010 * Medical Research, grant U54 GM072970. See https://simtk.org. * 00011 * * 00012 * Portions copyright (c) 2005-9 Stanford University and the Authors. * 00013 * Authors: Michael Sherman * 00014 * Contributors: * 00015 * * 00016 * Permission is hereby granted, free of charge, to any person obtaining a * 00017 * copy of this software and associated documentation files (the "Software"), * 00018 * to deal in the Software without restriction, including without limitation * 00019 * the rights to use, copy, modify, merge, publish, distribute, sublicense, * 00020 * and/or sell copies of the Software, and to permit persons to whom the * 00021 * Software is furnished to do so, subject to the following conditions: * 00022 * * 00023 * The above copyright notice and this permission notice shall be included in * 00024 * all copies or substantial portions of the Software. * 00025 * * 00026 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * 00027 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * 00028 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * 00029 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, * 00030 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR * 00031 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE * 00032 * USE OR OTHER DEALINGS IN THE SOFTWARE. * 00033 * -------------------------------------------------------------------------- */ 00034 00058 #include <cstddef> 00059 #include <cassert> 00060 #include <complex> 00061 #include <iostream> 00062 #include <limits> 00063 00064 using std::complex; 00065 00066 namespace SimTK { 00067 00068 // This is the 3rd type of number, conjugate. It is like complex except that 00069 // the represented value is the conjugate of the value represented by a complex 00070 // number containing the same bit pattern. That is, complex and conjugate both 00071 // contain two real numbers, re and im, with complex(re,im) meaning 00072 // (re + im*i) while conjugate(re,im) means (re - im*i). It is guaranteed that 00073 // our conjugate type and complex type have identical sizes and representations. 00074 // Together, these defininitions and guarantees permit conjugation 00075 // to be done by reinterpretation rather than be computation. 00076 template <class R> class conjugate; // Only defined for float, double, long double 00077 00078 // Specializations of this class provide information about Composite Numerical 00079 // Types in the style of std::numeric_limits<T>. It is specialized for the 00080 // numeric types but can be invoked on any composite numerical type as well. 00081 template <class T> class CNT; 00082 00083 // NTraits provides information like CNT, but for numeric types only. 00084 template <class N> class NTraits; // allowed only for N=<number> 00085 template <class R> class NTraits< complex<R> >; 00086 template <class R> class NTraits< conjugate<R> >; 00087 template <> class NTraits<float>; 00088 template <> class NTraits<double>; 00089 template <> class NTraits<long double>; 00090 00091 // This is an adaptor for numeric types which negates the apparent values. A 00092 // negator<N> has exactly the same internal representation as a numeric value N, 00093 // but it is to be interpreted has having the negative of the value it would 00094 // have if interpreted as an N. This permits negation to be done by 00095 // reinterpretation rather than compuation. A full set of arithmetic operators 00096 // are provided involving negator<N>'s and N's. Sometimes we can save an op or 00097 // two this way. For example negator<N>*negator<N> can be performed as an N*N 00098 // since the negations cancel, and we saved two floating point negations. 00099 template <class N> class negator; // Only defined for numbers 00100 00101 // This is here so we can provide references to 0's when needed, e.g. 00102 // when returning the imaginary part of a real number. These are local to 00103 // the compilation unit, so the returned addresses will differ in different 00104 // files. There are enough zeroes here for the widest number, 00105 // complex<long double> (or conjugate<long double>). 00106 static const complex<long double> zeroes(0); 00107 00115 template <class R1, class R2> struct Widest {/* Only defined for built-ins. */}; 00116 template <> struct Widest<float,float> {typedef float Type; typedef float Precision;}; 00117 template <> struct Widest<float,double> {typedef double Type; typedef double Precision;}; 00118 template <> struct Widest<float,long double> {typedef long double Type; typedef long double Precision;}; 00119 template <> struct Widest<double,float> {typedef double Type; typedef double Precision;}; 00120 template <> struct Widest<double,double> {typedef double Type; typedef double Precision;}; 00121 template <> struct Widest<double,long double> {typedef long double Type; typedef long double Precision;}; 00122 template <> struct Widest<long double,float> {typedef long double Type; typedef long double Precision;}; 00123 template <> struct Widest<long double,double> {typedef long double Type; typedef long double Precision;}; 00124 template <> struct Widest<long double,long double> {typedef long double Type; typedef long double Precision;}; 00125 template <class R1, class R2> struct Widest< complex<R1>,complex<R2> > { 00126 typedef complex< typename Widest<R1,R2>::Type > Type; 00127 typedef typename Widest<R1,R2>::Precision Precision; 00128 }; 00129 template <class R1, class R2> struct Widest< complex<R1>,R2 > { 00130 typedef complex< typename Widest<R1,R2>::Type > Type; 00131 typedef typename Widest<R1,R2>::Precision Precision; 00132 }; 00133 template <class R1, class R2> struct Widest< R1,complex<R2> > { 00134 typedef complex< typename Widest<R1,R2>::Type > Type; 00135 typedef typename Widest<R1,R2>::Precision Precision; 00136 }; 00137 00147 template <class R1, class R2> struct Narrowest {/* Only defined for built-ins. */}; 00148 template <> struct Narrowest<float,float> {typedef float Type;}; 00149 template <> struct Narrowest<float,double> {typedef float Type;}; 00150 template <> struct Narrowest<float,long double> {typedef float Type;}; 00151 template <> struct Narrowest<double,float> {typedef float Type;}; 00152 template <> struct Narrowest<double,double> {typedef double Type;}; 00153 template <> struct Narrowest<double,long double> {typedef double Type;}; 00154 template <> struct Narrowest<long double,float> {typedef float Type;}; 00155 template <> struct Narrowest<long double,double> {typedef double Type;}; 00156 template <> struct Narrowest<long double,long double> {typedef long double Type;}; 00157 template <class R1, class R2> struct Narrowest< complex<R1>,complex<R2> > { 00158 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00159 typedef typename Narrowest<R1,R2>::Precision Precision; 00160 }; 00161 template <class R1, class R2> struct Narrowest< complex<R1>,R2 > { 00162 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00163 typedef typename Narrowest<R1,R2>::Precision Precision; 00164 }; 00165 template <class R1, class R2> struct Narrowest< R1,complex<R2> > { 00166 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00167 typedef typename Narrowest<R1,R2>::Precision Precision; 00168 }; 00169 00171 template <class R> class RTraits {/* Only defined for real types */}; 00172 template <> class RTraits<float> { 00173 public: 00175 static const float& getEps() {static const float c=std::numeric_limits<float>::epsilon(); return c;} 00177 static const float& getSignificant() {static const float c=std::pow(getEps(), 0.875f); return c;} 00179 static double getDefaultTolerance() {return (double)getSignificant();} 00180 }; 00181 template <> class RTraits<double> { 00182 public: 00183 static const double& getEps() {static const double c=std::numeric_limits<double>::epsilon(); return c;} 00184 static const double& getSignificant() {static const double c=std::pow(getEps(), 0.875); return c;} 00185 static double getDefaultTolerance() {return getSignificant();} 00186 }; 00187 template <> class RTraits<long double> { 00188 public: 00189 static const long double& getEps() {static const long double c=std::numeric_limits<long double>::epsilon(); return c;} 00190 static const long double& getSignificant() {static const long double c=std::pow(getEps(), 0.875L); return c;} 00191 static double getDefaultTolerance() {return (double)getSignificant();} 00192 }; 00193 00210 // See negator.h for isNaN() applied to negated scalars. 00212 inline bool isNaN(const float& x) {return std::isnan(x);} 00213 inline bool isNaN(const double& x) {return std::isnan(x);} 00214 inline bool isNaN(const long double& x) {return std::isnan(x);} 00215 00216 template <class P> inline bool 00217 isNaN(const std::complex<P>& x) 00218 { return isNaN(x.real()) || isNaN(x.imag());} 00219 00220 template <class P> inline bool 00221 isNaN(const conjugate<P>& x) 00222 { return isNaN(x.real()) || isNaN(x.negImag());} 00224 00237 // See negator.h for isFinite() applied to negated scalars. 00239 inline bool isFinite(const float& x) {return std::isfinite(x);} 00240 inline bool isFinite(const double& x) {return std::isfinite(x);} 00241 inline bool isFinite(const long double& x) {return std::isfinite(x);} 00242 00243 template <class P> inline bool 00244 isFinite(const std::complex<P>& x) 00245 { return isFinite(x.real()) && isFinite(x.imag());} 00246 00247 template <class P> inline bool 00248 isFinite(const conjugate<P>& x) 00249 { return isFinite(x.real()) && isFinite(x.negImag());} 00251 00266 // See negator.h for isInf() applied to negated scalars. 00268 inline bool isInf(const float& x) {return std::isinf(x);} 00269 inline bool isInf(const double& x) {return std::isinf(x);} 00270 inline bool isInf(const long double& x) {return std::isinf(x);} 00271 00272 template <class P> inline bool 00273 isInf(const std::complex<P>& x) { 00274 return (isInf(x.real()) && !isNaN(x.imag())) 00275 || (isInf(x.imag()) && !isNaN(x.real())); 00276 } 00277 00278 template <class P> inline bool 00279 isInf(const conjugate<P>& x) { 00280 return (isInf(x.real()) && !isNaN(x.negImag())) 00281 || (isInf(x.negImag()) && !isNaN(x.real())); 00282 } 00284 00324 00325 inline bool isNumericallyEqual(const float& a, const float& b, 00326 double tol = RTraits<float>::getDefaultTolerance()) 00327 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00328 const float scale = std::max(std::max(std::abs(a),std::abs(b)), 1.f); 00329 return std::abs(a-b) <= scale*(float)tol; } 00331 inline bool isNumericallyEqual(const double& a, const double& b, 00332 double tol = RTraits<double>::getDefaultTolerance()) 00333 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00334 const double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.); 00335 return std::abs(a-b) <= scale*tol; } 00337 inline bool isNumericallyEqual(const long double& a, const long double& b, 00338 double tol = RTraits<long double>::getDefaultTolerance()) 00339 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00340 const long double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.L); 00341 return std::abs(a-b) <= scale*(long double)tol; } 00342 00344 inline bool isNumericallyEqual(const float& a, const double& b, 00345 double tol = RTraits<float>::getDefaultTolerance()) 00346 { return isNumericallyEqual((double)a, b, tol); } 00348 inline bool isNumericallyEqual(const double& a, const float& b, 00349 double tol = RTraits<float>::getDefaultTolerance()) 00350 { return isNumericallyEqual(a, (double)b, tol); } 00352 inline bool isNumericallyEqual(const float& a, const long double& b, 00353 double tol = RTraits<float>::getDefaultTolerance()) 00354 { return isNumericallyEqual((long double)a, b, tol); } 00356 inline bool isNumericallyEqual(const long double& a, const float& b, 00357 double tol = RTraits<float>::getDefaultTolerance()) 00358 { return isNumericallyEqual(a, (long double)b, tol); } 00360 inline bool isNumericallyEqual(const double& a, const long double& b, 00361 double tol = RTraits<double>::getDefaultTolerance()) 00362 { return isNumericallyEqual((long double)a, b, tol); } 00364 inline bool isNumericallyEqual(const long double& a, const double& b, 00365 double tol = RTraits<double>::getDefaultTolerance()) 00366 { return isNumericallyEqual(a, (long double)b, tol); } 00367 00369 inline bool isNumericallyEqual(const float& a, int b, 00370 double tol = RTraits<float>::getDefaultTolerance()) 00371 { return isNumericallyEqual(a, (double)b, tol); } 00373 inline bool isNumericallyEqual(int a, const float& b, 00374 double tol = RTraits<float>::getDefaultTolerance()) 00375 { return isNumericallyEqual((double)a, b, tol); } 00377 inline bool isNumericallyEqual(const double& a, int b, 00378 double tol = RTraits<double>::getDefaultTolerance()) 00379 { return isNumericallyEqual(a, (double)b, tol); } 00381 inline bool isNumericallyEqual(int a, const double& b, 00382 double tol = RTraits<double>::getDefaultTolerance()) 00383 { return isNumericallyEqual((double)a, b, tol); } 00385 inline bool isNumericallyEqual(const long double& a, int b, 00386 double tol = RTraits<long double>::getDefaultTolerance()) 00387 { return isNumericallyEqual(a, (long double)b, tol); } 00389 inline bool isNumericallyEqual(int a, const long double& b, 00390 double tol = RTraits<long double>::getDefaultTolerance()) 00391 { return isNumericallyEqual((long double)a, b, tol); } 00392 00396 template <class P, class Q> 00397 inline bool isNumericallyEqual 00398 ( const std::complex<P>& a, const std::complex<Q>& b, 00399 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00400 { return isNumericallyEqual(a.real(),b.real(),tol) 00401 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00402 00406 template <class P, class Q> 00407 inline bool isNumericallyEqual 00408 ( const conjugate<P>& a, const conjugate<Q>& b, 00409 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00410 { return isNumericallyEqual(a.real(),b.real(),tol) 00411 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00412 00416 template <class P, class Q> 00417 inline bool isNumericallyEqual 00418 ( const std::complex<P>& a, const conjugate<Q>& b, 00419 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00420 { return isNumericallyEqual(a.real(),b.real(),tol) 00421 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00422 00426 template <class P, class Q> 00427 inline bool isNumericallyEqual 00428 ( const conjugate<P>& a, const std::complex<Q>& b, 00429 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00430 { return isNumericallyEqual(a.real(),b.real(),tol) 00431 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00432 00434 template <class P> inline bool 00435 isNumericallyEqual(const std::complex<P>& a, const float& b, 00436 double tol = RTraits<float>::getDefaultTolerance()) 00437 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); } 00439 template <class P> inline bool 00440 isNumericallyEqual(const float& a, const std::complex<P>& b, 00441 double tol = RTraits<float>::getDefaultTolerance()) 00442 { return isNumericallyEqual(b,a,tol); } 00444 template <class P> inline bool 00445 isNumericallyEqual(const std::complex<P>& a, const double& b, 00446 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00447 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); } 00449 template <class P> inline bool 00450 isNumericallyEqual(const double& a, const std::complex<P>& b, 00451 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00452 { return isNumericallyEqual(b,a,tol); } 00454 template <class P> inline bool 00455 isNumericallyEqual(const std::complex<P>& a, const long double& b, 00456 double tol = RTraits<P>::getDefaultTolerance()) 00457 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); } 00459 template <class P> inline bool 00460 isNumericallyEqual(const long double& a, const std::complex<P>& b, 00461 double tol = RTraits<P>::getDefaultTolerance()) 00462 { return isNumericallyEqual(b,a,tol); } 00464 template <class P> inline bool 00465 isNumericallyEqual(const std::complex<P>& a, int b, 00466 double tol = RTraits<P>::getDefaultTolerance()) 00467 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); } 00469 template <class P> inline bool 00470 isNumericallyEqual(int a, const std::complex<P>& b, 00471 double tol = RTraits<P>::getDefaultTolerance()) 00472 { return isNumericallyEqual(b,a,tol); } 00473 00475 template <class P> inline bool 00476 isNumericallyEqual(const conjugate<P>& a, const float& b, 00477 double tol = RTraits<float>::getDefaultTolerance()) 00478 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); } 00480 template <class P> inline bool 00481 isNumericallyEqual(const float& a, const conjugate<P>& b, 00482 double tol = RTraits<float>::getDefaultTolerance()) 00483 { return isNumericallyEqual(b,a,tol); } 00485 template <class P> inline bool 00486 isNumericallyEqual(const conjugate<P>& a, const double& b, 00487 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00488 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); } 00490 template <class P> inline bool 00491 isNumericallyEqual(const double& a, const conjugate<P>& b, 00492 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00493 { return isNumericallyEqual(b,a,tol); } 00495 template <class P> inline bool 00496 isNumericallyEqual(const conjugate<P>& a, const long double& b, 00497 double tol = RTraits<P>::getDefaultTolerance()) 00498 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); } 00500 template <class P> inline bool 00501 isNumericallyEqual(const long double& a, const conjugate<P>& b, 00502 double tol = RTraits<P>::getDefaultTolerance()) 00503 { return isNumericallyEqual(b,a,tol); } 00505 template <class P> inline bool 00506 isNumericallyEqual(const conjugate<P>& a, int b, 00507 double tol = RTraits<P>::getDefaultTolerance()) 00508 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); } 00510 template <class P> inline bool 00511 isNumericallyEqual(int a, const conjugate<P>& b, 00512 double tol = RTraits<P>::getDefaultTolerance()) 00513 { return isNumericallyEqual(b,a,tol); } 00514 00516 00517 00518 template <class N> class NTraits { 00519 // only the specializations below are allowed 00520 }; 00521 00524 template <class R> class NTraits< complex<R> > { 00525 typedef complex<R> C; 00526 public: 00527 typedef C T; 00528 typedef negator<C> TNeg; // type of this after *cast* to its negative 00529 typedef C TWithoutNegator; // type of this ignoring negator (there isn't one!) 00530 00531 typedef R TReal; 00532 typedef R TImag; 00533 typedef C TComplex; 00534 typedef conjugate<R> THerm; // this is a recast 00535 typedef C TPosTrans; 00536 typedef R TSqHermT; // ~C*C 00537 typedef R TSqTHerm; // C*~C (same) 00538 typedef C TElement; 00539 typedef C TRow; 00540 typedef C TCol; 00541 00542 typedef C TSqrt; 00543 typedef R TAbs; 00544 typedef C TStandard; // complex is a standard type 00545 typedef C TInvert; // this is a calculation, so use standard number 00546 typedef C TNormalize; 00547 00548 typedef C Scalar; 00549 typedef C ULessScalar; 00550 typedef C Number; 00551 typedef C StdNumber; 00552 typedef R Precision; 00553 typedef R ScalarNormSq; 00554 00555 // For complex scalar C, op result types are: 00556 // Typeof(C*P) = Typeof(P*C) 00557 // Typeof(C/P) = Typeof(inv(P)*C) 00558 // Typeof(C+P) = Typeof(P+C) 00559 // typeof(C-P) = Typeof(P::TNeg + C) 00560 // These must be specialized for P=real, complex, conjugate. 00561 template <class P> struct Result { 00562 typedef typename CNT<P>::template Result<C>::Mul Mul; 00563 typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd; 00564 typedef typename CNT<P>::template Result<C>::Add Add; 00565 typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub; 00566 }; 00567 00568 // Shape-preserving element substitution (easy for scalars!) 00569 template <class P> struct Substitute { 00570 typedef P Type; 00571 }; 00572 00573 enum { 00574 NRows = 1, 00575 NCols = 1, 00576 RowSpacing = 1, 00577 ColSpacing = 1, 00578 NPackedElements = 1, // not two! 00579 NActualElements = 1, 00580 NActualScalars = 1, 00581 ImagOffset = 1, 00582 RealStrideFactor = 2, // double stride when casting to real or imaginary 00583 ArgDepth = SCALAR_DEPTH, 00584 IsScalar = 1, 00585 IsULessScalar = 1, 00586 IsNumber = 1, 00587 IsStdNumber = 1, 00588 IsPrecision = 0, 00589 SignInterpretation = 1 // after cast to Number, what is the sign? 00590 }; 00591 static const T* getData(const T& t) { return &t; } 00592 static T* updData(T& t) { return &t; } 00593 static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; } 00594 static R& real(T& t) { return (reinterpret_cast<R*>(&t))[0]; } 00595 static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; } 00596 static R& imag(T& t) { return (reinterpret_cast<R*>(&t))[1]; } 00597 00598 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} 00599 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} 00600 00601 static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} 00602 static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} 00603 00604 static const TPosTrans& positionalTranspose(const T& t) 00605 {return reinterpret_cast<const TPosTrans&>(t);} 00606 static TPosTrans& positionalTranspose(T& t) 00607 {return reinterpret_cast<TPosTrans&>(t);} 00608 00609 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) 00610 {return reinterpret_cast<const TWithoutNegator&>(t);} 00611 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) 00612 {return reinterpret_cast<TWithoutNegator&>(t);} 00613 00614 static ScalarNormSq scalarNormSqr(const T& t) 00615 { return t.real()*t.real() + t.imag()*t.imag(); } 00616 static TSqrt sqrt(const T& t) 00617 { return std::sqrt(t); } 00618 static TAbs abs(const T& t) 00619 { return std::abs(t); } // no, not just sqrt of scalarNormSqr()! 00620 static const TStandard& standardize(const T& t) {return t;} // already standard 00621 static TNormalize normalize(const T& t) {return t/abs(t);} 00622 static TInvert invert(const T& t) {return TReal(1)/t;} 00623 00624 static const T& getNaN() { 00625 static const T c=T(NTraits<R>::getNaN(), NTraits<R>::getNaN()); 00626 return c; 00627 } 00628 static const T& getInfinity() { 00629 static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity()); 00630 return c; 00631 } 00632 00633 static const T& getI() { 00634 static const T c = T(0,1); 00635 return c; 00636 } 00637 00638 static bool isFinite(const T& t) {return SimTK::isFinite(t);} 00639 static bool isNaN(const T& t) {return SimTK::isNaN(t);} 00640 static bool isInf(const T& t) {return SimTK::isInf(t);} 00641 00642 static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();} 00643 00644 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b) 00645 { return SimTK::isNumericallyEqual(a,b); } 00646 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol) 00647 { return SimTK::isNumericallyEqual(a,b,tol); } 00648 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b) 00649 { return SimTK::isNumericallyEqual(a,b); } 00650 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol) 00651 { return SimTK::isNumericallyEqual(a,b,tol); } 00652 00653 static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);} 00654 static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00655 static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);} 00656 static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00657 static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);} 00658 static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00659 static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);} 00660 static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00661 00662 // The rest are the same as the real equivalents, with zero imaginary part. 00663 static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;} 00664 static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;} 00665 static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;} 00666 static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;} 00667 static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;} 00668 static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;} 00669 static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;} 00670 static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;} 00671 static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;} 00672 static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;} 00673 static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;} 00674 static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;} 00675 static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;} 00676 static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;} 00677 static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;} 00678 static const T& getE() {static const T c(NTraits<R>::getE()); return c;} 00679 static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;} 00680 static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;} 00681 static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;} 00682 static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;} 00683 static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;} 00684 static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;} 00685 static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;} 00686 static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;} 00687 static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;} 00688 static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;} 00689 }; 00690 00691 00692 // Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is 00693 // always just the complex type of sufficient precision. 00694 #define SimTK_BNTCMPLX_SPEC(T1,T2) \ 00695 template<> template<> struct NTraits< complex<T1> >::Result<T2> { \ 00696 typedef Widest< complex<T1>,T2 >::Type W; \ 00697 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00698 }; \ 00699 template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > { \ 00700 typedef Widest< complex<T1>,complex<T2> >::Type W; \ 00701 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00702 }; \ 00703 template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > { \ 00704 typedef Widest< complex<T1>,complex<T2> >::Type W; \ 00705 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00706 } 00707 SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);SimTK_BNTCMPLX_SPEC(float,long double); 00708 SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);SimTK_BNTCMPLX_SPEC(double,long double); 00709 SimTK_BNTCMPLX_SPEC(long double,float);SimTK_BNTCMPLX_SPEC(long double,double);SimTK_BNTCMPLX_SPEC(long double,long double); 00710 #undef SimTK_BNTCMPLX_SPEC 00711 00712 00713 // conjugate -- should be instantiated only for float, double, long double. 00714 template <class R> class NTraits< conjugate<R> > { 00715 typedef complex<R> C; 00716 public: 00717 typedef conjugate<R> T; 00718 typedef negator<T> TNeg; // type of this after *cast* to its negative 00719 typedef conjugate<R> TWithoutNegator; // type of this ignoring negator (there isn't one!) 00720 typedef R TReal; 00721 typedef negator<R> TImag; 00722 typedef conjugate<R> TComplex; 00723 typedef complex<R> THerm; // conjugate evaporates 00724 typedef conjugate<R> TPosTrans; // Positional transpose of scalar does nothing 00725 typedef R TSqHermT; // C*C~ 00726 typedef R TSqTHerm; // ~C*C (same) 00727 typedef conjugate<R> TElement; 00728 typedef conjugate<R> TRow; 00729 typedef conjugate<R> TCol; 00730 00731 typedef complex<R> TSqrt; 00732 typedef R TAbs; 00733 typedef complex<R> TStandard; 00734 typedef conjugate<R> TInvert; 00735 typedef conjugate<R> TNormalize; 00736 00737 typedef conjugate<R> Scalar; 00738 typedef conjugate<R> ULessScalar; 00739 typedef conjugate<R> Number; 00740 typedef complex<R> StdNumber; 00741 typedef R Precision; 00742 typedef R ScalarNormSq; 00743 00744 // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>) 00745 // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>) 00746 // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>) 00747 // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>) 00748 // Must specialize for P=real or P=complex or P=conjugate 00749 template <class P> struct Result { 00750 typedef typename CNT<P>::template Result<T>::Mul Mul; 00751 typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd; 00752 typedef typename CNT<P>::template Result<T>::Add Add; 00753 typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub; 00754 }; 00755 00756 // Shape-preserving element substitution (easy for scalars!) 00757 template <class P> struct Substitute { 00758 typedef P Type; 00759 }; 00760 00761 enum { 00762 NRows = 1, 00763 NCols = 1, 00764 RowSpacing = 1, 00765 ColSpacing = 1, 00766 NPackedElements = 1, // not two! 00767 NActualElements = 1, 00768 NActualScalars = 1, 00769 ImagOffset = 1, 00770 RealStrideFactor = 2, // double stride when casting to real or imaginary 00771 ArgDepth = SCALAR_DEPTH, 00772 IsScalar = 1, 00773 IsULessScalar = 1, 00774 IsNumber = 1, 00775 IsStdNumber = 0, 00776 IsPrecision = 0, 00777 SignInterpretation = 1 // after cast to Number, what is the sign? 00778 }; 00779 00780 static const T* getData(const T& t) { return &t; } 00781 static T* updData(T& t) { return &t; } 00782 static const TReal& real(const T& t) { return t.real(); } 00783 static TReal& real(T& t) { return t.real(); } 00784 static const TImag& imag(const T& t) { return t.imag(); } 00785 static TImag& imag(T& t) { return t.imag(); } 00786 00787 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} 00788 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} 00789 00790 static const THerm& transpose(const T& t) {return t.conj();} 00791 static THerm& transpose(T& t) {return t.conj();} 00792 00793 static const TPosTrans& positionalTranspose(const T& t) 00794 {return reinterpret_cast<const TPosTrans&>(t);} 00795 static TPosTrans& positionalTranspose(T& t) 00796 {return reinterpret_cast<TPosTrans&>(t);} 00797 00798 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) 00799 {return reinterpret_cast<const TWithoutNegator&>(t);} 00800 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) 00801 {return reinterpret_cast<TWithoutNegator&>(t);} 00802 00803 static ScalarNormSq scalarNormSqr(const T& t) 00804 { return t.real()*t.real() + t.negImag()*t.negImag(); } 00805 static TSqrt sqrt(const T& t) 00806 { return std::sqrt(C(t)); } // cast to complex (one negation) 00807 static TAbs abs(const T& t) 00808 { return std::abs(t.conj()); } // no, not just sqrt of scalarNormSqr()! 00809 static TStandard standardize(const T& t) 00810 { return TStandard(t); } // i.e., convert to complex 00811 static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));} 00812 00813 // 1/conj(z) = conj(1/z), for complex z. 00814 static TInvert invert(const T& t) 00815 { const typename NTraits<THerm>::TInvert cmplx(NTraits<THerm>::invert(t.conj())); 00816 return reinterpret_cast<const TInvert&>(cmplx); } // recast complex to conjugate it 00817 00818 // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should 00819 // be positive NaN. 00820 static const T& getNaN() { 00821 static const T c=T(NTraits<R>::getNaN(),NTraits<R>::getNaN()); 00822 return c; 00823 } 00824 // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals 00825 // are positive Inf. 00826 static const T& getInfinity() { 00827 static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity()); 00828 return c; 00829 } 00830 // But we want the constant i (=sqrt(-1)) to be the same however we represent it, 00831 // so for conjugate i = 0 - (-1)i. 00832 static const T& getI() { 00833 static const T c = T(0,-1); 00834 return c; 00835 } 00836 00837 static bool isFinite(const T& t) {return SimTK::isFinite(t);} 00838 static bool isNaN(const T& t) {return SimTK::isNaN(t);} 00839 static bool isInf(const T& t) {return SimTK::isInf(t);} 00840 00841 static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();} 00842 00843 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b) 00844 { return SimTK::isNumericallyEqual(a,b); } 00845 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol) 00846 { return SimTK::isNumericallyEqual(a,b,tol); } 00847 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b) 00848 { return SimTK::isNumericallyEqual(a,b); } 00849 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol) 00850 { return SimTK::isNumericallyEqual(a,b,tol); } 00851 00852 static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);} 00853 static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00854 static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);} 00855 static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00856 static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);} 00857 static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00858 static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);} 00859 static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00860 00861 // The rest are the same as the real equivalents, with zero imaginary part. 00862 static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;} 00863 static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;} 00864 static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;} 00865 static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;} 00866 static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;} 00867 static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;} 00868 static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;} 00869 static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;} 00870 static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;} 00871 static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;} 00872 static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;} 00873 static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;} 00874 static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;} 00875 static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;} 00876 static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;} 00877 static const T& getE() {static const T c(NTraits<R>::getE()); return c;} 00878 static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;} 00879 static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;} 00880 static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;} 00881 static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;} 00882 static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;} 00883 static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;} 00884 static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;} 00885 static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;} 00886 static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;} 00887 static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;} 00888 }; 00889 00890 // Any op involving conjugate & a real is best left as a conjugate. However, 00891 // an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost 00892 // and return just a complex in some cases. Also, we prefer negator<complex> to conjugate. 00893 // 00894 // Conj op complex 00895 // a-bi * r+si = (ar+bs) + (as-br)i (complex) 00896 // a-bi / r+si = hairy and slow anyway; we'll convert to complex 00897 // a-bi + r+si = (a+r) + (s-b)i (complex) 00898 // a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>) 00899 // 00900 // Conj op Conj 00901 // a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>) 00902 // a-bi / r-si = hairy and slow anyway; we'll convert to complex 00903 // a-bi + r-si = (a+r) - (b+s)i (conjugate) 00904 // a-bi - r-si = (a-r) + (s-b)i (complex) 00905 00906 #define SimTK_NTRAITS_CONJ_SPEC(T1,T2) \ 00907 template<> template<> struct NTraits< conjugate<T1> >::Result<T2> { \ 00908 typedef conjugate<Widest<T1,T2>::Type> W; \ 00909 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00910 }; \ 00911 template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\ 00912 typedef Widest<complex<T1>,complex<T2> >::Type W; \ 00913 typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub; \ 00914 }; \ 00915 template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\ 00916 typedef Widest<T1,T2>::Type W; typedef complex<W> WC; \ 00917 typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\ 00918 } 00919 SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double); 00920 SimTK_NTRAITS_CONJ_SPEC(float,long double); 00921 SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double); 00922 SimTK_NTRAITS_CONJ_SPEC(double,long double); 00923 SimTK_NTRAITS_CONJ_SPEC(long double,float);SimTK_NTRAITS_CONJ_SPEC(long double,double); 00924 SimTK_NTRAITS_CONJ_SPEC(long double,long double); 00925 #undef SimTK_NTRAITS_CONJ_SPEC 00926 00927 00928 // Specializations for real numbers. 00929 // For real scalar R, op result types are: 00930 // Typeof(R*P) = Typeof(P*R) 00931 // Typeof(R/P) = Typeof(inv(P)*R) 00932 // Typeof(R+P) = Typeof(P+R) 00933 // typeof(R-P) = Typeof(P::TNeg + R) 00934 // These must be specialized for P=Real and P=Complex. 00935 #define SimTK_DEFINE_REAL_NTRAITS(R) \ 00936 template <> class NTraits<R> { \ 00937 public: \ 00938 typedef R T; \ 00939 typedef negator<T> TNeg; \ 00940 typedef T TWithoutNegator; \ 00941 typedef T TReal; \ 00942 typedef T TImag; \ 00943 typedef complex<T> TComplex; \ 00944 typedef T THerm; \ 00945 typedef T TPosTrans; \ 00946 typedef T TSqHermT; \ 00947 typedef T TSqTHerm; \ 00948 typedef T TElement; \ 00949 typedef T TRow; \ 00950 typedef T TCol; \ 00951 typedef T TSqrt; \ 00952 typedef T TAbs; \ 00953 typedef T TStandard; \ 00954 typedef T TInvert; \ 00955 typedef T TNormalize; \ 00956 typedef T Scalar; \ 00957 typedef T ULessScalar; \ 00958 typedef T Number; \ 00959 typedef T StdNumber; \ 00960 typedef T Precision; \ 00961 typedef T ScalarNormSq; \ 00962 template <class P> struct Result { \ 00963 typedef typename CNT<P>::template Result<R>::Mul Mul; \ 00964 typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd; \ 00965 typedef typename CNT<P>::template Result<R>::Add Add; \ 00966 typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub; \ 00967 }; \ 00968 template <class P> struct Substitute { \ 00969 typedef P Type; \ 00970 }; \ 00971 enum { \ 00972 NRows = 1, \ 00973 NCols = 1, \ 00974 RowSpacing = 1, \ 00975 ColSpacing = 1, \ 00976 NPackedElements = 1, \ 00977 NActualElements = 1, \ 00978 NActualScalars = 1, \ 00979 ImagOffset = 0, \ 00980 RealStrideFactor = 1, \ 00981 ArgDepth = SCALAR_DEPTH, \ 00982 IsScalar = 1, \ 00983 IsULessScalar = 1, \ 00984 IsNumber = 1, \ 00985 IsStdNumber = 1, \ 00986 IsPrecision = 1, \ 00987 SignInterpretation = 1 \ 00988 }; \ 00989 static const T* getData(const T& t) { return &t; } \ 00990 static T* updData(T& t) { return &t; } \ 00991 static const T& real(const T& t) { return t; } \ 00992 static T& real(T& t) { return t; } \ 00993 static const T& imag(const T&) { return reinterpret_cast<const T&>(zeroes); } \ 00994 static T& imag(T&) { assert(false); return *reinterpret_cast<T*>(0); } \ 00995 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} \ 00996 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} \ 00997 static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} \ 00998 static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} \ 00999 static const TPosTrans& positionalTranspose(const T& t) \ 01000 {return reinterpret_cast<const TPosTrans&>(t);} \ 01001 static TPosTrans& positionalTranspose(T& t) \ 01002 {return reinterpret_cast<TPosTrans&>(t);} \ 01003 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) \ 01004 {return reinterpret_cast<const TWithoutNegator&>(t);} \ 01005 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) \ 01006 {return reinterpret_cast<TWithoutNegator&>(t);} \ 01007 static ScalarNormSq scalarNormSqr(const T& t) {return t*t;} \ 01008 static TSqrt sqrt(const T& t) {return std::sqrt(t);} \ 01009 static TAbs abs(const T& t) {return std::abs(t);} \ 01010 static const TStandard& standardize(const T& t) {return t;} \ 01011 static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));} \ 01012 static TInvert invert(const T& t) {return T(1)/t;} \ 01013 /* properties of this floating point representation, with memory addresses */ \ 01014 static const T& getEps() {return RTraits<T>::getEps();} \ 01015 static const T& getSignificant() {return RTraits<T>::getSignificant();} \ 01016 static const T& getNaN() {static const T c=std::numeric_limits<T>::quiet_NaN(); return c;} \ 01017 static const T& getInfinity() {static const T c=std::numeric_limits<T>::infinity(); return c;} \ 01018 static const T& getLeastPositive(){static const T c=std::numeric_limits<T>::min(); return c;} \ 01019 static const T& getMostPositive() {static const T c=std::numeric_limits<T>::max(); return c;} \ 01020 static const T& getLeastNegative(){static const T c=-std::numeric_limits<T>::min(); return c;} \ 01021 static const T& getMostNegative() {static const T c=-std::numeric_limits<T>::max(); return c;} \ 01022 static const T& getSqrtEps() {static const T c=std::sqrt(getEps()); return c;} \ 01023 static const T& getTiny() {static const T c=std::pow(getEps(), (T)1.25L); return c;} \ 01024 static bool isFinite(const T& t) {return SimTK::isFinite(t);} \ 01025 static bool isNaN (const T& t) {return SimTK::isNaN(t);} \ 01026 static bool isInf (const T& t) {return SimTK::isInf(t);} \ 01027 /* Methods to use for approximate comparisons. Perform comparison in the wider of the two */ \ 01028 /* precisions, using the default tolerance from the narrower of the two precisions. */ \ 01029 static double getDefaultTolerance() {return RTraits<T>::getDefaultTolerance();} \ 01030 static bool isNumericallyEqual(const T& t, const float& f) {return SimTK::isNumericallyEqual(t,f);} \ 01031 static bool isNumericallyEqual(const T& t, const double& d) {return SimTK::isNumericallyEqual(t,d);} \ 01032 static bool isNumericallyEqual(const T& t, const long double& l) {return SimTK::isNumericallyEqual(t,l);} \ 01033 static bool isNumericallyEqual(const T& t, int i) {return SimTK::isNumericallyEqual(t,i);} \ 01034 /* Here the tolerance is given so we don't have to figure it out. */ \ 01035 static bool isNumericallyEqual(const T& t, const float& f, double tol){return SimTK::isNumericallyEqual(t,f,tol);} \ 01036 static bool isNumericallyEqual(const T& t, const double& d, double tol){return SimTK::isNumericallyEqual(t,d,tol);} \ 01037 static bool isNumericallyEqual(const T& t, const long double& l, double tol){return SimTK::isNumericallyEqual(t,l,tol);} \ 01038 static bool isNumericallyEqual(const T& t, int i, double tol){return SimTK::isNumericallyEqual(t,i,tol);} \ 01039 /* Carefully calculated constants with convenient memory addresses. */ \ 01040 static const T& getZero() {static const T c=(T)(0); return c;} \ 01041 static const T& getOne() {static const T c=(T)(1); return c;} \ 01042 static const T& getMinusOne() {static const T c=(T)(-1); return c;} \ 01043 static const T& getTwo() {static const T c=(T)(2); return c;} \ 01044 static const T& getThree() {static const T c=(T)(3); return c;} \ 01045 static const T& getOneHalf() {static const T c=(T)(0.5L); return c;} \ 01046 static const T& getOneThird() {static const T c=(T)(1.L/3.L); return c;} \ 01047 static const T& getOneFourth() {static const T c=(T)(0.25L); return c;} \ 01048 static const T& getOneFifth() {static const T c=(T)(0.2L); return c;} \ 01049 static const T& getOneSixth() {static const T c=(T)(1.L/6.L); return c;} \ 01050 static const T& getOneSeventh() {static const T c=(T)(1.L/7.L); return c;} \ 01051 static const T& getOneEighth() {static const T c=(T)(0.125L); return c;} \ 01052 static const T& getOneNinth() {static const T c=(T)(1.L/9.L); return c;} \ 01053 static const T& getPi() {static const T c=(T)(SimTK_PI); return c;} \ 01054 static const T& getOneOverPi() {static const T c=(T)(1.L/SimTK_PI); return c;} \ 01055 static const T& getE() {static const T c=(T)(SimTK_E); return c;} \ 01056 static const T& getLog2E() {static const T c=(T)(SimTK_LOG2E); return c;} \ 01057 static const T& getLog10E() {static const T c=(T)(SimTK_LOG10E); return c;} \ 01058 static const T& getSqrt2() {static const T c=(T)(SimTK_SQRT2); return c;} \ 01059 static const T& getOneOverSqrt2() {static const T c=(T)(1.L/SimTK_SQRT2); return c;} \ 01060 static const T& getSqrt3() {static const T c=(T)(SimTK_SQRT3); return c;} \ 01061 static const T& getOneOverSqrt3() {static const T c=(T)(1.L/SimTK_SQRT3); return c;} \ 01062 static const T& getCubeRoot2() {static const T c=(T)(SimTK_CBRT2); return c;} \ 01063 static const T& getCubeRoot3() {static const T c=(T)(SimTK_CBRT3); return c;} \ 01064 static const T& getLn2() {static const T c=(T)(SimTK_LN2); return c;} \ 01065 static const T& getLn10() {static const T c=(T)(SimTK_LN10); return c;} \ 01066 /* integer digit counts useful for formatted input and output */ \ 01067 static const int getNumDigits() {static const int c=(int)(std::log10(1/getEps()) -0.5); return c;} \ 01068 static const int getLosslessNumDigits() {static const int c=(int)(std::log10(1/getTiny())+0.5); return c;} \ 01069 }; \ 01070 template<> struct NTraits<R>::Result<float> \ 01071 {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01072 template<> struct NTraits<R>::Result<double> \ 01073 {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01074 template<> struct NTraits<R>::Result<long double> \ 01075 {typedef Widest<R,long double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01076 template<> struct NTraits<R>::Result<complex<float> > \ 01077 {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01078 template<> struct NTraits<R>::Result<complex<double> > \ 01079 {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01080 template<> struct NTraits<R>::Result<complex<long double> > \ 01081 {typedef Widest<R,complex<long double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01082 template<> struct NTraits<R>::Result<conjugate<float> > \ 01083 {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01084 template<> struct NTraits<R>::Result<conjugate<double> > \ 01085 {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01086 template<> struct NTraits<R>::Result<conjugate<long double> > \ 01087 {typedef conjugate<Widest<R,long double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;} 01088 SimTK_DEFINE_REAL_NTRAITS(float); 01089 SimTK_DEFINE_REAL_NTRAITS(double); 01090 SimTK_DEFINE_REAL_NTRAITS(long double); 01091 #undef SimTK_DEFINE_REAL_NTRAITS 01092 01094 template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { }; 01095 template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { }; 01096 template <> class CNT<float> : public NTraits<float> { }; 01097 template <> class CNT<double> : public NTraits<double> { }; 01098 template <> class CNT<long double> : public NTraits<long double> { }; 01099 01100 01101 } // namespace SimTK 01102 01103 #endif //SimTK_SIMMATRIX_NTRAITS_H_
