Simbody
3.5
|
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 00103 template <class R1, class R2> struct Widest {/* Only defined for built-ins. */}; 00104 template <> struct Widest<float,float> {typedef float Type; typedef float Precision;}; 00105 template <> struct Widest<float,double> {typedef double Type; typedef double Precision;}; 00106 template <> struct Widest<float,long double> {typedef long double Type; typedef long double Precision;}; 00107 template <> struct Widest<double,float> {typedef double Type; typedef double Precision;}; 00108 template <> struct Widest<double,double> {typedef double Type; typedef double Precision;}; 00109 template <> struct Widest<double,long double> {typedef long double Type; typedef long double Precision;}; 00110 template <> struct Widest<long double,float> {typedef long double Type; typedef long double Precision;}; 00111 template <> struct Widest<long double,double> {typedef long double Type; typedef long double Precision;}; 00112 template <> struct Widest<long double,long double> {typedef long double Type; typedef long double Precision;}; 00113 template <class R1, class R2> struct Widest< complex<R1>,complex<R2> > { 00114 typedef complex< typename Widest<R1,R2>::Type > Type; 00115 typedef typename Widest<R1,R2>::Precision Precision; 00116 }; 00117 template <class R1, class R2> struct Widest< complex<R1>,R2 > { 00118 typedef complex< typename Widest<R1,R2>::Type > Type; 00119 typedef typename Widest<R1,R2>::Precision Precision; 00120 }; 00121 template <class R1, class R2> struct Widest< R1,complex<R2> > { 00122 typedef complex< typename Widest<R1,R2>::Type > Type; 00123 typedef typename Widest<R1,R2>::Precision Precision; 00124 }; 00125 00135 template <class R1, class R2> struct Narrowest {/* Only defined for built-ins. */}; 00136 template <> struct Narrowest<float,float> {typedef float Type; typedef float Precision;}; 00137 template <> struct Narrowest<float,double> {typedef float Type; typedef float Precision;}; 00138 template <> struct Narrowest<float,long double> {typedef float Type; typedef float Precision;}; 00139 template <> struct Narrowest<double,float> {typedef float Type; typedef float Precision;}; 00140 template <> struct Narrowest<double,double> {typedef double Type; typedef double Precision;}; 00141 template <> struct Narrowest<double,long double> {typedef double Type; typedef double Precision;}; 00142 template <> struct Narrowest<long double,float> {typedef float Type; typedef float Precision;}; 00143 template <> struct Narrowest<long double,double> {typedef double Type; typedef double Precision;}; 00144 template <> struct Narrowest<long double,long double> {typedef long double Type; typedef long double Precision;}; 00145 template <class R1, class R2> struct Narrowest< complex<R1>,complex<R2> > { 00146 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00147 typedef typename Narrowest<R1,R2>::Precision Precision; 00148 }; 00149 template <class R1, class R2> struct Narrowest< complex<R1>,R2 > { 00150 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00151 typedef typename Narrowest<R1,R2>::Precision Precision; 00152 }; 00153 template <class R1, class R2> struct Narrowest< R1,complex<R2> > { 00154 typedef complex< typename Narrowest<R1,R2>::Type > Type; 00155 typedef typename Narrowest<R1,R2>::Precision Precision; 00156 }; 00157 00159 template <class R> class RTraits {/* Only defined for real types */}; 00160 template <> class RTraits<float> { 00161 public: 00163 static const float& getEps() {static const float c=std::numeric_limits<float>::epsilon(); return c;} 00165 static const float& getSignificant() {static const float c=std::pow(getEps(), 0.875f); return c;} 00167 static double getDefaultTolerance() {return (double)getSignificant();} 00168 }; 00169 template <> class RTraits<double> { 00170 public: 00171 static const double& getEps() {static const double c=std::numeric_limits<double>::epsilon(); return c;} 00172 static const double& getSignificant() {static const double c=std::pow(getEps(), 0.875); return c;} 00173 static double getDefaultTolerance() {return getSignificant();} 00174 }; 00175 template <> class RTraits<long double> { 00176 public: 00177 static const long double& getEps() {static const long double c=std::numeric_limits<long double>::epsilon(); return c;} 00178 static const long double& getSignificant() {static const long double c=std::pow(getEps(), 0.875L); return c;} 00179 static double getDefaultTolerance() {return (double)getSignificant();} 00180 }; 00181 00198 // See negator.h for isNaN() applied to negated scalars. 00200 inline bool isNaN(const float& x) {return std::isnan(x);} 00201 inline bool isNaN(const double& x) {return std::isnan(x);} 00202 inline bool isNaN(const long double& x) {return std::isnan(x);} 00203 00204 template <class P> inline bool 00205 isNaN(const std::complex<P>& x) 00206 { return isNaN(x.real()) || isNaN(x.imag());} 00207 00208 template <class P> inline bool 00209 isNaN(const conjugate<P>& x) 00210 { return isNaN(x.real()) || isNaN(x.negImag());} 00212 00225 // See negator.h for isFinite() applied to negated scalars. 00227 inline bool isFinite(const float& x) {return std::isfinite(x);} 00228 inline bool isFinite(const double& x) {return std::isfinite(x);} 00229 inline bool isFinite(const long double& x) {return std::isfinite(x);} 00230 00231 template <class P> inline bool 00232 isFinite(const std::complex<P>& x) 00233 { return isFinite(x.real()) && isFinite(x.imag());} 00234 00235 template <class P> inline bool 00236 isFinite(const conjugate<P>& x) 00237 { return isFinite(x.real()) && isFinite(x.negImag());} 00239 00254 // See negator.h for isInf() applied to negated scalars. 00256 inline bool isInf(const float& x) {return std::isinf(x);} 00257 inline bool isInf(const double& x) {return std::isinf(x);} 00258 inline bool isInf(const long double& x) {return std::isinf(x);} 00259 00260 template <class P> inline bool 00261 isInf(const std::complex<P>& x) { 00262 return (isInf(x.real()) && !isNaN(x.imag())) 00263 || (isInf(x.imag()) && !isNaN(x.real())); 00264 } 00265 00266 template <class P> inline bool 00267 isInf(const conjugate<P>& x) { 00268 return (isInf(x.real()) && !isNaN(x.negImag())) 00269 || (isInf(x.negImag()) && !isNaN(x.real())); 00270 } 00272 00312 00313 inline bool isNumericallyEqual(const float& a, const float& b, 00314 double tol = RTraits<float>::getDefaultTolerance()) 00315 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00316 const float scale = std::max(std::max(std::abs(a),std::abs(b)), 1.f); 00317 return std::abs(a-b) <= scale*(float)tol; } 00319 inline bool isNumericallyEqual(const double& a, const double& b, 00320 double tol = RTraits<double>::getDefaultTolerance()) 00321 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00322 const double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.); 00323 return std::abs(a-b) <= scale*tol; } 00325 inline bool isNumericallyEqual(const long double& a, const long double& b, 00326 double tol = RTraits<long double>::getDefaultTolerance()) 00327 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false; 00328 const long double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.L); 00329 return std::abs(a-b) <= scale*(long double)tol; } 00330 00332 inline bool isNumericallyEqual(const float& a, const double& b, 00333 double tol = RTraits<float>::getDefaultTolerance()) 00334 { return isNumericallyEqual((double)a, b, tol); } 00336 inline bool isNumericallyEqual(const double& a, const float& b, 00337 double tol = RTraits<float>::getDefaultTolerance()) 00338 { return isNumericallyEqual(a, (double)b, tol); } 00340 inline bool isNumericallyEqual(const float& a, const long double& b, 00341 double tol = RTraits<float>::getDefaultTolerance()) 00342 { return isNumericallyEqual((long double)a, b, tol); } 00344 inline bool isNumericallyEqual(const long double& a, const float& b, 00345 double tol = RTraits<float>::getDefaultTolerance()) 00346 { return isNumericallyEqual(a, (long double)b, tol); } 00348 inline bool isNumericallyEqual(const double& a, const long double& b, 00349 double tol = RTraits<double>::getDefaultTolerance()) 00350 { return isNumericallyEqual((long double)a, b, tol); } 00352 inline bool isNumericallyEqual(const long double& a, const double& b, 00353 double tol = RTraits<double>::getDefaultTolerance()) 00354 { return isNumericallyEqual(a, (long double)b, tol); } 00355 00357 inline bool isNumericallyEqual(const float& a, int b, 00358 double tol = RTraits<float>::getDefaultTolerance()) 00359 { return isNumericallyEqual(a, (double)b, tol); } 00361 inline bool isNumericallyEqual(int a, const float& b, 00362 double tol = RTraits<float>::getDefaultTolerance()) 00363 { return isNumericallyEqual((double)a, b, tol); } 00365 inline bool isNumericallyEqual(const double& a, int b, 00366 double tol = RTraits<double>::getDefaultTolerance()) 00367 { return isNumericallyEqual(a, (double)b, tol); } 00369 inline bool isNumericallyEqual(int a, const double& b, 00370 double tol = RTraits<double>::getDefaultTolerance()) 00371 { return isNumericallyEqual((double)a, b, tol); } 00373 inline bool isNumericallyEqual(const long double& a, int b, 00374 double tol = RTraits<long double>::getDefaultTolerance()) 00375 { return isNumericallyEqual(a, (long double)b, tol); } 00377 inline bool isNumericallyEqual(int a, const long double& b, 00378 double tol = RTraits<long double>::getDefaultTolerance()) 00379 { return isNumericallyEqual((long double)a, b, tol); } 00380 00384 template <class P, class Q> 00385 inline bool isNumericallyEqual 00386 ( const std::complex<P>& a, const std::complex<Q>& b, 00387 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00388 { return isNumericallyEqual(a.real(),b.real(),tol) 00389 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00390 00394 template <class P, class Q> 00395 inline bool isNumericallyEqual 00396 ( const conjugate<P>& a, const conjugate<Q>& b, 00397 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00398 { return isNumericallyEqual(a.real(),b.real(),tol) 00399 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00400 00404 template <class P, class Q> 00405 inline bool isNumericallyEqual 00406 ( const std::complex<P>& a, const conjugate<Q>& b, 00407 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00408 { return isNumericallyEqual(a.real(),b.real(),tol) 00409 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00410 00414 template <class P, class Q> 00415 inline bool isNumericallyEqual 00416 ( const conjugate<P>& a, const std::complex<Q>& b, 00417 double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance()) 00418 { return isNumericallyEqual(a.real(),b.real(),tol) 00419 && isNumericallyEqual(a.imag(),b.imag(),tol); } 00420 00422 template <class P> inline bool 00423 isNumericallyEqual(const std::complex<P>& a, const float& b, 00424 double tol = RTraits<float>::getDefaultTolerance()) 00425 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); } 00427 template <class P> inline bool 00428 isNumericallyEqual(const float& a, const std::complex<P>& b, 00429 double tol = RTraits<float>::getDefaultTolerance()) 00430 { return isNumericallyEqual(b,a,tol); } 00432 template <class P> inline bool 00433 isNumericallyEqual(const std::complex<P>& a, const double& b, 00434 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00435 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); } 00437 template <class P> inline bool 00438 isNumericallyEqual(const double& a, const std::complex<P>& b, 00439 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00440 { return isNumericallyEqual(b,a,tol); } 00442 template <class P> inline bool 00443 isNumericallyEqual(const std::complex<P>& a, const long double& b, 00444 double tol = RTraits<P>::getDefaultTolerance()) 00445 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); } 00447 template <class P> inline bool 00448 isNumericallyEqual(const long double& a, const std::complex<P>& b, 00449 double tol = RTraits<P>::getDefaultTolerance()) 00450 { return isNumericallyEqual(b,a,tol); } 00452 template <class P> inline bool 00453 isNumericallyEqual(const std::complex<P>& a, int b, 00454 double tol = RTraits<P>::getDefaultTolerance()) 00455 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); } 00457 template <class P> inline bool 00458 isNumericallyEqual(int a, const std::complex<P>& b, 00459 double tol = RTraits<P>::getDefaultTolerance()) 00460 { return isNumericallyEqual(b,a,tol); } 00461 00463 template <class P> inline bool 00464 isNumericallyEqual(const conjugate<P>& a, const float& b, 00465 double tol = RTraits<float>::getDefaultTolerance()) 00466 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); } 00468 template <class P> inline bool 00469 isNumericallyEqual(const float& a, const conjugate<P>& b, 00470 double tol = RTraits<float>::getDefaultTolerance()) 00471 { return isNumericallyEqual(b,a,tol); } 00473 template <class P> inline bool 00474 isNumericallyEqual(const conjugate<P>& a, const double& b, 00475 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00476 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); } 00478 template <class P> inline bool 00479 isNumericallyEqual(const double& a, const conjugate<P>& b, 00480 double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance()) 00481 { return isNumericallyEqual(b,a,tol); } 00483 template <class P> inline bool 00484 isNumericallyEqual(const conjugate<P>& a, const long double& b, 00485 double tol = RTraits<P>::getDefaultTolerance()) 00486 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); } 00488 template <class P> inline bool 00489 isNumericallyEqual(const long double& a, const conjugate<P>& b, 00490 double tol = RTraits<P>::getDefaultTolerance()) 00491 { return isNumericallyEqual(b,a,tol); } 00493 template <class P> inline bool 00494 isNumericallyEqual(const conjugate<P>& a, int b, 00495 double tol = RTraits<P>::getDefaultTolerance()) 00496 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); } 00498 template <class P> inline bool 00499 isNumericallyEqual(int a, const conjugate<P>& b, 00500 double tol = RTraits<P>::getDefaultTolerance()) 00501 { return isNumericallyEqual(b,a,tol); } 00502 00504 00505 00506 template <class N> class NTraits { 00507 // only the specializations below are allowed 00508 }; 00509 00512 template <class R> class NTraits< complex<R> > { 00513 typedef complex<R> C; 00514 public: 00515 typedef C T; 00516 typedef negator<C> TNeg; // type of this after *cast* to its negative 00517 typedef C TWithoutNegator; // type of this ignoring negator (there isn't one!) 00518 00519 typedef R TReal; 00520 typedef R TImag; 00521 typedef C TComplex; 00522 typedef conjugate<R> THerm; // this is a recast 00523 typedef C TPosTrans; 00524 typedef R TSqHermT; // ~C*C 00525 typedef R TSqTHerm; // C*~C (same) 00526 typedef C TElement; 00527 typedef C TRow; 00528 typedef C TCol; 00529 00530 typedef C TSqrt; 00531 typedef R TAbs; 00532 typedef C TStandard; // complex is a standard type 00533 typedef C TInvert; // this is a calculation, so use standard number 00534 typedef C TNormalize; 00535 00536 typedef C Scalar; 00537 typedef C ULessScalar; 00538 typedef C Number; 00539 typedef C StdNumber; 00540 typedef R Precision; 00541 typedef R ScalarNormSq; 00542 00543 // For complex scalar C, op result types are: 00544 // Typeof(C*P) = Typeof(P*C) 00545 // Typeof(C/P) = Typeof(inv(P)*C) 00546 // Typeof(C+P) = Typeof(P+C) 00547 // typeof(C-P) = Typeof(P::TNeg + C) 00548 // These must be specialized for P=real, complex, conjugate. 00549 template <class P> struct Result { 00550 typedef typename CNT<P>::template Result<C>::Mul Mul; 00551 typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd; 00552 typedef typename CNT<P>::template Result<C>::Add Add; 00553 typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub; 00554 }; 00555 00556 // Shape-preserving element substitution (easy for scalars!) 00557 template <class P> struct Substitute { 00558 typedef P Type; 00559 }; 00560 00561 enum { 00562 NRows = 1, 00563 NCols = 1, 00564 RowSpacing = 1, 00565 ColSpacing = 1, 00566 NPackedElements = 1, // not two! 00567 NActualElements = 1, 00568 NActualScalars = 1, 00569 ImagOffset = 1, 00570 RealStrideFactor = 2, // double stride when casting to real or imaginary 00571 ArgDepth = SCALAR_DEPTH, 00572 IsScalar = 1, 00573 IsULessScalar = 1, 00574 IsNumber = 1, 00575 IsStdNumber = 1, 00576 IsPrecision = 0, 00577 SignInterpretation = 1 // after cast to Number, what is the sign? 00578 }; 00579 static const T* getData(const T& t) { return &t; } 00580 static T* updData(T& t) { return &t; } 00581 static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; } 00582 static R& real(T& t) { return (reinterpret_cast<R*>(&t))[0]; } 00583 static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; } 00584 static R& imag(T& t) { return (reinterpret_cast<R*>(&t))[1]; } 00585 00586 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} 00587 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} 00588 00589 static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} 00590 static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} 00591 00592 static const TPosTrans& positionalTranspose(const T& t) 00593 {return reinterpret_cast<const TPosTrans&>(t);} 00594 static TPosTrans& positionalTranspose(T& t) 00595 {return reinterpret_cast<TPosTrans&>(t);} 00596 00597 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) 00598 {return reinterpret_cast<const TWithoutNegator&>(t);} 00599 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) 00600 {return reinterpret_cast<TWithoutNegator&>(t);} 00601 00602 static ScalarNormSq scalarNormSqr(const T& t) 00603 { return t.real()*t.real() + t.imag()*t.imag(); } 00604 static TSqrt sqrt(const T& t) 00605 { return std::sqrt(t); } 00606 static TAbs abs(const T& t) 00607 { return std::abs(t); } // no, not just sqrt of scalarNormSqr()! 00608 static const TStandard& standardize(const T& t) {return t;} // already standard 00609 static TNormalize normalize(const T& t) {return t/abs(t);} 00610 static TInvert invert(const T& t) {return TReal(1)/t;} 00611 00612 static const T& getNaN() { 00613 static const T c=T(NTraits<R>::getNaN(), NTraits<R>::getNaN()); 00614 return c; 00615 } 00616 static const T& getInfinity() { 00617 static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity()); 00618 return c; 00619 } 00620 00621 static const T& getI() { 00622 static const T c = T(0,1); 00623 return c; 00624 } 00625 00626 static bool isFinite(const T& t) {return SimTK::isFinite(t);} 00627 static bool isNaN(const T& t) {return SimTK::isNaN(t);} 00628 static bool isInf(const T& t) {return SimTK::isInf(t);} 00629 00630 static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();} 00631 00632 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b) 00633 { return SimTK::isNumericallyEqual(a,b); } 00634 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol) 00635 { return SimTK::isNumericallyEqual(a,b,tol); } 00636 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b) 00637 { return SimTK::isNumericallyEqual(a,b); } 00638 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol) 00639 { return SimTK::isNumericallyEqual(a,b,tol); } 00640 00641 static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);} 00642 static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00643 static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);} 00644 static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00645 static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);} 00646 static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00647 static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);} 00648 static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00649 00650 // The rest are the same as the real equivalents, with zero imaginary part. 00651 static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;} 00652 static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;} 00653 static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;} 00654 static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;} 00655 static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;} 00656 static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;} 00657 static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;} 00658 static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;} 00659 static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;} 00660 static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;} 00661 static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;} 00662 static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;} 00663 static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;} 00664 static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;} 00665 static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;} 00666 static const T& getE() {static const T c(NTraits<R>::getE()); return c;} 00667 static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;} 00668 static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;} 00669 static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;} 00670 static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;} 00671 static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;} 00672 static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;} 00673 static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;} 00674 static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;} 00675 static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;} 00676 static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;} 00677 }; 00678 00679 00680 // Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is 00681 // always just the complex type of sufficient precision. 00682 #define SimTK_BNTCMPLX_SPEC(T1,T2) \ 00683 template<> template<> struct NTraits< complex<T1> >::Result<T2> { \ 00684 typedef Widest< complex<T1>,T2 >::Type W; \ 00685 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00686 }; \ 00687 template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > { \ 00688 typedef Widest< complex<T1>,complex<T2> >::Type W; \ 00689 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00690 }; \ 00691 template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > { \ 00692 typedef Widest< complex<T1>,complex<T2> >::Type W; \ 00693 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00694 } 00695 SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);SimTK_BNTCMPLX_SPEC(float,long double); 00696 SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);SimTK_BNTCMPLX_SPEC(double,long double); 00697 SimTK_BNTCMPLX_SPEC(long double,float);SimTK_BNTCMPLX_SPEC(long double,double);SimTK_BNTCMPLX_SPEC(long double,long double); 00698 #undef SimTK_BNTCMPLX_SPEC 00699 00700 00701 // conjugate -- should be instantiated only for float, double, long double. 00702 template <class R> class NTraits< conjugate<R> > { 00703 typedef complex<R> C; 00704 public: 00705 typedef conjugate<R> T; 00706 typedef negator<T> TNeg; // type of this after *cast* to its negative 00707 typedef conjugate<R> TWithoutNegator; // type of this ignoring negator (there isn't one!) 00708 typedef R TReal; 00709 typedef negator<R> TImag; 00710 typedef conjugate<R> TComplex; 00711 typedef complex<R> THerm; // conjugate evaporates 00712 typedef conjugate<R> TPosTrans; // Positional transpose of scalar does nothing 00713 typedef R TSqHermT; // C*C~ 00714 typedef R TSqTHerm; // ~C*C (same) 00715 typedef conjugate<R> TElement; 00716 typedef conjugate<R> TRow; 00717 typedef conjugate<R> TCol; 00718 00719 typedef complex<R> TSqrt; 00720 typedef R TAbs; 00721 typedef complex<R> TStandard; 00722 typedef conjugate<R> TInvert; 00723 typedef conjugate<R> TNormalize; 00724 00725 typedef conjugate<R> Scalar; 00726 typedef conjugate<R> ULessScalar; 00727 typedef conjugate<R> Number; 00728 typedef complex<R> StdNumber; 00729 typedef R Precision; 00730 typedef R ScalarNormSq; 00731 00732 // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>) 00733 // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>) 00734 // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>) 00735 // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>) 00736 // Must specialize for P=real or P=complex or P=conjugate 00737 template <class P> struct Result { 00738 typedef typename CNT<P>::template Result<T>::Mul Mul; 00739 typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd; 00740 typedef typename CNT<P>::template Result<T>::Add Add; 00741 typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub; 00742 }; 00743 00744 // Shape-preserving element substitution (easy for scalars!) 00745 template <class P> struct Substitute { 00746 typedef P Type; 00747 }; 00748 00749 enum { 00750 NRows = 1, 00751 NCols = 1, 00752 RowSpacing = 1, 00753 ColSpacing = 1, 00754 NPackedElements = 1, // not two! 00755 NActualElements = 1, 00756 NActualScalars = 1, 00757 ImagOffset = 1, 00758 RealStrideFactor = 2, // double stride when casting to real or imaginary 00759 ArgDepth = SCALAR_DEPTH, 00760 IsScalar = 1, 00761 IsULessScalar = 1, 00762 IsNumber = 1, 00763 IsStdNumber = 0, 00764 IsPrecision = 0, 00765 SignInterpretation = 1 // after cast to Number, what is the sign? 00766 }; 00767 00768 static const T* getData(const T& t) { return &t; } 00769 static T* updData(T& t) { return &t; } 00770 static const TReal& real(const T& t) { return t.real(); } 00771 static TReal& real(T& t) { return t.real(); } 00772 static const TImag& imag(const T& t) { return t.imag(); } 00773 static TImag& imag(T& t) { return t.imag(); } 00774 00775 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} 00776 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} 00777 00778 static const THerm& transpose(const T& t) {return t.conj();} 00779 static THerm& transpose(T& t) {return t.conj();} 00780 00781 static const TPosTrans& positionalTranspose(const T& t) 00782 {return reinterpret_cast<const TPosTrans&>(t);} 00783 static TPosTrans& positionalTranspose(T& t) 00784 {return reinterpret_cast<TPosTrans&>(t);} 00785 00786 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) 00787 {return reinterpret_cast<const TWithoutNegator&>(t);} 00788 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) 00789 {return reinterpret_cast<TWithoutNegator&>(t);} 00790 00791 static ScalarNormSq scalarNormSqr(const T& t) 00792 { return t.real()*t.real() + t.negImag()*t.negImag(); } 00793 static TSqrt sqrt(const T& t) 00794 { return std::sqrt(C(t)); } // cast to complex (one negation) 00795 static TAbs abs(const T& t) 00796 { return std::abs(t.conj()); } // no, not just sqrt of scalarNormSqr()! 00797 static TStandard standardize(const T& t) 00798 { return TStandard(t); } // i.e., convert to complex 00799 static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));} 00800 00801 // 1/conj(z) = conj(1/z), for complex z. 00802 static TInvert invert(const T& t) 00803 { const typename NTraits<THerm>::TInvert cmplx(NTraits<THerm>::invert(t.conj())); 00804 return reinterpret_cast<const TInvert&>(cmplx); } // recast complex to conjugate it 00805 00806 // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should 00807 // be positive NaN. 00808 static const T& getNaN() { 00809 static const T c=T(NTraits<R>::getNaN(),NTraits<R>::getNaN()); 00810 return c; 00811 } 00812 // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals 00813 // are positive Inf. 00814 static const T& getInfinity() { 00815 static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity()); 00816 return c; 00817 } 00818 // But we want the constant i (=sqrt(-1)) to be the same however we represent it, 00819 // so for conjugate i = 0 - (-1)i. 00820 static const T& getI() { 00821 static const T c = T(0,-1); 00822 return c; 00823 } 00824 00825 static bool isFinite(const T& t) {return SimTK::isFinite(t);} 00826 static bool isNaN(const T& t) {return SimTK::isNaN(t);} 00827 static bool isInf(const T& t) {return SimTK::isInf(t);} 00828 00829 static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();} 00830 00831 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b) 00832 { return SimTK::isNumericallyEqual(a,b); } 00833 template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol) 00834 { return SimTK::isNumericallyEqual(a,b,tol); } 00835 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b) 00836 { return SimTK::isNumericallyEqual(a,b); } 00837 template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol) 00838 { return SimTK::isNumericallyEqual(a,b,tol); } 00839 00840 static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);} 00841 static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00842 static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);} 00843 static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00844 static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);} 00845 static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00846 static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);} 00847 static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);} 00848 00849 // The rest are the same as the real equivalents, with zero imaginary part. 00850 static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;} 00851 static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;} 00852 static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;} 00853 static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;} 00854 static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;} 00855 static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;} 00856 static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;} 00857 static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;} 00858 static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;} 00859 static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;} 00860 static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;} 00861 static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;} 00862 static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;} 00863 static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;} 00864 static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;} 00865 static const T& getE() {static const T c(NTraits<R>::getE()); return c;} 00866 static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;} 00867 static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;} 00868 static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;} 00869 static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;} 00870 static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;} 00871 static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;} 00872 static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;} 00873 static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;} 00874 static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;} 00875 static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;} 00876 }; 00877 00878 // Any op involving conjugate & a real is best left as a conjugate. However, 00879 // an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost 00880 // and return just a complex in some cases. Also, we prefer negator<complex> to conjugate. 00881 // 00882 // Conj op complex 00883 // a-bi * r+si = (ar+bs) + (as-br)i (complex) 00884 // a-bi / r+si = hairy and slow anyway; we'll convert to complex 00885 // a-bi + r+si = (a+r) + (s-b)i (complex) 00886 // a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>) 00887 // 00888 // Conj op Conj 00889 // a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>) 00890 // a-bi / r-si = hairy and slow anyway; we'll convert to complex 00891 // a-bi + r-si = (a+r) - (b+s)i (conjugate) 00892 // a-bi - r-si = (a-r) + (s-b)i (complex) 00893 00894 #define SimTK_NTRAITS_CONJ_SPEC(T1,T2) \ 00895 template<> template<> struct NTraits< conjugate<T1> >::Result<T2> { \ 00896 typedef conjugate<Widest<T1,T2>::Type> W; \ 00897 typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \ 00898 }; \ 00899 template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\ 00900 typedef Widest<complex<T1>,complex<T2> >::Type W; \ 00901 typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub; \ 00902 }; \ 00903 template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\ 00904 typedef Widest<T1,T2>::Type W; typedef complex<W> WC; \ 00905 typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\ 00906 } 00907 SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double); 00908 SimTK_NTRAITS_CONJ_SPEC(float,long double); 00909 SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double); 00910 SimTK_NTRAITS_CONJ_SPEC(double,long double); 00911 SimTK_NTRAITS_CONJ_SPEC(long double,float);SimTK_NTRAITS_CONJ_SPEC(long double,double); 00912 SimTK_NTRAITS_CONJ_SPEC(long double,long double); 00913 #undef SimTK_NTRAITS_CONJ_SPEC 00914 00915 00916 // Specializations for real numbers. 00917 // For real scalar R, op result types are: 00918 // Typeof(R*P) = Typeof(P*R) 00919 // Typeof(R/P) = Typeof(inv(P)*R) 00920 // Typeof(R+P) = Typeof(P+R) 00921 // typeof(R-P) = Typeof(P::TNeg + R) 00922 // These must be specialized for P=Real and P=Complex. 00923 #define SimTK_DEFINE_REAL_NTRAITS(R) \ 00924 template <> class NTraits<R> { \ 00925 public: \ 00926 typedef R T; \ 00927 typedef negator<T> TNeg; \ 00928 typedef T TWithoutNegator; \ 00929 typedef T TReal; \ 00930 typedef T TImag; \ 00931 typedef complex<T> TComplex; \ 00932 typedef T THerm; \ 00933 typedef T TPosTrans; \ 00934 typedef T TSqHermT; \ 00935 typedef T TSqTHerm; \ 00936 typedef T TElement; \ 00937 typedef T TRow; \ 00938 typedef T TCol; \ 00939 typedef T TSqrt; \ 00940 typedef T TAbs; \ 00941 typedef T TStandard; \ 00942 typedef T TInvert; \ 00943 typedef T TNormalize; \ 00944 typedef T Scalar; \ 00945 typedef T ULessScalar; \ 00946 typedef T Number; \ 00947 typedef T StdNumber; \ 00948 typedef T Precision; \ 00949 typedef T ScalarNormSq; \ 00950 template <class P> struct Result { \ 00951 typedef typename CNT<P>::template Result<R>::Mul Mul; \ 00952 typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd; \ 00953 typedef typename CNT<P>::template Result<R>::Add Add; \ 00954 typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub; \ 00955 }; \ 00956 template <class P> struct Substitute { \ 00957 typedef P Type; \ 00958 }; \ 00959 enum { \ 00960 NRows = 1, \ 00961 NCols = 1, \ 00962 RowSpacing = 1, \ 00963 ColSpacing = 1, \ 00964 NPackedElements = 1, \ 00965 NActualElements = 1, \ 00966 NActualScalars = 1, \ 00967 ImagOffset = 0, \ 00968 RealStrideFactor = 1, \ 00969 ArgDepth = SCALAR_DEPTH, \ 00970 IsScalar = 1, \ 00971 IsULessScalar = 1, \ 00972 IsNumber = 1, \ 00973 IsStdNumber = 1, \ 00974 IsPrecision = 1, \ 00975 SignInterpretation = 1 \ 00976 }; \ 00977 static const T* getData(const T& t) { return &t; } \ 00978 static T* updData(T& t) { return &t; } \ 00979 static const T& real(const T& t) { return t; } \ 00980 static T& real(T& t) { return t; } \ 00981 static const T& imag(const T&) { return getZero(); } \ 00982 static T& imag(T&) { assert(false); return *reinterpret_cast<T*>(0); } \ 00983 static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} \ 00984 static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} \ 00985 static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} \ 00986 static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} \ 00987 static const TPosTrans& positionalTranspose(const T& t) \ 00988 {return reinterpret_cast<const TPosTrans&>(t);} \ 00989 static TPosTrans& positionalTranspose(T& t) \ 00990 {return reinterpret_cast<TPosTrans&>(t);} \ 00991 static const TWithoutNegator& castAwayNegatorIfAny(const T& t) \ 00992 {return reinterpret_cast<const TWithoutNegator&>(t);} \ 00993 static TWithoutNegator& updCastAwayNegatorIfAny(T& t) \ 00994 {return reinterpret_cast<TWithoutNegator&>(t);} \ 00995 static ScalarNormSq scalarNormSqr(const T& t) {return t*t;} \ 00996 static TSqrt sqrt(const T& t) {return std::sqrt(t);} \ 00997 static TAbs abs(const T& t) {return std::abs(t);} \ 00998 static const TStandard& standardize(const T& t) {return t;} \ 00999 static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));} \ 01000 static TInvert invert(const T& t) {return T(1)/t;} \ 01001 /* properties of this floating point representation, with memory addresses */ \ 01002 static const T& getEps() {return RTraits<T>::getEps();} \ 01003 static const T& getSignificant() {return RTraits<T>::getSignificant();} \ 01004 static const T& getNaN() {static const T c=std::numeric_limits<T>::quiet_NaN(); return c;} \ 01005 static const T& getInfinity() {static const T c=std::numeric_limits<T>::infinity(); return c;} \ 01006 static const T& getLeastPositive(){static const T c=std::numeric_limits<T>::min(); return c;} \ 01007 static const T& getMostPositive() {static const T c=std::numeric_limits<T>::max(); return c;} \ 01008 static const T& getLeastNegative(){static const T c=-std::numeric_limits<T>::min(); return c;} \ 01009 static const T& getMostNegative() {static const T c=-std::numeric_limits<T>::max(); return c;} \ 01010 static const T& getSqrtEps() {static const T c=std::sqrt(getEps()); return c;} \ 01011 static const T& getTiny() {static const T c=std::pow(getEps(), (T)1.25L); return c;} \ 01012 static bool isFinite(const T& t) {return SimTK::isFinite(t);} \ 01013 static bool isNaN (const T& t) {return SimTK::isNaN(t);} \ 01014 static bool isInf (const T& t) {return SimTK::isInf(t);} \ 01015 /* Methods to use for approximate comparisons. Perform comparison in the wider of the two */ \ 01016 /* precisions, using the default tolerance from the narrower of the two precisions. */ \ 01017 static double getDefaultTolerance() {return RTraits<T>::getDefaultTolerance();} \ 01018 static bool isNumericallyEqual(const T& t, const float& f) {return SimTK::isNumericallyEqual(t,f);} \ 01019 static bool isNumericallyEqual(const T& t, const double& d) {return SimTK::isNumericallyEqual(t,d);} \ 01020 static bool isNumericallyEqual(const T& t, const long double& l) {return SimTK::isNumericallyEqual(t,l);} \ 01021 static bool isNumericallyEqual(const T& t, int i) {return SimTK::isNumericallyEqual(t,i);} \ 01022 /* Here the tolerance is given so we don't have to figure it out. */ \ 01023 static bool isNumericallyEqual(const T& t, const float& f, double tol){return SimTK::isNumericallyEqual(t,f,tol);} \ 01024 static bool isNumericallyEqual(const T& t, const double& d, double tol){return SimTK::isNumericallyEqual(t,d,tol);} \ 01025 static bool isNumericallyEqual(const T& t, const long double& l, double tol){return SimTK::isNumericallyEqual(t,l,tol);} \ 01026 static bool isNumericallyEqual(const T& t, int i, double tol){return SimTK::isNumericallyEqual(t,i,tol);} \ 01027 /* Carefully calculated constants with convenient memory addresses. */ \ 01028 static const T& getZero() {static const T c=(T)(0); return c;} \ 01029 static const T& getOne() {static const T c=(T)(1); return c;} \ 01030 static const T& getMinusOne() {static const T c=(T)(-1); return c;} \ 01031 static const T& getTwo() {static const T c=(T)(2); return c;} \ 01032 static const T& getThree() {static const T c=(T)(3); return c;} \ 01033 static const T& getOneHalf() {static const T c=(T)(0.5L); return c;} \ 01034 static const T& getOneThird() {static const T c=(T)(1.L/3.L); return c;} \ 01035 static const T& getOneFourth() {static const T c=(T)(0.25L); return c;} \ 01036 static const T& getOneFifth() {static const T c=(T)(0.2L); return c;} \ 01037 static const T& getOneSixth() {static const T c=(T)(1.L/6.L); return c;} \ 01038 static const T& getOneSeventh() {static const T c=(T)(1.L/7.L); return c;} \ 01039 static const T& getOneEighth() {static const T c=(T)(0.125L); return c;} \ 01040 static const T& getOneNinth() {static const T c=(T)(1.L/9.L); return c;} \ 01041 static const T& getPi() {static const T c=(T)(SimTK_PI); return c;} \ 01042 static const T& getOneOverPi() {static const T c=(T)(1.L/SimTK_PI); return c;} \ 01043 static const T& getE() {static const T c=(T)(SimTK_E); return c;} \ 01044 static const T& getLog2E() {static const T c=(T)(SimTK_LOG2E); return c;} \ 01045 static const T& getLog10E() {static const T c=(T)(SimTK_LOG10E); return c;} \ 01046 static const T& getSqrt2() {static const T c=(T)(SimTK_SQRT2); return c;} \ 01047 static const T& getOneOverSqrt2() {static const T c=(T)(1.L/SimTK_SQRT2); return c;} \ 01048 static const T& getSqrt3() {static const T c=(T)(SimTK_SQRT3); return c;} \ 01049 static const T& getOneOverSqrt3() {static const T c=(T)(1.L/SimTK_SQRT3); return c;} \ 01050 static const T& getCubeRoot2() {static const T c=(T)(SimTK_CBRT2); return c;} \ 01051 static const T& getCubeRoot3() {static const T c=(T)(SimTK_CBRT3); return c;} \ 01052 static const T& getLn2() {static const T c=(T)(SimTK_LN2); return c;} \ 01053 static const T& getLn10() {static const T c=(T)(SimTK_LN10); return c;} \ 01054 /* integer digit counts useful for formatted input and output */ \ 01055 static int getNumDigits() {static const int c=(int)(std::log10(1/getEps()) -0.5); return c;} \ 01056 static int getLosslessNumDigits() {static const int c=(int)(std::log10(1/getTiny())+0.5); return c;} \ 01057 }; \ 01058 template<> struct NTraits<R>::Result<float> \ 01059 {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01060 template<> struct NTraits<R>::Result<double> \ 01061 {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01062 template<> struct NTraits<R>::Result<long double> \ 01063 {typedef Widest<R,long double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01064 template<> struct NTraits<R>::Result<complex<float> > \ 01065 {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01066 template<> struct NTraits<R>::Result<complex<double> > \ 01067 {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01068 template<> struct NTraits<R>::Result<complex<long double> > \ 01069 {typedef Widest<R,complex<long double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01070 template<> struct NTraits<R>::Result<conjugate<float> > \ 01071 {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01072 template<> struct NTraits<R>::Result<conjugate<double> > \ 01073 {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \ 01074 template<> struct NTraits<R>::Result<conjugate<long double> > \ 01075 {typedef conjugate<Widest<R,long double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;} 01076 SimTK_DEFINE_REAL_NTRAITS(float); 01077 SimTK_DEFINE_REAL_NTRAITS(double); 01078 SimTK_DEFINE_REAL_NTRAITS(long double); 01079 #undef SimTK_DEFINE_REAL_NTRAITS 01080 01082 template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { }; 01083 template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { }; 01084 template <> class CNT<float> : public NTraits<float> { }; 01085 template <> class CNT<double> : public NTraits<double> { }; 01086 template <> class CNT<long double> : public NTraits<long double> { }; 01087 01088 01089 } // namespace SimTK 01090 01091 #endif //SimTK_SIMMATRIX_NTRAITS_H_