Simbody  3.5
NTraits.h
Go to the documentation of this file.
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_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines