NTraits.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_NTRAITS_H_
00002 #define SimTK_SIMMATRIX_NTRAITS_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Simmatrix(tm)                       *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core biosimulation toolkit originating from      *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-7 Stanford University and the Authors.         *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00058 #include <cstddef>
00059 #include <cassert>
00060 #include <complex>
00061 #include <iostream>
00062 #include <limits>
00063 
00064 using std::complex;
00065     
00066 namespace SimTK {
00067 
00068 // This is the 3rd type of number, conjugate. It is like complex except that
00069 // the represented value is the conjugate of the value represented by a complex
00070 // number containing the same bit pattern. That is, complex and conjugate both
00071 // contain two real numbers, re and im, with complex(re,im) meaning
00072 // (re + im*i) while conjugate(re,im) means (re - im*i). It is guaranteed that
00073 // our conjugate type and complex type have identical sizes and representations.
00074 // Together, these defininitions and guarantees permit conjugation
00075 // to be done by reinterpretation rather than be computation.
00076 template <class R> class conjugate;    // Only defined for float, double, long double
00077 
00078 // Specializations of this class provide information about Composite Numerical Types
00079 // in the style of std::numeric_limits<T>. It is specialized for the numeric types
00080 // but can be invoked on any composite numerical type as well.
00081 template <class T> class CNT;
00082 
00083 // NTraits provides information like CNT, but for numeric types only.
00084 template <class N> class NTraits; // allowed only for N=<number>
00085 template <class R> class NTraits< complex<R> >;
00086 template <class R> class NTraits< conjugate<R> >;
00087 template <> class NTraits<float>;
00088 template <> class NTraits<double>;
00089 template <> class NTraits<long double>;
00090 
00091 // This is an adaptor for numeric types which negates the apparent values. A
00092 // negator<N> has exactly the same internal representation as a numeric
00093 // value N, but it is to be interpreted has having the negative of the value
00094 // it would have if interpreted as an N. This permits negation to be done
00095 // by reinterpretation rather than compuation. A full set of arithmetic operators
00096 // are provided involving negator<N>'s and N's. Sometimes we can save an op or
00097 // two this way. For example negator<N>*negator<N> can be performed as an N*N
00098 // since the negations cancel, and we saved two floating point negations.
00099 template <class N> class negator;      // Only defined for numbers
00100 
00101 // This is here so we can provide references to 0's when needed, e.g.
00102 // when returning the imaginary part of a real number. These are local to
00103 // the compilation unit, so the returned addresses will differ in different
00104 // files. There are enough zeroes here for the widest number, complex<long double>
00105 // (or conjugate<long double>).
00106 static const complex<long double> zeroes(0);
00107 
00108 // This class is specialized for all 36 combinations of standard types
00109 // (that is, real and complex types in each of three precisions)
00110 // and has a single typedef "Type" which is the appropriate "widened"
00111 // type for use when R1 & R2 appear in an operation together. For example,
00112 // if R1=complex<float> and R2=long double, Widest<R1,R2>::Type is
00113 // complex<long double>.
00114 template <class R1, class R2> struct Widest {/* Only defined for built-ins. */};
00115 template <> struct Widest<float,float>              {typedef float Type;};
00116 template <> struct Widest<float,double>             {typedef double Type;};
00117 template <> struct Widest<float,long double>        {typedef long double Type;};
00118 template <> struct Widest<double,float>             {typedef double Type;};
00119 template <> struct Widest<double,double>            {typedef double Type;};
00120 template <> struct Widest<double,long double>       {typedef long double Type;};
00121 template <> struct Widest<long double,float>        {typedef long double Type;};
00122 template <> struct Widest<long double,double>       {typedef long double Type;};
00123 template <> struct Widest<long double,long double>  {typedef long double Type;};
00124 template <class R1, class R2> struct Widest< complex<R1>,complex<R2> >
00125   { typedef complex< typename Widest<R1,R2>::Type > Type; };
00126 template <class R1, class R2> struct Widest< complex<R1>,R2 >
00127   { typedef complex< typename Widest<R1,R2>::Type > Type; };
00128 template <class R1, class R2> struct Widest< R1,complex<R2> >
00129   { typedef complex< typename Widest<R1,R2>::Type > Type; };
00130 
00131 template <class N> class NTraits { 
00132     // only the specializations below are allowed 
00133 };
00134 
00137 template <class R> class NTraits< complex<R> > {
00138     typedef complex<R>  C;
00139 public:
00140     typedef C                T;
00141     typedef negator<C>       TNeg;            // type of this after *cast* to its negative
00142     typedef C                TWithoutNegator; // type of this ignoring negator (there isn't one!)
00143 
00144     typedef R                TReal;
00145     typedef R                TImag;
00146     typedef C                TComplex;    
00147     typedef conjugate<R>     THerm;     // this is a recast
00148     typedef C                TPosTrans;
00149     typedef R                TSqHermT;  // ~C*C
00150     typedef R                TSqTHerm;  // C*~C (same)
00151     typedef C                TElement;
00152     typedef C                TRow;
00153     typedef C                TCol;
00154 
00155     typedef R                TAbs;
00156     typedef C                TStandard; // complex is a standard type
00157     typedef C                TInvert;   // this is a calculation, so use standard number
00158     typedef C                TNormalize;
00159 
00160     typedef C                Scalar;
00161     typedef C                Number;
00162     typedef C                StdNumber;
00163     typedef R                Precision;
00164     typedef R                ScalarSq;
00165 
00166     // For complex scalar C, op result types are:
00167     //   Typeof(C*P) = Typeof(P*C)
00168     //   Typeof(C/P) = Typeof(inv(P)*C)
00169     //   Typeof(C+P) = Typeof(P+C)
00170     //   typeof(C-P) = Typeof(P::TNeg + C)
00171     // These must be specialized for P=real, complex, conjugate.
00172     template <class P> struct Result { 
00173         typedef typename CNT<P>::template Result<C>::Mul Mul;
00174         typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd;
00175         typedef typename CNT<P>::template Result<C>::Add Add;
00176         typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub;
00177     };
00178 
00179     // Shape-preserving element substitution (easy for scalars!)
00180     template <class P> struct Substitute {
00181         typedef P Type;
00182     };
00183 
00184     enum {
00185         NRows               = 1,
00186         NCols               = 1,
00187         RowSpacing          = 1,
00188         ColSpacing          = 1,
00189         NPackedElements     = 1,      // not two!
00190         NActualElements     = 1,
00191         NActualScalars      = 1,
00192         ImagOffset          = 1,
00193         RealStrideFactor    = 2,      // double stride when casting to real or imaginary
00194         ArgDepth            = SCALAR_DEPTH,
00195         IsScalar            = 1,
00196         IsNumber            = 1,
00197         IsStdNumber         = 1,
00198         IsPrecision         = 0,
00199         SignInterpretation  = 1       // after cast to Number, what is the sign?
00200     }; 
00201     static const T* getData(const T& t) { return &t; } 
00202     static T*       updData(T& t)       { return &t; }
00203     static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; }
00204     static R&       real(T& t)       { return (reinterpret_cast<R*>(&t))[0]; }
00205     static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; }
00206     static R&       imag(T& t)       { return (reinterpret_cast<R*>(&t))[1]; }
00207 
00208     static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
00209     static       TNeg& negate(T& t)       {return reinterpret_cast<TNeg&>(t);}
00210 
00211     static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}
00212     static       THerm& transpose(T& t)       {return reinterpret_cast<THerm&>(t);}
00213 
00214     static const TPosTrans& positionalTranspose(const T& t)
00215         {return reinterpret_cast<const TPosTrans&>(t);}
00216     static       TPosTrans& positionalTranspose(T& t)
00217         {return reinterpret_cast<TPosTrans&>(t);} 
00218 
00219     static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
00220         {return reinterpret_cast<const TWithoutNegator&>(t);}
00221     static       TWithoutNegator& updCastAwayNegatorIfAny(const T& t)
00222         {return reinterpret_cast<TWithoutNegator&>(t);}
00223 
00224     static ScalarSq scalarNormSqr(const T& t)
00225         { return t.real()*t.real() + t.imag()*t.imag(); }
00226     static TAbs     abs(const T& t)
00227         { return std::abs(t); } // no, not just sqrt of scalarNormSqr()!
00228     static const TStandard& standardize(const T& t) {return t;} // already standard
00229     static TNormalize normalize(const T& t) {return t/abs(t);}
00230     static TInvert    invert(const T& t)    {return TReal(1)/t;}
00231 
00232     static const T& getNaN() {
00233         static const T c=T(NTraits<R>::getNaN(), NTraits<R>::getNaN());
00234         return c;
00235     }
00236     static const T& getInfinity() {
00237         static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
00238         return c;
00239     }
00240     static const T& getI() {
00241         static const T c = T(0,1);
00242         return c;
00243     }
00244 
00245     // The rest are the same as the real equivalents, with zero imaginary part.              
00246     static const T& getZero()         {static const T c(NTraits<R>::getZero());         return c;}
00247     static const T& getOne()          {static const T c(NTraits<R>::getOne());          return c;}
00248     static const T& getMinusOne()     {static const T c(NTraits<R>::getMinusOne());     return c;}
00249     static const T& getTwo()          {static const T c(NTraits<R>::getTwo());          return c;}
00250     static const T& getThree()        {static const T c(NTraits<R>::getThree());        return c;}
00251     static const T& getOneHalf()      {static const T c(NTraits<R>::getOneHalf());      return c;}
00252     static const T& getOneThird()     {static const T c(NTraits<R>::getOneThird());     return c;}
00253     static const T& getOneFourth()    {static const T c(NTraits<R>::getOneFourth());    return c;}
00254     static const T& getOneFifth()     {static const T c(NTraits<R>::getOneFifth());     return c;}
00255     static const T& getOneSixth()     {static const T c(NTraits<R>::getOneSixth());     return c;}
00256     static const T& getOneSeventh()   {static const T c(NTraits<R>::getOneSeventh());   return c;}
00257     static const T& getOneEighth()    {static const T c(NTraits<R>::getOneEighth());    return c;}
00258     static const T& getOneNinth()     {static const T c(NTraits<R>::getOneNinth());     return c;}
00259     static const T& getPi()           {static const T c(NTraits<R>::getPi());           return c;}
00260     static const T& getOneOverPi()    {static const T c(NTraits<R>::getOneOverPi());    return c;}
00261     static const T& getE()            {static const T c(NTraits<R>::getE());            return c;}
00262     static const T& getLog2E()        {static const T c(NTraits<R>::getLog2E());        return c;}
00263     static const T& getLog10E()       {static const T c(NTraits<R>::getLog10E());       return c;}
00264     static const T& getSqrt2()        {static const T c(NTraits<R>::getSqrt2());        return c;}
00265     static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
00266     static const T& getSqrt3()        {static const T c(NTraits<R>::getSqrt3());        return c;}
00267     static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
00268     static const T& getCubeRoot2()    {static const T c(NTraits<R>::getCubeRoot2());    return c;}
00269     static const T& getCubeRoot3()    {static const T c(NTraits<R>::getCubeRoot3());    return c;}
00270     static const T& getLn2()          {static const T c(NTraits<R>::getLn2());          return c;}
00271     static const T& getLn10()         {static const T c(NTraits<R>::getLn10());         return c;}
00272 };
00273 
00274 
00275 // Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is
00276 // always just the complex type of sufficient precision.
00277 #define SimTK_BNTCMPLX_SPEC(T1,T2)  \
00278 template<> template<> struct NTraits< complex<T1> >::Result<T2> {      \
00279     typedef Widest< complex<T1>,T2 >::Type W;                      \
00280     typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
00281 };                                                                      \
00282 template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > {  \
00283     typedef Widest< complex<T1>,complex<T2> >::Type W;        \
00284     typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
00285 };                                                                      \
00286 template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > {  \
00287     typedef Widest< complex<T1>,complex<T2> >::Type W;        \
00288     typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
00289 }
00290 SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);SimTK_BNTCMPLX_SPEC(float,long double);
00291 SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);SimTK_BNTCMPLX_SPEC(double,long double);
00292 SimTK_BNTCMPLX_SPEC(long double,float);SimTK_BNTCMPLX_SPEC(long double,double);SimTK_BNTCMPLX_SPEC(long double,long double);
00293 #undef SimTK_BNTCMPLX_SPEC
00294 
00295 
00296 // conjugate -- should be instantiated only for float, double, long double.
00297 template <class R> class NTraits< conjugate<R> > {
00298     typedef complex<R>          C;
00299 public:
00300     typedef conjugate<R>        T;
00301     typedef negator<T>          TNeg;            // type of this after *cast* to its negative
00302     typedef T                   TWithoutNegator; // type of this ignoring negator (there isn't one!)
00303     typedef R                   TReal;
00304     typedef negator<R>          TImag;
00305     typedef conjugate<R>        TComplex;     
00306     typedef complex<R>          THerm;      // conjugate evaporates
00307     typedef conjugate<R>        TPosTrans;  // Positional transpose of scalar does nothing
00308     typedef R                   TSqHermT;   // C*C~
00309     typedef R                   TSqTHerm;   // ~C*C (same)
00310     typedef conjugate<R>        TElement;
00311     typedef conjugate<R>        TRow;
00312     typedef conjugate<R>        TCol;
00313 
00314     typedef R                   TAbs;
00315     typedef C                   TStandard;
00316     typedef complex<R>          TInvert;
00317     typedef T                   TNormalize;
00318 
00319     typedef conjugate<R>        Scalar;
00320     typedef conjugate<R>        Number;
00321     typedef C                   StdNumber;
00322     typedef R                   Precision;
00323     typedef R                   ScalarSq;
00324 
00325     // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>)
00326     // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>)
00327     // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>)
00328     // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>)
00329     // Must specialize for P=real or P=complex or P=conjugate
00330     template <class P> struct Result {
00331         typedef typename CNT<P>::template Result<T>::Mul Mul;
00332         typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd;
00333         typedef typename CNT<P>::template Result<T>::Add Add;
00334         typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub;
00335     };
00336 
00337     // Shape-preserving element substitution (easy for scalars!)
00338     template <class P> struct Substitute {
00339         typedef P Type;
00340     };
00341 
00342     enum {
00343         NRows               = 1,
00344         NCols               = 1,
00345         RowSpacing          = 1,
00346         ColSpacing          = 1,
00347         NPackedElements     = 1,      // not two!
00348         NActualElements     = 1,
00349         NActualScalars      = 1,
00350         ImagOffset          = 1,
00351         RealStrideFactor    = 2,      // double stride when casting to real or imaginary
00352         ArgDepth            = SCALAR_DEPTH,
00353         IsScalar            = 1,
00354         IsNumber            = 1,
00355         IsStdNumber         = 0,
00356         IsPrecision         = 0,
00357         SignInterpretation  = 1       // after cast to Number, what is the sign?
00358     }; 
00359 
00360     static const T*     getData(const T& t) { return &t; } 
00361     static T*           updData(T& t)       { return &t; }
00362     static const TReal& real(const T& t) { return t.real(); }
00363     static TReal&       real(T& t)       { return t.real(); }
00364     static const TImag& imag(const T& t) { return t.imag(); }
00365     static TImag&       imag(T& t)       { return t.imag(); }
00366 
00367     static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
00368     static       TNeg& negate(T& t)       {return reinterpret_cast<TNeg&>(t);}
00369 
00370     static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}
00371     static       THerm& transpose(T& t)       {return reinterpret_cast<THerm&>(t);}
00372 
00373     static const TPosTrans& positionalTranspose(const T& t)
00374         {return reinterpret_cast<const TPosTrans&>(t);}
00375     static       TPosTrans& positionalTranspose(T& t)
00376         {return reinterpret_cast<TPosTrans&>(t);} 
00377 
00378     static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
00379         {return reinterpret_cast<const TWithoutNegator&>(t);}
00380     static       TWithoutNegator& updCastAwayNegatorIfAny(const T& t)
00381         {return reinterpret_cast<TWithoutNegator&>(t);}
00382 
00383     static ScalarSq scalarNormSqr(const T& t)
00384         { return t.real()*t.real() + t.negImag()*t.negImag(); }
00385     static TAbs     abs(const T& t)
00386         { return std::abs(t.conj()); }  // no, not just sqrt of scalarNormSqr()!
00387     static TStandard standardize(const T& t)
00388         { return TStandard(t); }        // i.e., convert to complex
00389     static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));}
00390     static TInvert    invert(const T& t)    {return TReal(1)/t;}
00391 
00392     // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should
00393     // be positive NaN.
00394     static const T& getNaN() { 
00395         static const T c=T(NTraits<R>::getNaN(),NTraits<R>::getNaN());
00396         return c;
00397     }
00398     // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals
00399     // are positive Inf.
00400     static const T& getInfinity() {
00401         static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
00402         return c;
00403     }
00404     // But we want the constant i (=sqrt(-1)) to be the same however we represent it,
00405     // so for conjugate i = 0 - (-1)i.
00406     static const T& getI() {
00407         static const T c = T(0,-1);
00408         return c;
00409     }
00410 
00411     // The rest are the same as the real equivalents, with zero imaginary part.              
00412     static const T& getZero()         {static const T c(NTraits<R>::getZero());         return c;}
00413     static const T& getOne()          {static const T c(NTraits<R>::getOne());          return c;}
00414     static const T& getMinusOne()     {static const T c(NTraits<R>::getMinusOne());     return c;}
00415     static const T& getTwo()          {static const T c(NTraits<R>::getTwo());          return c;}
00416     static const T& getThree()        {static const T c(NTraits<R>::getThree());        return c;}
00417     static const T& getOneHalf()      {static const T c(NTraits<R>::getOneHalf());      return c;}
00418     static const T& getOneThird()     {static const T c(NTraits<R>::getOneThird());     return c;}
00419     static const T& getOneFourth()    {static const T c(NTraits<R>::getOneFourth());    return c;}
00420     static const T& getOneFifth()     {static const T c(NTraits<R>::getOneFifth());     return c;}
00421     static const T& getOneSixth()     {static const T c(NTraits<R>::getOneSixth());     return c;}
00422     static const T& getOneSeventh()   {static const T c(NTraits<R>::getOneSeventh());   return c;}
00423     static const T& getOneEighth()    {static const T c(NTraits<R>::getOneEighth());    return c;}
00424     static const T& getOneNinth()     {static const T c(NTraits<R>::getOneNinth());     return c;}
00425     static const T& getPi()           {static const T c(NTraits<R>::getPi());           return c;}
00426     static const T& getOneOverPi()    {static const T c(NTraits<R>::getOneOverPi());    return c;}
00427     static const T& getE()            {static const T c(NTraits<R>::getE());            return c;}
00428     static const T& getLog2E()        {static const T c(NTraits<R>::getLog2E());        return c;}
00429     static const T& getLog10E()       {static const T c(NTraits<R>::getLog10E());       return c;}
00430     static const T& getSqrt2()        {static const T c(NTraits<R>::getSqrt2());        return c;}
00431     static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
00432     static const T& getSqrt3()        {static const T c(NTraits<R>::getSqrt3());        return c;}
00433     static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
00434     static const T& getCubeRoot2()    {static const T c(NTraits<R>::getCubeRoot2());    return c;}
00435     static const T& getCubeRoot3()    {static const T c(NTraits<R>::getCubeRoot3());    return c;}
00436     static const T& getLn2()          {static const T c(NTraits<R>::getLn2());          return c;}
00437     static const T& getLn10()         {static const T c(NTraits<R>::getLn10());         return c;}
00438 };
00439 
00440 // Any op involving conjugate & a real is best left as a conjugate. However,
00441 // an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost
00442 // and return just a complex in some cases. Also, we prefer negator<complex> to conjugate.
00443 //
00444 // Conj op complex
00445 //   a-bi * r+si = (ar+bs) + (as-br)i               (complex)
00446 //   a-bi / r+si = hairy and slow anyway; we'll convert to complex
00447 //   a-bi + r+si = (a+r) + (s-b)i                   (complex)
00448 //   a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>)
00449 //
00450 // Conj op Conj
00451 //   a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>)
00452 //   a-bi / r-si = hairy and slow anyway; we'll convert to complex
00453 //   a-bi + r-si = (a+r) - (b+s)i  (conjugate)
00454 //   a-bi - r-si = (a-r) + (s-b)i  (complex)
00455 
00456 #define SimTK_NTRAITS_CONJ_SPEC(T1,T2)                                      \
00457 template<> template<> struct NTraits< conjugate<T1> >::Result<T2> {         \
00458   typedef conjugate<Widest<T1,T2>::Type> W;                                 \
00459   typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;               \
00460 };                                                                          \
00461 template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\
00462   typedef Widest<complex<T1>,complex<T2> >::Type W;               \
00463   typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub;      \
00464 };                                                                          \
00465 template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\
00466     typedef Widest<T1,T2>::Type W; typedef complex<W> WC;              \
00467     typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\
00468 }
00469 SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double);
00470 SimTK_NTRAITS_CONJ_SPEC(float,long double);
00471 SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double);
00472 SimTK_NTRAITS_CONJ_SPEC(double,long double);
00473 SimTK_NTRAITS_CONJ_SPEC(long double,float);SimTK_NTRAITS_CONJ_SPEC(long double,double);
00474 SimTK_NTRAITS_CONJ_SPEC(long double,long double);
00475 #undef SimTK_NTRAITS_CONJ_SPEC 
00476 
00477 
00478 // Specializations for real numbers.
00479 // For real scalar R, op result types are:
00480 //   Typeof(R*P) = Typeof(P*R)
00481 //   Typeof(R/P) = Typeof(inv(P)*R)
00482 //   Typeof(R+P) = Typeof(P+R)
00483 //   typeof(R-P) = Typeof(P::TNeg + R)
00484 // These must be specialized for P=Real and P=Complex.
00485 #define SimTK_DEFINE_REAL_NTRAITS(R)            \
00486 template <> class NTraits<R> {                  \
00487 public:                                         \
00488     typedef R                T;                 \
00489     typedef negator<T>       TNeg;              \
00490     typedef T                TWithoutNegator;   \
00491     typedef T                TReal;             \
00492     typedef T                TImag;             \
00493     typedef complex<T>       TComplex;          \
00494     typedef T                THerm;             \
00495     typedef T                TPosTrans;         \
00496     typedef T                TSqHermT;          \
00497     typedef T                TSqTHerm;          \
00498     typedef T                TElement;          \
00499     typedef T                TRow;              \
00500     typedef T                TCol;              \
00501     typedef T                TAbs;              \
00502     typedef T                TStandard;         \
00503     typedef T                TInvert;           \
00504     typedef T                TNormalize;        \
00505     typedef T                Scalar;            \
00506     typedef T                Number;            \
00507     typedef T                StdNumber;         \
00508     typedef T                Precision;         \
00509     typedef T                ScalarSq;          \
00510     template <class P> struct Result {          \
00511         typedef typename CNT<P>::template Result<R>::Mul Mul;   \
00512         typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd;    \
00513         typedef typename CNT<P>::template Result<R>::Add Add;   \
00514         typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub;     \
00515     };                                          \
00516     template <class P> struct Substitute {      \
00517         typedef P Type;                         \
00518     };                                          \
00519     enum {                                      \
00520         NRows               = 1,                \
00521         NCols               = 1,                \
00522         RowSpacing          = 1,                \
00523         ColSpacing          = 1,                \
00524         NPackedElements     = 1,                \
00525         NActualElements     = 1,                \
00526         NActualScalars      = 1,                \
00527         ImagOffset          = 0,                \
00528         RealStrideFactor    = 1,                \
00529         ArgDepth            = SCALAR_DEPTH,     \
00530         IsScalar            = 1,                \
00531         IsNumber            = 1,                \
00532         IsStdNumber         = 1,                \
00533         IsPrecision         = 1,                \
00534         SignInterpretation  = 1                 \
00535     };                                          \
00536     static const T* getData(const T& t) { return &t; }  \
00537     static T*       updData(T& t)       { return &t; }  \
00538     static const T& real(const T& t) { return t; }      \
00539     static T&       real(T& t)       { return t; }      \
00540     static const T& imag(const T&)   { return reinterpret_cast<const T&>(zeroes); }   \
00541     static T&       imag(T&)         { assert(false); return *reinterpret_cast<T*>(0); } \
00542     static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}        \
00543     static       TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}                    \
00544     static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}   \
00545     static       THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);}               \
00546     static const TPosTrans& positionalTranspose(const T& t)                 \
00547         {return reinterpret_cast<const TPosTrans&>(t);}                     \
00548     static       TPosTrans& positionalTranspose(T& t)                       \
00549         {return reinterpret_cast<TPosTrans&>(t);}                           \
00550     static const TWithoutNegator& castAwayNegatorIfAny(const T& t)          \
00551         {return reinterpret_cast<const TWithoutNegator&>(t);}               \
00552     static       TWithoutNegator& updCastAwayNegatorIfAny(T& t)             \
00553         {return reinterpret_cast<TWithoutNegator&>(t);}                     \
00554     static ScalarSq   scalarNormSqr(const T& t) {return t*t;}               \
00555     static TAbs       abs(const T& t) {return std::abs(t);}                 \
00556     static const TStandard& standardize(const T& t) {return t;}             \
00557     static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));} \
00558     static TInvert invert(const T& t) {return T(1)/t;}                      \
00559     /* properties of this floating point representation, with memory addresses */     \
00560     static const T& getNaN()          {static const T c=std::numeric_limits<T>::quiet_NaN(); return c;} \
00561     static const T& getInfinity()     {static const T c=std::numeric_limits<T>::infinity();  return c;} \
00562     static const T& getEps()          {static const T c=std::numeric_limits<T>::epsilon();   return c;} \
00563     static const T& getLeastPositive(){static const T c=std::numeric_limits<T>::min();       return c;} \
00564     static const T& getMostPositive() {static const T c=std::numeric_limits<T>::max();       return c;} \
00565     static const T& getLeastNegative(){static const T c=-std::numeric_limits<T>::min();      return c;} \
00566     static const T& getMostNegative() {static const T c=-std::numeric_limits<T>::max();      return c;} \
00567     static const T& getSqrtEps()      {static const T c=std::sqrt(getEps());                 return c;} \
00568     static const T& getTiny()         {static const T c=std::pow(getEps(), (T)1.25L);        return c;} \
00569     static const T& getSignificant()  {static const T c=std::pow(getEps(), (T)0.875L);       return c;} \
00570     /* carefully calculated constants with convenient memory addresses */                \
00571     static const T& getZero()         {static const T c=(T)(0);               return c;} \
00572     static const T& getOne()          {static const T c=(T)(1);               return c;} \
00573     static const T& getMinusOne()     {static const T c=(T)(-1);              return c;} \
00574     static const T& getTwo()          {static const T c=(T)(2);               return c;} \
00575     static const T& getThree()        {static const T c=(T)(3);               return c;} \
00576     static const T& getOneHalf()      {static const T c=(T)(0.5L);            return c;} \
00577     static const T& getOneThird()     {static const T c=(T)(1.L/3.L);         return c;} \
00578     static const T& getOneFourth()    {static const T c=(T)(0.25L);           return c;} \
00579     static const T& getOneFifth()     {static const T c=(T)(0.2L);            return c;} \
00580     static const T& getOneSixth()     {static const T c=(T)(1.L/6.L);         return c;} \
00581     static const T& getOneSeventh()   {static const T c=(T)(1.L/7.L);         return c;} \
00582     static const T& getOneEighth()    {static const T c=(T)(0.125L);          return c;} \
00583     static const T& getOneNinth()     {static const T c=(T)(1.L/9.L);         return c;} \
00584     static const T& getPi()           {static const T c=(T)(SimTK_PI);        return c;} \
00585     static const T& getOneOverPi()    {static const T c=(T)(1.L/SimTK_PI);    return c;} \
00586     static const T& getE()            {static const T c=(T)(SimTK_E);         return c;} \
00587     static const T& getLog2E()        {static const T c=(T)(SimTK_LOG2E);     return c;} \
00588     static const T& getLog10E()       {static const T c=(T)(SimTK_LOG10E);    return c;} \
00589     static const T& getSqrt2()        {static const T c=(T)(SimTK_SQRT2);     return c;} \
00590     static const T& getOneOverSqrt2() {static const T c=(T)(1.L/SimTK_SQRT2); return c;} \
00591     static const T& getSqrt3()        {static const T c=(T)(SimTK_SQRT3);     return c;} \
00592     static const T& getOneOverSqrt3() {static const T c=(T)(1.L/SimTK_SQRT3); return c;} \
00593     static const T& getCubeRoot2()    {static const T c=(T)(SimTK_CBRT2);     return c;} \
00594     static const T& getCubeRoot3()    {static const T c=(T)(SimTK_CBRT3);     return c;} \
00595     static const T& getLn2()          {static const T c=(T)(SimTK_LN2);       return c;} \
00596     static const T& getLn10()         {static const T c=(T)(SimTK_LN10);      return c;} \
00597     /* integer digit counts useful for formatted input and output */                     \
00598     static const int getNumDigits()         {static const int c=(int)(std::log10(1/getEps()) -0.5); return c;} \
00599     static const int getLosslessNumDigits() {static const int c=(int)(std::log10(1/getTiny())+0.5); return c;} \
00600 }; \
00601 template<> struct NTraits<R>::Result<float> \
00602   {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;};    \
00603 template<> struct NTraits<R>::Result<double> \
00604   {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;};    \
00605 template<> struct NTraits<R>::Result<long double> \
00606   {typedef Widest<R,long double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;};    \
00607 template<> struct NTraits<R>::Result<complex<float> > \
00608   {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
00609 template<> struct NTraits<R>::Result<complex<double> > \
00610   {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
00611 template<> struct NTraits<R>::Result<complex<long double> > \
00612   {typedef Widest<R,complex<long double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
00613 template<> struct NTraits<R>::Result<conjugate<float> > \
00614   {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
00615 template<> struct NTraits<R>::Result<conjugate<double> > \
00616   {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
00617 template<> struct NTraits<R>::Result<conjugate<long double> > \
00618   {typedef conjugate<Widest<R,long double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}
00619 SimTK_DEFINE_REAL_NTRAITS(float);
00620 SimTK_DEFINE_REAL_NTRAITS(double);
00621 SimTK_DEFINE_REAL_NTRAITS(long double);
00622 #undef SimTK_DEFINE_REAL_NTRAITS
00623 
00625 template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { };
00626 template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { };
00627 template <> class CNT<float> : public NTraits<float> { };
00628 template <> class CNT<double> : public NTraits<double> { };
00629 template <> class CNT<long double> : public NTraits<long double> { };
00630 
00631 
00632 } // namespace SimTK
00633 
00634 #endif //SimTK_SIMMATRIX_NTRAITS_H_

Generated on Fri Sep 26 07:44:15 2008 for SimTKcore by  doxygen 1.5.6