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