Simbody
3.5
|
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_