Simbody  3.5
Vec.h
Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
00002 #define SimTK_SIMMATRIX_SMALLMATRIX_VEC_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: Peter Eastman                                                *
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 
00031 #include "SimTKcommon/internal/common.h"
00032 
00033 namespace SimTK {
00034 
00035 
00036 // The following functions are used internally by Vec.
00037 
00038 // Hide from Doxygen.
00040 namespace Impl {
00041 
00042 // For those wimpy compilers that don't unroll short, constant-limit loops, 
00043 // Peter Eastman added these recursive template implementations of 
00044 // elementwise add, subtract, and copy. Sherm added multiply and divide.
00045 
00046 template <class E1, int S1, class E2, int S2> void
00047 conformingAdd(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2, 
00048               Vec<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00049     result[0] = r1[0] + r2[0];
00050 }
00051 template <int N, class E1, int S1, class E2, int S2> void
00052 conformingAdd(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2, 
00053               Vec<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00054     conformingAdd(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1), 
00055                   reinterpret_cast<const Vec<N-1,E2,S2>&>(r2), 
00056                   reinterpret_cast<Vec<N-1,typename CNT<E1>::
00057                               template Result<E2>::Add>&>(result));
00058     result[N-1] = r1[N-1] + r2[N-1];
00059 }
00060 
00061 template <class E1, int S1, class E2, int S2> void
00062 conformingSubtract(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2, 
00063                    Vec<1,typename CNT<E1>::template Result<E2>::Sub>& result) {
00064     result[0] = r1[0] - r2[0];
00065 }
00066 template <int N, class E1, int S1, class E2, int S2> void
00067 conformingSubtract(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2, 
00068                    Vec<N,typename CNT<E1>::template Result<E2>::Sub>& result) {
00069     conformingSubtract(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1), 
00070                        reinterpret_cast<const Vec<N-1,E2,S2>&>(r2), 
00071                        reinterpret_cast<Vec<N-1,typename CNT<E1>::
00072                                    template Result<E2>::Sub>&>(result));
00073     result[N-1] = r1[N-1] - r2[N-1];
00074 }
00075 
00076 template <class E1, int S1, class E2, int S2> void
00077 elementwiseMultiply(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2, 
00078               Vec<1,typename CNT<E1>::template Result<E2>::Mul>& result) {
00079     result[0] = r1[0] * r2[0];
00080 }
00081 template <int N, class E1, int S1, class E2, int S2> void
00082 elementwiseMultiply(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2, 
00083               Vec<N,typename CNT<E1>::template Result<E2>::Mul>& result) {
00084     elementwiseMultiply(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1), 
00085                         reinterpret_cast<const Vec<N-1,E2,S2>&>(r2), 
00086                         reinterpret_cast<Vec<N-1,typename CNT<E1>::
00087                                     template Result<E2>::Mul>&>(result));
00088     result[N-1] = r1[N-1] * r2[N-1];
00089 }
00090 
00091 template <class E1, int S1, class E2, int S2> void
00092 elementwiseDivide(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2, 
00093               Vec<1,typename CNT<E1>::template Result<E2>::Dvd>& result) {
00094     result[0] = r1[0] / r2[0];
00095 }
00096 template <int N, class E1, int S1, class E2, int S2> void
00097 elementwiseDivide(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2, 
00098               Vec<N,typename CNT<E1>::template Result<E2>::Dvd>& result) {
00099     elementwiseDivide(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1), 
00100                       reinterpret_cast<const Vec<N-1,E2,S2>&>(r2), 
00101                       reinterpret_cast<Vec<N-1,typename CNT<E1>::
00102                                   template Result<E2>::Dvd>&>(result));
00103     result[N-1] = r1[N-1] / r2[N-1];
00104 }
00105 
00106 template <class E1, int S1, class E2, int S2> void
00107 copy(Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2) {
00108     r1[0] = r2[0];
00109 }
00110 template <int N, class E1, int S1, class E2, int S2> void
00111 copy(Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2) {
00112     copy(reinterpret_cast<Vec<N-1,E1,S1>&>(r1), 
00113          reinterpret_cast<const Vec<N-1,E2,S2>&>(r2));
00114     r1[N-1] = r2[N-1];
00115 }
00116 
00117 }
00183 template <int M, class ELT, int STRIDE>
00184 class Vec {
00185 public:
00191     typedef ELT                                 E;
00193     typedef typename CNT<E>::TNeg               ENeg;
00195     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00198     typedef typename CNT<E>::TReal              EReal;
00202     typedef typename CNT<E>::TImag              EImag;
00205     typedef typename CNT<E>::TComplex           EComplex;
00207     typedef typename CNT<E>::THerm              EHerm;
00209     typedef typename CNT<E>::TPosTrans          EPosTrans;
00212     typedef typename CNT<E>::TSqHermT           ESqHermT;
00214     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00216     typedef typename CNT<E>::TSqrt              ESqrt;
00218     typedef typename CNT<E>::TAbs               EAbs;
00221     typedef typename CNT<E>::TStandard          EStandard;
00224     typedef typename CNT<E>::TInvert            EInvert;
00226     typedef typename CNT<E>::TNormalize         ENormalize;
00227 
00228     typedef typename CNT<E>::Scalar             EScalar;
00229     typedef typename CNT<E>::ULessScalar        EULessScalar;
00230     typedef typename CNT<E>::Number             ENumber;
00231     typedef typename CNT<E>::StdNumber          EStdNumber;
00232     typedef typename CNT<E>::Precision          EPrecision;
00233     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00234 
00237     enum {
00238         NRows               = M,
00239         NCols               = 1,
00240         NPackedElements     = M,
00241         NActualElements     = M * STRIDE,   // includes trailing gap
00242         NActualScalars      = CNT<E>::NActualScalars * NActualElements,
00243         RowSpacing          = STRIDE,
00244         ColSpacing          = NActualElements,
00245         ImagOffset          = NTraits<ENumber>::ImagOffset,
00246         RealStrideFactor    = 1, // composite types don't change size when
00247                                  // cast from complex to real or imaginary
00248         ArgDepth            = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH 
00249                                 ? CNT<E>::ArgDepth + 1 
00250                                 : MAX_RESOLVED_DEPTH),
00251         IsScalar            = 0,
00252         IsULessScalar       = 0,
00253         IsNumber            = 0,
00254         IsStdNumber         = 0,
00255         IsPrecision         = 0,
00256         SignInterpretation  = CNT<E>::SignInterpretation
00257     };
00258 
00259     // These are reinterpretations of the current data, so have the
00260     // same packing (stride).
00261 
00263     typedef Vec<M,E,STRIDE>                 T;
00266     typedef Vec<M,ENeg,STRIDE>              TNeg;
00269     typedef Vec<M,EWithoutNegator,STRIDE>   TWithoutNegator;
00272     typedef Vec<M,EReal,STRIDE*CNT<E>::RealStrideFactor>         
00273                                             TReal;
00276     typedef Vec<M,EImag,STRIDE*CNT<E>::RealStrideFactor>         
00277                                             TImag;
00278     typedef Vec<M,EComplex,STRIDE>          TComplex;
00282     typedef Row<M,EHerm,STRIDE>             THerm;
00285     typedef Row<M,E,STRIDE>                 TPosTrans;
00287     typedef E                               TElement;
00289     typedef E                               TRow;
00291     typedef Vec                             TCol;
00292 
00293     // These are the results of calculations, so are returned in new, packed
00294     // memory. Be sure to refer to element types here which are also packed.
00295     typedef Vec<M,ESqrt,1>                  TSqrt;      // Note stride
00296     typedef Vec<M,EAbs,1>                   TAbs;       // Note stride
00297     typedef Vec<M,EStandard,1>              TStandard;
00298     typedef Row<M,EInvert,1>                TInvert;
00299     typedef Vec<M,ENormalize,1>             TNormalize;
00300 
00301     typedef ESqHermT                        TSqHermT;   // result of self dot product
00302     typedef SymMat<M,ESqTHerm>              TSqTHerm;   // result of self outer product
00303 
00304     // These recurse right down to the underlying scalar type no matter how
00305     // deep the elements are.
00306     typedef EScalar                         Scalar;
00307     typedef EULessScalar                    ULessScalar;
00308     typedef ENumber                         Number;
00309     typedef EStdNumber                      StdNumber;
00310     typedef EPrecision                      Precision;
00311     typedef EScalarNormSq                   ScalarNormSq;
00316     static int size() { return M; }
00318     static int nrow() { return M; }
00320     static int ncol() { return 1; }
00321 
00322 
00325     ScalarNormSq scalarNormSqr() const { 
00326         ScalarNormSq sum(0);
00327         for(int i=0;i<M;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
00328         return sum;
00329     }
00330 
00335     TSqrt sqrt() const {
00336         TSqrt vsqrt;
00337         for(int i=0;i<M;++i) vsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
00338         return vsqrt;
00339     }
00340 
00345     TAbs abs() const {
00346         TAbs vabs;
00347         for(int i=0;i<M;++i) vabs[i] = CNT<E>::abs(d[i*STRIDE]);
00348         return vabs;
00349     }
00350 
00355     TStandard standardize() const {
00356         TStandard vstd;
00357         for(int i=0;i<M;++i) vstd[i] = CNT<E>::standardize(d[i*STRIDE]);
00358         return vstd;
00359     }
00360 
00364     EStandard sum() const {
00365         E sum(0);
00366         for (int i=0;i<M;++i) sum += d[i*STRIDE];
00367         return CNT<E>::standardize(sum);
00368     }
00369 
00370 
00371     // This gives the resulting vector type when (v[i] op P) is applied to 
00372     // each element of v. It is a vector of length M, stride 1, and element 
00373     // types which are the regular composite result of E op P. Typically P is 
00374     // a scalar type but it doesn't have to be.
00375     template <class P> struct EltResult { 
00376         typedef Vec<M, typename CNT<E>::template Result<P>::Mul, 1> Mul;
00377         typedef Vec<M, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
00378         typedef Vec<M, typename CNT<E>::template Result<P>::Add, 1> Add;
00379         typedef Vec<M, typename CNT<E>::template Result<P>::Sub, 1> Sub;
00380     };
00381 
00382     // This is the composite result for v op P where P is some kind of 
00383     // appropriately shaped non-scalar type.
00384     template <class P> struct Result { 
00385         typedef MulCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
00386             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00387             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
00388         typedef typename MulOp::Type Mul;
00389 
00390         typedef MulCNTsNonConforming<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
00391             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00392             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
00393         typedef typename MulOpNonConforming::Type MulNon;
00394 
00395         typedef DvdCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
00396             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00397             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
00398         typedef typename DvdOp::Type Dvd;
00399 
00400         typedef AddCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
00401             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00402             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
00403         typedef typename AddOp::Type Add;
00404 
00405         typedef SubCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
00406             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00407             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
00408         typedef typename SubOp::Type Sub;
00409     };
00410 
00416     template <class P> struct Substitute {
00417         typedef Vec<M,P> Type;
00418     };
00419 
00423     Vec(){ 
00424     #ifndef NDEBUG
00425         setToNaN();
00426     #endif
00427     }
00428 
00429     // It's important not to use the default copy constructor or copy
00430     // assignment because the compiler doesn't understand that we may
00431     // have noncontiguous storage and will try to copy the whole array.
00432 
00436     Vec(const Vec& src) {
00437         Impl::copy(*this, src);
00438     }
00443     Vec& operator=(const Vec& src) {    
00444         Impl::copy(*this, src);
00445         return *this;
00446     }
00447 
00450     template <int SS> Vec(const Vec<M,E,SS>& src) {
00451         Impl::copy(*this, src);
00452     }
00453 
00456     template <int SS> Vec(const Vec<M,ENeg,SS>& src) {
00457         Impl::copy(*this, src);
00458     }
00459 
00462     template <class EE, int SS> explicit Vec(const Vec<M,EE,SS>& src) {
00463         Impl::copy(*this, src);
00464     }
00465 
00468     explicit Vec(const E& e) {for (int i=0;i<M;++i) d[i*STRIDE]=e;}
00469 
00474     explicit Vec(const ENeg& ne) {
00475         const E e = ne; // requires floating point negation
00476         for (int i=0;i<M;++i) d[i*STRIDE]=e;
00477     }
00478 
00483     explicit Vec(int i) {new (this) Vec(E(Precision(i)));}
00484 
00485     // A bevy of constructors for Vecs up to length 9.
00486 
00488     Vec(const E& e0,const E& e1)
00489       { assert(M==2);(*this)[0]=e0;(*this)[1]=e1; }
00490     Vec(const E& e0,const E& e1,const E& e2)
00491       { assert(M==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
00492     Vec(const E& e0,const E& e1,const E& e2,const E& e3)
00493       { assert(M==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
00494     Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
00495       { assert(M==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00496         (*this)[3]=e3;(*this)[4]=e4; }
00497     Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
00498       { assert(M==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00499         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
00500     Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6)
00501       { assert(M==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00502         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
00503     Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6, const E& e7)
00504       { assert(M==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00505         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
00506     Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6, const E& e7, const E& e8)
00507       { assert(M==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00508         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
00509 
00514     template <class EE> explicit Vec(const EE* p)
00515     {   assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; }
00516 
00521     template <class EE> Vec& operator=(const EE* p)
00522     {   assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; return *this; }
00523 
00526     template <class EE, int SS> Vec& operator=(const Vec<M,EE,SS>& vv) 
00527     {   Impl::copy(*this, vv); return *this; }
00528 
00531     template <class EE, int SS> Vec& operator+=(const Vec<M,EE,SS>& r)
00532     {   for(int i=0;i<M;++i) d[i*STRIDE] += r[i]; return *this; }
00536     template <class EE, int SS> Vec& operator+=(const Vec<M,negator<EE>,SS>& r)
00537     {   for(int i=0;i<M;++i) d[i*STRIDE] -= -(r[i]); return *this; }
00538 
00541     template <class EE, int SS> Vec& operator-=(const Vec<M,EE,SS>& r)
00542     {   for(int i=0;i<M;++i) d[i*STRIDE] -= r[i]; return *this; }
00546     template <class EE, int SS> Vec& operator-=(const Vec<M,negator<EE>,SS>& r)
00547     {   for(int i=0;i<M;++i) d[i*STRIDE] += -(r[i]); return *this; }
00548 
00549     // Conforming binary ops with 'this' on left, producing new packed result.
00550     // Cases: v=v+v, v=v-v, m=v*r
00551 
00553     template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Add>
00554     conformingAdd(const Vec<M,EE,SS>& r) const {
00555         Vec<M,typename CNT<E>::template Result<EE>::Add> result;
00556         Impl::conformingAdd(*this, r, result);
00557         return result;
00558     }
00560     template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Sub>
00561     conformingSubtract(const Vec<M,EE,SS>& r) const {
00562         Vec<M,typename CNT<E>::template Result<EE>::Sub> result;
00563         Impl::conformingSubtract(*this, r, result);
00564         return result;
00565     }
00566 
00569     template <class EE, int SS> Mat<M,M,typename CNT<E>::template Result<EE>::Mul>
00570     conformingMultiply(const Row<M,EE,SS>& r) const {
00571         Mat<M,M,typename CNT<E>::template Result<EE>::Mul> result;
00572         for (int j=0;j<M;++j) result(j) = scalarMultiply(r(j));
00573         return result;
00574     }
00575 
00577     template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Mul>
00578     elementwiseMultiply(const Vec<M,EE,SS>& r) const {
00579         Vec<M,typename CNT<E>::template Result<EE>::Mul> result;
00580         Impl::elementwiseMultiply(*this, r, result);
00581         return result;
00582     }
00584     template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Dvd>
00585     elementwiseDivide(const Vec<M,EE,SS>& r) const {
00586         Vec<M,typename CNT<E>::template Result<EE>::Dvd> result;
00587         Impl::elementwiseDivide(*this, r, result);
00588         return result;
00589     }
00590 
00594     const E& operator[](int i) const 
00595     {   assert(0 <= i && i < M); return d[i*STRIDE]; }
00597     const E& operator()(int i) const {return (*this)[i];}
00598 
00602     E& operator[](int i) {assert(0 <= i && i < M); return d[i*STRIDE];}
00604     E& operator()(int i) {return (*this)[i];}
00605 
00606     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00607     typename CNT<ScalarNormSq>::TSqrt 
00608         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00609 
00621     TNormalize normalize() const {
00622         if (CNT<E>::IsScalar) {
00623             return castAwayNegatorIfAny() / (SignInterpretation*norm());
00624         } else {
00625             TNormalize elementwiseNormalized;
00626             for (int i=0; i<M; ++i) 
00627                 elementwiseNormalized[i] = CNT<E>::normalize((*this)[i]);
00628             return elementwiseNormalized;
00629         }
00630     }
00631 
00633     TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
00634 
00636     const Vec&   operator+() const { return *this; }
00640     const TNeg&  operator-() const { return negate(); }
00643     TNeg&        operator-()       { return updNegate(); }
00647     const THerm& operator~() const { return transpose(); }
00651     THerm&       operator~()       { return updTranspose(); }
00652 
00654     const TNeg&  negate() const { return *reinterpret_cast<const TNeg*>(this); }
00657     TNeg&        updNegate()    { return *reinterpret_cast<      TNeg*>(this); }
00658 
00660     const THerm& transpose()    const { return *reinterpret_cast<const THerm*>(this); }
00663     THerm&       updTranspose()       { return *reinterpret_cast<      THerm*>(this); }
00664 
00669     const TPosTrans& positionalTranspose() const
00670         { return *reinterpret_cast<const TPosTrans*>(this); }
00672     TPosTrans&       updPositionalTranspose()
00673         { return *reinterpret_cast<TPosTrans*>(this); }
00674 
00679     const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
00682     TReal&       real()       { return *reinterpret_cast<      TReal*>(this); }
00683 
00684     // Had to contort these next two routines to get them through VC++ 7.net
00685 
00690     const TImag& imag()    const { 
00691         const int offs = ImagOffset;
00692         const EImag* p = reinterpret_cast<const EImag*>(this);
00693         return *reinterpret_cast<const TImag*>(p+offs);
00694     }
00697     TImag& imag() { 
00698         const int offs = ImagOffset;
00699         EImag* p = reinterpret_cast<EImag*>(this);
00700         return *reinterpret_cast<TImag*>(p+offs);
00701     }
00702 
00706     const TWithoutNegator& castAwayNegatorIfAny() const 
00707     {   return *reinterpret_cast<const TWithoutNegator*>(this); }
00710     TWithoutNegator&       updCastAwayNegatorIfAny()    
00711     {   return *reinterpret_cast<TWithoutNegator*>(this); }
00712 
00713     // These are elementwise binary operators, (this op ee) by default but 
00714     // (ee op this) if 'FromLeft' appears in the name. The result is a packed 
00715     // Vec<M> but the element type may change. These are mostly used to 
00716     // implement global operators. We call these "scalar" operators but 
00717     // actually the "scalar" can be a composite type.
00718 
00719     //TODO: consider converting 'e' to Standard Numbers as precalculation and 
00720     // changing return type appropriately.
00721     template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Mul>
00722     scalarMultiply(const EE& e) const {
00723         Vec<M, typename CNT<E>::template Result<EE>::Mul> result;
00724         for (int i=0; i<M; ++i) result[i] = (*this)[i] * e;
00725         return result;
00726     }
00727     template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Mul>
00728     scalarMultiplyFromLeft(const EE& e) const {
00729         Vec<M, typename CNT<EE>::template Result<E>::Mul> result;
00730         for (int i=0; i<M; ++i) result[i] = e * (*this)[i];
00731         return result;
00732     }
00733 
00734     // TODO: should precalculate and store 1/e, while converting to Standard 
00735     // Numbers. Note that return type should change appropriately.
00736     template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Dvd>
00737     scalarDivide(const EE& e) const {
00738         Vec<M, typename CNT<E>::template Result<EE>::Dvd> result;
00739         for (int i=0; i<M; ++i) result[i] = (*this)[i] / e;
00740         return result;
00741     }
00742     template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Dvd>
00743     scalarDivideFromLeft(const EE& e) const {
00744         Vec<M, typename CNT<EE>::template Result<E>::Dvd> result;
00745         for (int i=0; i<M; ++i) result[i] = e / (*this)[i];
00746         return result;
00747     }
00748 
00749     template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Add>
00750     scalarAdd(const EE& e) const {
00751         Vec<M, typename CNT<E>::template Result<EE>::Add> result;
00752         for (int i=0; i<M; ++i) result[i] = (*this)[i] + e;
00753         return result;
00754     }
00755     // Add is commutative, so no 'FromLeft'.
00756 
00757     template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Sub>
00758     scalarSubtract(const EE& e) const {
00759         Vec<M, typename CNT<E>::template Result<EE>::Sub> result;
00760         for (int i=0; i<M; ++i) result[i] = (*this)[i] - e;
00761         return result;
00762     }
00763     template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Sub>
00764     scalarSubtractFromLeft(const EE& e) const {
00765         Vec<M, typename CNT<EE>::template Result<E>::Sub> result;
00766         for (int i=0; i<M; ++i) result[i] = e - (*this)[i];
00767         return result;
00768     }
00769 
00770     // Generic assignments for any element type not listed explicitly, including scalars.
00771     // These are done repeatedly for each element and only work if the operation can
00772     // be performed leaving the original element type.
00773     template <class EE> Vec& operator =(const EE& e) {return scalarEq(e);}
00774     template <class EE> Vec& operator+=(const EE& e) {return scalarPlusEq(e);}
00775     template <class EE> Vec& operator-=(const EE& e) {return scalarMinusEq(e);}
00776     template <class EE> Vec& operator*=(const EE& e) {return scalarTimesEq(e);}
00777     template <class EE> Vec& operator/=(const EE& e) {return scalarDivideEq(e);}
00778 
00779     // Generalized element assignment & computed assignment methods. These will work
00780     // for any assignment-compatible element, not just scalars.
00781     template <class EE> Vec& scalarEq(const EE& ee)
00782       { for(int i=0;i<M;++i) d[i*STRIDE] = ee; return *this; }
00783     template <class EE> Vec& scalarPlusEq(const EE& ee)
00784       { for(int i=0;i<M;++i) d[i*STRIDE] += ee; return *this; }
00785     template <class EE> Vec& scalarMinusEq(const EE& ee)
00786       { for(int i=0;i<M;++i) d[i*STRIDE] -= ee; return *this; }
00787     template <class EE> Vec& scalarMinusEqFromLeft(const EE& ee)
00788       { for(int i=0;i<M;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
00789     template <class EE> Vec& scalarTimesEq(const EE& ee)
00790       { for(int i=0;i<M;++i) d[i*STRIDE] *= ee; return *this; }
00791     template <class EE> Vec& scalarTimesEqFromLeft(const EE& ee)
00792       { for(int i=0;i<M;++i) d[i*STRIDE] = ee * d[i*STRIDE]; return *this; }
00793     template <class EE> Vec& scalarDivideEq(const EE& ee)
00794       { for(int i=0;i<M;++i) d[i*STRIDE] /= ee; return *this; }
00795     template <class EE> Vec& scalarDivideEqFromLeft(const EE& ee)
00796       { for(int i=0;i<M;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
00797 
00798     // Specialize for int to avoid warnings and ambiguities.
00799     Vec& scalarEq(int ee)       {return scalarEq(Precision(ee));}
00800     Vec& scalarPlusEq(int ee)   {return scalarPlusEq(Precision(ee));}
00801     Vec& scalarMinusEq(int ee)  {return scalarMinusEq(Precision(ee));}
00802     Vec& scalarTimesEq(int ee)  {return scalarTimesEq(Precision(ee));}
00803     Vec& scalarDivideEq(int ee) {return scalarDivideEq(Precision(ee));}
00804     Vec& scalarMinusEqFromLeft(int ee)  {return scalarMinusEqFromLeft(Precision(ee));}
00805     Vec& scalarTimesEqFromLeft(int ee)  {return scalarTimesEqFromLeft(Precision(ee));}
00806     Vec& scalarDivideEqFromLeft(int ee) {return scalarDivideEqFromLeft(Precision(ee));}
00807 
00810     void setToNaN() {
00811         (*this) = CNT<ELT>::getNaN();
00812     }
00813 
00815     void setToZero() {
00816         (*this) = ELT(0);
00817     }
00818 
00824     template <int MM>
00825     const Vec<MM,ELT,STRIDE>& getSubVec(int i) const {
00826         assert(0 <= i && i + MM <= M);
00827         return Vec<MM,ELT,STRIDE>::getAs(&(*this)[i]);
00828     }
00834     template <int MM>
00835     Vec<MM,ELT,STRIDE>& updSubVec(int i) {
00836         assert(0 <= i && i + MM <= M);
00837         return Vec<MM,ELT,STRIDE>::updAs(&(*this)[i]);
00838     }
00839 
00840 
00844     template <int MM>
00845     static const Vec& getSubVec(const Vec<MM,ELT,STRIDE>& v, int i) {
00846         assert(0 <= i && i + M <= MM);
00847         return getAs(&v[i]);
00848     }
00852     template <int MM>
00853     static Vec& updSubVec(Vec<MM,ELT,STRIDE>& v, int i) {
00854         assert(0 <= i && i + M <= MM);
00855         return updAs(&v[i]);
00856     }
00857 
00861     Vec<M-1,ELT,1> drop1(int p) const {
00862         assert(0 <= p && p < M);
00863         Vec<M-1,ELT,1> out;
00864         int nxt=0;
00865         for (int i=0; i<M-1; ++i, ++nxt) {
00866             if (nxt==p) ++nxt;  // skip the loser
00867             out[i] = (*this)[nxt];
00868         }
00869         return out;
00870     }
00871 
00875     template <class EE> Vec<M+1,ELT,1> append1(const EE& v) const {
00876         Vec<M+1,ELT,1> out;
00877         Vec<M,ELT,1>::updAs(&out[0]) = (*this);
00878         out[M] = v;
00879         return out;
00880     }
00881 
00882 
00888     template <class EE> Vec<M+1,ELT,1> insert1(int p, const EE& v) const {
00889         assert(0 <= p && p <= M);
00890         if (p==M) return append1(v);
00891         Vec<M+1,ELT,1> out;
00892         int nxt=0;
00893         for (int i=0; i<M; ++i, ++nxt) {
00894             if (i==p) out[nxt++] = v;
00895             out[nxt] = (*this)[i];
00896         }
00897         return out;
00898     }
00899             
00902     static const Vec& getAs(const ELT* p)  
00903     {   return *reinterpret_cast<const Vec*>(p); }
00906     static Vec&       updAs(ELT* p)
00907     {   return *reinterpret_cast<Vec*>(p); }
00908 
00909 
00913     static Vec<M,ELT,1> getNaN() { return Vec<M,ELT,1>(CNT<ELT>::getNaN()); }
00914 
00916     bool isNaN() const {
00917         for (int i=0; i<M; ++i)
00918             if (CNT<ELT>::isNaN((*this)[i]))
00919                 return true;
00920         return false;
00921     }
00922 
00925     bool isInf() const {
00926         bool seenInf = false;
00927         for (int i=0; i<M; ++i) {
00928             const ELT& e = (*this)[i];
00929             if (!CNT<ELT>::isFinite(e)) {
00930                 if (!CNT<ELT>::isInf(e)) 
00931                     return false; // something bad was found
00932                 seenInf = true; 
00933             }
00934         }
00935         return seenInf;
00936     }
00937 
00940     bool isFinite() const {
00941         for (int i=0; i<M; ++i)
00942             if (!CNT<ELT>::isFinite((*this)[i]))
00943                 return false;
00944         return true;
00945     }
00946 
00949     static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
00950 
00953     template <class E2, int RS2>
00954     bool isNumericallyEqual(const Vec<M,E2,RS2>& v, double tol) const {
00955         for (int i=0; i<M; ++i)
00956             if (!CNT<ELT>::isNumericallyEqual((*this)[i], v[i], tol))
00957                 return false;
00958         return true;
00959     }
00960 
00964     template <class E2, int RS2>
00965     bool isNumericallyEqual(const Vec<M,E2,RS2>& v) const {
00966         const double tol = std::max(getDefaultTolerance(),v.getDefaultTolerance());
00967         return isNumericallyEqual(v, tol);
00968     }
00969 
00974     bool isNumericallyEqual
00975        (const ELT& e,
00976         double     tol = getDefaultTolerance()) const 
00977     {
00978         for (int i=0; i<M; ++i)
00979             if (!CNT<ELT>::isNumericallyEqual((*this)[i], e, tol))
00980                 return false;
00981         return true;
00982     }
00983 
00984     // Functions to be used for Scripting in MATLAB and languages that do not support operator overloading
00986     std::string toString() const {
00987         std::stringstream stream;
00988         stream <<  (*this);
00989         return stream.str(); 
00990     }
00991 
00993     void set(int i, const E& value)  
00994     {   (*this)[i] = value; }
00995 
00997     const E& get(int i) const 
00998     {   return operator[](i); }
00999 
01000 private:
01001     // TODO: should be an array of scalars rather than elements to control
01002     // packing more carefully.
01003     ELT d[NActualElements];    // data
01004 };
01005 
01007 // Global operators involving two vectors. //
01008 //   v+v, v-v, v==v, v!=v                  //
01010 
01011 // v3 = v1 + v2 where all v's have the same length M. 
01012 template <int M, class E1, int S1, class E2, int S2> inline
01013 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Add
01014 operator+(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {
01015     return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
01016         ::AddOp::perform(l,r);
01017 }
01018 
01019 // v3 = v1 - v2, similar to +
01020 template <int M, class E1, int S1, class E2, int S2> inline
01021 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Sub
01022 operator-(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) { 
01023     return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
01024         ::SubOp::perform(l,r);
01025 }
01026 
01028 template <int M, class E1, int S1, class E2, int S2> inline bool
01029 operator==(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) 
01030 {   for (int i=0; i < M; ++i) if (l[i] != r[i]) return false;
01031     return true; }
01033 template <int M, class E1, int S1, class E2, int S2> inline bool
01034 operator!=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {return !(l==r);} 
01035 
01037 template <int M, class E1, int S1, class E2> inline bool
01038 operator==(const Vec<M,E1,S1>& v, const E2& e) 
01039 {   for (int i=0; i < M; ++i) if (v[i] != e) return false;
01040     return true; }
01042 template <int M, class E1, int S1, class E2> inline bool
01043 operator!=(const Vec<M,E1,S1>& v, const E2& e) {return !(v==e);} 
01044 
01046 template <int M, class E1, int S1, class E2, int S2> inline bool
01047 operator<(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) 
01048 {   for (int i=0; i < M; ++i) if (l[i] >= r[i]) return false;
01049     return true; }
01051 template <int M, class E1, int S1, class E2> inline bool
01052 operator<(const Vec<M,E1,S1>& v, const E2& e) 
01053 {   for (int i=0; i < M; ++i) if (v[i] >= e) return false;
01054     return true; }
01055 
01057 template <int M, class E1, int S1, class E2, int S2> inline bool
01058 operator>(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) 
01059 {   for (int i=0; i < M; ++i) if (l[i] <= r[i]) return false;
01060     return true; }
01062 template <int M, class E1, int S1, class E2> inline bool
01063 operator>(const Vec<M,E1,S1>& v, const E2& e) 
01064 {   for (int i=0; i < M; ++i) if (v[i] <= e) return false;
01065     return true; }
01066 
01069 template <int M, class E1, int S1, class E2, int S2> inline bool
01070 operator<=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) 
01071 {   for (int i=0; i < M; ++i) if (l[i] > r[i]) return false;
01072     return true; }
01075 template <int M, class E1, int S1, class E2> inline bool
01076 operator<=(const Vec<M,E1,S1>& v, const E2& e) 
01077 {   for (int i=0; i < M; ++i) if (v[i] > e) return false;
01078     return true; }
01079 
01082 template <int M, class E1, int S1, class E2, int S2> inline bool
01083 operator>=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) 
01084 {   for (int i=0; i < M; ++i) if (l[i] < r[i]) return false;
01085     return true; }
01088 template <int M, class E1, int S1, class E2> inline bool
01089 operator>=(const Vec<M,E1,S1>& v, const E2& e) 
01090 {   for (int i=0; i < M; ++i) if (v[i] < e) return false;
01091     return true; }
01092 
01094 // Global operators involving a vector and a scalar. //
01096 
01097 // I haven't been able to figure out a nice way to templatize for the
01098 // built-in reals without introducing a lot of unwanted type matches
01099 // as well. So we'll just grind them out explicitly here.
01100 
01101 // SCALAR MULTIPLY
01102 
01103 // v = v*real, real*v 
01104 template <int M, class E, int S> inline
01105 typename Vec<M,E,S>::template Result<float>::Mul
01106 operator*(const Vec<M,E,S>& l, const float& r)
01107   { return Vec<M,E,S>::template Result<float>::MulOp::perform(l,r); }
01108 template <int M, class E, int S> inline
01109 typename Vec<M,E,S>::template Result<float>::Mul
01110 operator*(const float& l, const Vec<M,E,S>& r) {return r*l;}
01111 
01112 template <int M, class E, int S> inline
01113 typename Vec<M,E,S>::template Result<double>::Mul
01114 operator*(const Vec<M,E,S>& l, const double& r)
01115   { return Vec<M,E,S>::template Result<double>::MulOp::perform(l,r); }
01116 template <int M, class E, int S> inline
01117 typename Vec<M,E,S>::template Result<double>::Mul
01118 operator*(const double& l, const Vec<M,E,S>& r) {return r*l;}
01119 
01120 template <int M, class E, int S> inline
01121 typename Vec<M,E,S>::template Result<long double>::Mul
01122 operator*(const Vec<M,E,S>& l, const long double& r)
01123   { return Vec<M,E,S>::template Result<long double>::MulOp::perform(l,r); }
01124 template <int M, class E, int S> inline
01125 typename Vec<M,E,S>::template Result<long double>::Mul
01126 operator*(const long double& l, const Vec<M,E,S>& r) {return r*l;}
01127 
01128 // v = v*int, int*v -- just convert int to v's precision float
01129 template <int M, class E, int S> inline
01130 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
01131 operator*(const Vec<M,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
01132 template <int M, class E, int S> inline
01133 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
01134 operator*(int l, const Vec<M,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
01135 
01136 // Complex, conjugate, and negator are all easy to templatize.
01137 
01138 // v = v*complex, complex*v
01139 template <int M, class E, int S, class R> inline
01140 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
01141 operator*(const Vec<M,E,S>& l, const std::complex<R>& r)
01142   { return Vec<M,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
01143 template <int M, class E, int S, class R> inline
01144 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
01145 operator*(const std::complex<R>& l, const Vec<M,E,S>& r) {return r*l;}
01146 
01147 // v = v*conjugate, conjugate*v (convert conjugate->complex)
01148 template <int M, class E, int S, class R> inline
01149 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
01150 operator*(const Vec<M,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
01151 template <int M, class E, int S, class R> inline
01152 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
01153 operator*(const conjugate<R>& l, const Vec<M,E,S>& r) {return r*(std::complex<R>)l;}
01154 
01155 // v = v*negator, negator*v: convert negator to standard number
01156 template <int M, class E, int S, class R> inline
01157 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
01158 operator*(const Vec<M,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
01159 template <int M, class E, int S, class R> inline
01160 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
01161 operator*(const negator<R>& l, const Vec<M,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
01162 
01163 
01164 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
01165 // but when it is on the left it means scalar * pseudoInverse(vec), which is 
01166 // a row.
01167 
01168 // v = v/real, real/v 
01169 template <int M, class E, int S> inline
01170 typename Vec<M,E,S>::template Result<float>::Dvd
01171 operator/(const Vec<M,E,S>& l, const float& r)
01172   { return Vec<M,E,S>::template Result<float>::DvdOp::perform(l,r); }
01173 template <int M, class E, int S> inline
01174 typename CNT<float>::template Result<Vec<M,E,S> >::Dvd
01175 operator/(const float& l, const Vec<M,E,S>& r)
01176   { return CNT<float>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
01177 
01178 template <int M, class E, int S> inline
01179 typename Vec<M,E,S>::template Result<double>::Dvd
01180 operator/(const Vec<M,E,S>& l, const double& r)
01181   { return Vec<M,E,S>::template Result<double>::DvdOp::perform(l,r); }
01182 template <int M, class E, int S> inline
01183 typename CNT<double>::template Result<Vec<M,E,S> >::Dvd
01184 operator/(const double& l, const Vec<M,E,S>& r)
01185   { return CNT<double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
01186 
01187 template <int M, class E, int S> inline
01188 typename Vec<M,E,S>::template Result<long double>::Dvd
01189 operator/(const Vec<M,E,S>& l, const long double& r)
01190   { return Vec<M,E,S>::template Result<long double>::DvdOp::perform(l,r); }
01191 template <int M, class E, int S> inline
01192 typename CNT<long double>::template Result<Vec<M,E,S> >::Dvd
01193 operator/(const long double& l, const Vec<M,E,S>& r)
01194   { return CNT<long double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
01195 
01196 // v = v/int, int/v -- just convert int to v's precision float
01197 template <int M, class E, int S> inline
01198 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Dvd
01199 operator/(const Vec<M,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
01200 template <int M, class E, int S> inline
01201 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Dvd
01202 operator/(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
01203 
01204 
01205 // Complex, conjugate, and negator are all easy to templatize.
01206 
01207 // v = v/complex, complex/v
01208 template <int M, class E, int S, class R> inline
01209 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
01210 operator/(const Vec<M,E,S>& l, const std::complex<R>& r)
01211   { return Vec<M,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
01212 template <int M, class E, int S, class R> inline
01213 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
01214 operator/(const std::complex<R>& l, const Vec<M,E,S>& r)
01215   { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
01216 
01217 // v = v/conjugate, conjugate/v (convert conjugate->complex)
01218 template <int M, class E, int S, class R> inline
01219 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
01220 operator/(const Vec<M,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
01221 template <int M, class E, int S, class R> inline
01222 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
01223 operator/(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l/r;}
01224 
01225 // v = v/negator, negator/v: convert negator to number
01226 template <int M, class E, int S, class R> inline
01227 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
01228 operator/(const Vec<M,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
01229 template <int M, class E, int S, class R> inline
01230 typename CNT<R>::template Result<Vec<M,E,S> >::Dvd
01231 operator/(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
01232 
01233 
01234 // Add and subtract are odd as scalar ops. They behave as though the
01235 // scalar stands for a vector each of whose elements is that scalar,
01236 // and then a normal vector add or subtract is done.
01237 
01238 // SCALAR ADD
01239 
01240 // v = v+real, real+v 
01241 template <int M, class E, int S> inline
01242 typename Vec<M,E,S>::template Result<float>::Add
01243 operator+(const Vec<M,E,S>& l, const float& r)
01244   { return Vec<M,E,S>::template Result<float>::AddOp::perform(l,r); }
01245 template <int M, class E, int S> inline
01246 typename Vec<M,E,S>::template Result<float>::Add
01247 operator+(const float& l, const Vec<M,E,S>& r) {return r+l;}
01248 
01249 template <int M, class E, int S> inline
01250 typename Vec<M,E,S>::template Result<double>::Add
01251 operator+(const Vec<M,E,S>& l, const double& r)
01252   { return Vec<M,E,S>::template Result<double>::AddOp::perform(l,r); }
01253 template <int M, class E, int S> inline
01254 typename Vec<M,E,S>::template Result<double>::Add
01255 operator+(const double& l, const Vec<M,E,S>& r) {return r+l;}
01256 
01257 template <int M, class E, int S> inline
01258 typename Vec<M,E,S>::template Result<long double>::Add
01259 operator+(const Vec<M,E,S>& l, const long double& r)
01260   { return Vec<M,E,S>::template Result<long double>::AddOp::perform(l,r); }
01261 template <int M, class E, int S> inline
01262 typename Vec<M,E,S>::template Result<long double>::Add
01263 operator+(const long double& l, const Vec<M,E,S>& r) {return r+l;}
01264 
01265 // v = v+int, int+v -- just convert int to v's precision float
01266 template <int M, class E, int S> inline
01267 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
01268 operator+(const Vec<M,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
01269 template <int M, class E, int S> inline
01270 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
01271 operator+(int l, const Vec<M,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
01272 
01273 // Complex, conjugate, and negator are all easy to templatize.
01274 
01275 // v = v+complex, complex+v
01276 template <int M, class E, int S, class R> inline
01277 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
01278 operator+(const Vec<M,E,S>& l, const std::complex<R>& r)
01279   { return Vec<M,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
01280 template <int M, class E, int S, class R> inline
01281 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
01282 operator+(const std::complex<R>& l, const Vec<M,E,S>& r) {return r+l;}
01283 
01284 // v = v+conjugate, conjugate+v (convert conjugate->complex)
01285 template <int M, class E, int S, class R> inline
01286 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
01287 operator+(const Vec<M,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
01288 template <int M, class E, int S, class R> inline
01289 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
01290 operator+(const conjugate<R>& l, const Vec<M,E,S>& r) {return r+(std::complex<R>)l;}
01291 
01292 // v = v+negator, negator+v: convert negator to standard number
01293 template <int M, class E, int S, class R> inline
01294 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
01295 operator+(const Vec<M,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
01296 template <int M, class E, int S, class R> inline
01297 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
01298 operator+(const negator<R>& l, const Vec<M,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
01299 
01300 // SCALAR SUBTRACT -- careful, not commutative.
01301 
01302 // v = v-real, real-v 
01303 template <int M, class E, int S> inline
01304 typename Vec<M,E,S>::template Result<float>::Sub
01305 operator-(const Vec<M,E,S>& l, const float& r)
01306   { return Vec<M,E,S>::template Result<float>::SubOp::perform(l,r); }
01307 template <int M, class E, int S> inline
01308 typename CNT<float>::template Result<Vec<M,E,S> >::Sub
01309 operator-(const float& l, const Vec<M,E,S>& r)
01310   { return CNT<float>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
01311 
01312 template <int M, class E, int S> inline
01313 typename Vec<M,E,S>::template Result<double>::Sub
01314 operator-(const Vec<M,E,S>& l, const double& r)
01315   { return Vec<M,E,S>::template Result<double>::SubOp::perform(l,r); }
01316 template <int M, class E, int S> inline
01317 typename CNT<double>::template Result<Vec<M,E,S> >::Sub
01318 operator-(const double& l, const Vec<M,E,S>& r)
01319   { return CNT<double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
01320 
01321 template <int M, class E, int S> inline
01322 typename Vec<M,E,S>::template Result<long double>::Sub
01323 operator-(const Vec<M,E,S>& l, const long double& r)
01324   { return Vec<M,E,S>::template Result<long double>::SubOp::perform(l,r); }
01325 template <int M, class E, int S> inline
01326 typename CNT<long double>::template Result<Vec<M,E,S> >::Sub
01327 operator-(const long double& l, const Vec<M,E,S>& r)
01328   { return CNT<long double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
01329 
01330 // v = v-int, int-v // just convert int to v's precision float
01331 template <int M, class E, int S> inline
01332 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Sub
01333 operator-(const Vec<M,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
01334 template <int M, class E, int S> inline
01335 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Sub
01336 operator-(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
01337 
01338 
01339 // Complex, conjugate, and negator are all easy to templatize.
01340 
01341 // v = v-complex, complex-v
01342 template <int M, class E, int S, class R> inline
01343 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
01344 operator-(const Vec<M,E,S>& l, const std::complex<R>& r)
01345   { return Vec<M,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
01346 template <int M, class E, int S, class R> inline
01347 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
01348 operator-(const std::complex<R>& l, const Vec<M,E,S>& r)
01349   { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
01350 
01351 // v = v-conjugate, conjugate-v (convert conjugate->complex)
01352 template <int M, class E, int S, class R> inline
01353 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
01354 operator-(const Vec<M,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
01355 template <int M, class E, int S, class R> inline
01356 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
01357 operator-(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l-r;}
01358 
01359 // v = v-negator, negator-v: convert negator to standard number
01360 template <int M, class E, int S, class R> inline
01361 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Sub
01362 operator-(const Vec<M,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
01363 template <int M, class E, int S, class R> inline
01364 typename CNT<R>::template Result<Vec<M,E,S> >::Sub
01365 operator-(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
01366 
01367 // Vec I/O
01368 template <int M, class E, int S, class CHAR, class TRAITS> inline
01369 std::basic_ostream<CHAR,TRAITS>&
01370 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Vec<M,E,S>& v) {
01371     o << "~[" << v[0]; for(int i=1;i<M;++i) o<<','<<v[i]; o<<']'; return o;
01372 }
01373 
01376 template <int M, class E, int S, class CHAR, class TRAITS> inline
01377 std::basic_istream<CHAR,TRAITS>&
01378 operator>>(std::basic_istream<CHAR,TRAITS>& is, Vec<M,E,S>& v) {
01379     CHAR tilde;
01380     is >> tilde; if (is.fail()) return is;
01381     if (tilde != CHAR('~')) {
01382         tilde = CHAR(0);
01383         is.unget(); if (is.fail()) return is;
01384     }
01385 
01386     CHAR openBracket, closeBracket;
01387     is >> openBracket; if (is.fail()) return is;
01388     if (openBracket==CHAR('('))
01389         closeBracket = CHAR(')');
01390     else if (openBracket==CHAR('['))
01391         closeBracket = CHAR(']');
01392     else {
01393         closeBracket = CHAR(0);
01394         is.unget(); if (is.fail()) return is;
01395     }
01396 
01397     // If we saw a "~" but then we didn't see any brackets, that's an
01398     // error. Set the fail bit and return.
01399     if (tilde != CHAR(0) && closeBracket == CHAR(0)) {
01400         is.setstate( std::ios::failbit );
01401         return is;
01402     }
01403 
01404     for (int i=0; i < M; ++i) {
01405         is >> v[i];
01406         if (is.fail()) return is;
01407         if (i != M-1) {
01408             CHAR c; is >> c; if (is.fail()) return is;
01409             if (c != ',') is.unget();
01410             if (is.fail()) return is;
01411         }
01412     }
01413 
01414     // Get the closing bracket if there was an opening one. If we don't
01415     // see the expected character we'll set the fail bit in the istream.
01416     if (closeBracket != CHAR(0)) {
01417         CHAR closer; is >> closer; if (is.fail()) return is;
01418         if (closer != closeBracket) {
01419             is.unget(); if (is.fail()) return is;
01420             is.setstate( std::ios::failbit );
01421         }
01422     }
01423 
01424     return is;
01425 }
01426 
01427 } //namespace SimTK
01428 
01429 
01430 #endif //SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines