Simbody  3.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Vec.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
2 #define SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: Peter Eastman *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
31 #include "SimTKcommon/internal/common.h"
32 
33 namespace SimTK {
34 
35 
36 // The following functions are used internally by Vec.
37 
38 // Hide from Doxygen.
40 namespace Impl {
41 
42 // For those wimpy compilers that don't unroll short, constant-limit loops,
43 // Peter Eastman added these recursive template implementations of
44 // elementwise add, subtract, and copy. Sherm added multiply and divide.
45 
46 template <class E1, int S1, class E2, int S2> void
47 conformingAdd(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
48  Vec<1,typename CNT<E1>::template Result<E2>::Add>& result) {
49  result[0] = r1[0] + r2[0];
50 }
51 template <int N, class E1, int S1, class E2, int S2> void
52 conformingAdd(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
53  Vec<N,typename CNT<E1>::template Result<E2>::Add>& result) {
54  conformingAdd(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
55  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
56  reinterpret_cast<Vec<N-1,typename CNT<E1>::
57  template Result<E2>::Add>&>(result));
58  result[N-1] = r1[N-1] + r2[N-1];
59 }
60 
61 template <class E1, int S1, class E2, int S2> void
62 conformingSubtract(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
63  Vec<1,typename CNT<E1>::template Result<E2>::Sub>& result) {
64  result[0] = r1[0] - r2[0];
65 }
66 template <int N, class E1, int S1, class E2, int S2> void
67 conformingSubtract(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
68  Vec<N,typename CNT<E1>::template Result<E2>::Sub>& result) {
69  conformingSubtract(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
70  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
71  reinterpret_cast<Vec<N-1,typename CNT<E1>::
72  template Result<E2>::Sub>&>(result));
73  result[N-1] = r1[N-1] - r2[N-1];
74 }
75 
76 template <class E1, int S1, class E2, int S2> void
77 elementwiseMultiply(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
78  Vec<1,typename CNT<E1>::template Result<E2>::Mul>& result) {
79  result[0] = r1[0] * r2[0];
80 }
81 template <int N, class E1, int S1, class E2, int S2> void
82 elementwiseMultiply(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
83  Vec<N,typename CNT<E1>::template Result<E2>::Mul>& result) {
84  elementwiseMultiply(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
85  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
86  reinterpret_cast<Vec<N-1,typename CNT<E1>::
87  template Result<E2>::Mul>&>(result));
88  result[N-1] = r1[N-1] * r2[N-1];
89 }
90 
91 template <class E1, int S1, class E2, int S2> void
92 elementwiseDivide(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
93  Vec<1,typename CNT<E1>::template Result<E2>::Dvd>& result) {
94  result[0] = r1[0] / r2[0];
95 }
96 template <int N, class E1, int S1, class E2, int S2> void
97 elementwiseDivide(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
98  Vec<N,typename CNT<E1>::template Result<E2>::Dvd>& result) {
99  elementwiseDivide(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
100  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
101  reinterpret_cast<Vec<N-1,typename CNT<E1>::
102  template Result<E2>::Dvd>&>(result));
103  result[N-1] = r1[N-1] / r2[N-1];
104 }
105 
106 template <class E1, int S1, class E2, int S2> void
107 copy(Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2) {
108  r1[0] = r2[0];
109 }
110 template <int N, class E1, int S1, class E2, int S2> void
111 copy(Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2) {
112  copy(reinterpret_cast<Vec<N-1,E1,S1>&>(r1),
113  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2));
114  r1[N-1] = r2[N-1];
115 }
116 
117 }
183 template <int M, class ELT, int STRIDE>
184 class Vec {
185 public:
191  typedef ELT E;
193  typedef typename CNT<E>::TNeg ENeg;
195  typedef typename CNT<E>::TWithoutNegator EWithoutNegator;
198  typedef typename CNT<E>::TReal EReal;
202  typedef typename CNT<E>::TImag EImag;
205  typedef typename CNT<E>::TComplex EComplex;
207  typedef typename CNT<E>::THerm EHerm;
209  typedef typename CNT<E>::TPosTrans EPosTrans;
212  typedef typename CNT<E>::TSqHermT ESqHermT;
214  typedef typename CNT<E>::TSqTHerm ESqTHerm;
216  typedef typename CNT<E>::TSqrt ESqrt;
218  typedef typename CNT<E>::TAbs EAbs;
221  typedef typename CNT<E>::TStandard EStandard;
224  typedef typename CNT<E>::TInvert EInvert;
226  typedef typename CNT<E>::TNormalize ENormalize;
227 
228  typedef typename CNT<E>::Scalar EScalar;
229  typedef typename CNT<E>::ULessScalar EULessScalar;
230  typedef typename CNT<E>::Number ENumber;
231  typedef typename CNT<E>::StdNumber EStdNumber;
232  typedef typename CNT<E>::Precision EPrecision;
233  typedef typename CNT<E>::ScalarNormSq EScalarNormSq;
234 
237  enum {
238  NRows = M,
239  NCols = 1,
240  NPackedElements = M,
241  NActualElements = M * STRIDE, // includes trailing gap
242  NActualScalars = CNT<E>::NActualScalars * NActualElements,
243  RowSpacing = STRIDE,
244  ColSpacing = NActualElements,
245  ImagOffset = NTraits<ENumber>::ImagOffset,
246  RealStrideFactor = 1, // composite types don't change size when
247  // cast from complex to real or imaginary
248  ArgDepth = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH
249  ? CNT<E>::ArgDepth + 1
250  : MAX_RESOLVED_DEPTH),
251  IsScalar = 0,
252  IsULessScalar = 0,
253  IsNumber = 0,
254  IsStdNumber = 0,
255  IsPrecision = 0,
256  SignInterpretation = CNT<E>::SignInterpretation
257  };
258 
259  // These are reinterpretations of the current data, so have the
260  // same packing (stride).
261 
263  typedef Vec<M,E,STRIDE> T;
266  typedef Vec<M,ENeg,STRIDE> TNeg;
269  typedef Vec<M,EWithoutNegator,STRIDE> TWithoutNegator;
272  typedef Vec<M,EReal,STRIDE*CNT<E>::RealStrideFactor>
273  TReal;
276  typedef Vec<M,EImag,STRIDE*CNT<E>::RealStrideFactor>
277  TImag;
278  typedef Vec<M,EComplex,STRIDE> TComplex;
282  typedef Row<M,EHerm,STRIDE> THerm;
285  typedef Row<M,E,STRIDE> TPosTrans;
287  typedef E TElement;
289  typedef E TRow;
291  typedef Vec TCol;
292 
293  // These are the results of calculations, so are returned in new, packed
294  // memory. Be sure to refer to element types here which are also packed.
295  typedef Vec<M,ESqrt,1> TSqrt; // Note stride
296  typedef Vec<M,EAbs,1> TAbs; // Note stride
297  typedef Vec<M,EStandard,1> TStandard;
298  typedef Row<M,EInvert,1> TInvert;
299  typedef Vec<M,ENormalize,1> TNormalize;
300 
301  typedef ESqHermT TSqHermT; // result of self dot product
302  typedef SymMat<M,ESqTHerm> TSqTHerm; // result of self outer product
303 
304  // These recurse right down to the underlying scalar type no matter how
305  // deep the elements are.
306  typedef EScalar Scalar;
307  typedef EULessScalar ULessScalar;
308  typedef ENumber Number;
309  typedef EStdNumber StdNumber;
310  typedef EPrecision Precision;
311  typedef EScalarNormSq ScalarNormSq;
316  static int size() { return M; }
318  static int nrow() { return M; }
320  static int ncol() { return 1; }
321 
322 
325  ScalarNormSq scalarNormSqr() const {
326  ScalarNormSq sum(0);
327  for(int i=0;i<M;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
328  return sum;
329  }
330 
335  TSqrt sqrt() const {
336  TSqrt vsqrt;
337  for(int i=0;i<M;++i) vsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
338  return vsqrt;
339  }
340 
345  TAbs abs() const {
346  TAbs vabs;
347  for(int i=0;i<M;++i) vabs[i] = CNT<E>::abs(d[i*STRIDE]);
348  return vabs;
349  }
350 
355  TStandard standardize() const {
356  TStandard vstd;
357  for(int i=0;i<M;++i) vstd[i] = CNT<E>::standardize(d[i*STRIDE]);
358  return vstd;
359  }
360 
364  EStandard sum() const {
365  E sum(0);
366  for (int i=0;i<M;++i) sum += d[i*STRIDE];
367  return CNT<E>::standardize(sum);
368  }
369 
370 
371  // This gives the resulting vector type when (v[i] op P) is applied to
372  // each element of v. It is a vector of length M, stride 1, and element
373  // types which are the regular composite result of E op P. Typically P is
374  // a scalar type but it doesn't have to be.
375  template <class P> struct EltResult {
376  typedef Vec<M, typename CNT<E>::template Result<P>::Mul, 1> Mul;
377  typedef Vec<M, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
378  typedef Vec<M, typename CNT<E>::template Result<P>::Add, 1> Add;
379  typedef Vec<M, typename CNT<E>::template Result<P>::Sub, 1> Sub;
380  };
381 
382  // This is the composite result for v op P where P is some kind of
383  // appropriately shaped non-scalar type.
384  template <class P> struct Result {
385  typedef MulCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
386  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
387  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
388  typedef typename MulOp::Type Mul;
389 
390  typedef MulCNTsNonConforming<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
391  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
392  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
393  typedef typename MulOpNonConforming::Type MulNon;
394 
395  typedef DvdCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
396  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
397  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
398  typedef typename DvdOp::Type Dvd;
399 
400  typedef AddCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
401  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
402  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
403  typedef typename AddOp::Type Add;
404 
405  typedef SubCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
406  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
407  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
408  typedef typename SubOp::Type Sub;
409  };
410 
416  template <class P> struct Substitute {
417  typedef Vec<M,P> Type;
418  };
419 
423  Vec(){
424  #ifndef NDEBUG
425  setToNaN();
426  #endif
427  }
428 
429  // It's important not to use the default copy constructor or copy
430  // assignment because the compiler doesn't understand that we may
431  // have noncontiguous storage and will try to copy the whole array.
432 
436  Vec(const Vec& src) {
437  Impl::copy(*this, src);
438  }
443  Vec& operator=(const Vec& src) {
444  Impl::copy(*this, src);
445  return *this;
446  }
447 
450  template <int SS> Vec(const Vec<M,E,SS>& src) {
451  Impl::copy(*this, src);
452  }
453 
456  template <int SS> Vec(const Vec<M,ENeg,SS>& src) {
457  Impl::copy(*this, src);
458  }
459 
462  template <class EE, int SS> explicit Vec(const Vec<M,EE,SS>& src) {
463  Impl::copy(*this, src);
464  }
465 
468  explicit Vec(const E& e) {for (int i=0;i<M;++i) d[i*STRIDE]=e;}
469 
474  explicit Vec(const ENeg& ne) {
475  const E e = ne; // requires floating point negation
476  for (int i=0;i<M;++i) d[i*STRIDE]=e;
477  }
478 
483  explicit Vec(int i) {new (this) Vec(E(Precision(i)));}
484 
485  // A bevy of constructors for Vecs up to length 9.
486 
488  Vec(const E& e0,const E& e1)
489  { assert(M==2);(*this)[0]=e0;(*this)[1]=e1; }
490  Vec(const E& e0,const E& e1,const E& e2)
491  { assert(M==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
492  Vec(const E& e0,const E& e1,const E& e2,const E& e3)
493  { assert(M==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
494  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
495  { assert(M==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
496  (*this)[3]=e3;(*this)[4]=e4; }
497  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
498  { assert(M==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
499  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
500  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6)
501  { assert(M==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
502  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
503  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)
504  { assert(M==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
505  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
506  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)
507  { assert(M==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
508  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
509 
514  template <class EE> explicit Vec(const EE* p)
515  { assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; }
516 
521  template <class EE> Vec& operator=(const EE* p)
522  { assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; return *this; }
523 
526  template <class EE, int SS> Vec& operator=(const Vec<M,EE,SS>& vv)
527  { Impl::copy(*this, vv); return *this; }
528 
531  template <class EE, int SS> Vec& operator+=(const Vec<M,EE,SS>& r)
532  { for(int i=0;i<M;++i) d[i*STRIDE] += r[i]; return *this; }
536  template <class EE, int SS> Vec& operator+=(const Vec<M,negator<EE>,SS>& r)
537  { for(int i=0;i<M;++i) d[i*STRIDE] -= -(r[i]); return *this; }
538 
541  template <class EE, int SS> Vec& operator-=(const Vec<M,EE,SS>& r)
542  { for(int i=0;i<M;++i) d[i*STRIDE] -= r[i]; return *this; }
546  template <class EE, int SS> Vec& operator-=(const Vec<M,negator<EE>,SS>& r)
547  { for(int i=0;i<M;++i) d[i*STRIDE] += -(r[i]); return *this; }
548 
549  // Conforming binary ops with 'this' on left, producing new packed result.
550  // Cases: v=v+v, v=v-v, m=v*r
551 
553  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Add>
554  conformingAdd(const Vec<M,EE,SS>& r) const {
555  Vec<M,typename CNT<E>::template Result<EE>::Add> result;
556  Impl::conformingAdd(*this, r, result);
557  return result;
558  }
560  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Sub>
561  conformingSubtract(const Vec<M,EE,SS>& r) const {
562  Vec<M,typename CNT<E>::template Result<EE>::Sub> result;
563  Impl::conformingSubtract(*this, r, result);
564  return result;
565  }
566 
569  template <class EE, int SS> Mat<M,M,typename CNT<E>::template Result<EE>::Mul>
570  conformingMultiply(const Row<M,EE,SS>& r) const {
571  Mat<M,M,typename CNT<E>::template Result<EE>::Mul> result;
572  for (int j=0;j<M;++j) result(j) = scalarMultiply(r(j));
573  return result;
574  }
575 
577  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Mul>
578  elementwiseMultiply(const Vec<M,EE,SS>& r) const {
579  Vec<M,typename CNT<E>::template Result<EE>::Mul> result;
580  Impl::elementwiseMultiply(*this, r, result);
581  return result;
582  }
584  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Dvd>
585  elementwiseDivide(const Vec<M,EE,SS>& r) const {
586  Vec<M,typename CNT<E>::template Result<EE>::Dvd> result;
587  Impl::elementwiseDivide(*this, r, result);
588  return result;
589  }
590 
594  const E& operator[](int i) const
595  { assert(0 <= i && i < M); return d[i*STRIDE]; }
597  const E& operator()(int i) const {return (*this)[i];}
598 
602  E& operator[](int i) {assert(0 <= i && i < M); return d[i*STRIDE];}
604  E& operator()(int i) {return (*this)[i];}
605 
606  ScalarNormSq normSqr() const { return scalarNormSqr(); }
607  typename CNT<ScalarNormSq>::TSqrt
608  norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
609 
621  TNormalize normalize() const {
622  if (CNT<E>::IsScalar) {
623  return castAwayNegatorIfAny() / (SignInterpretation*norm());
624  } else {
625  TNormalize elementwiseNormalized;
626  for (int i=0; i<M; ++i)
627  elementwiseNormalized[i] = CNT<E>::normalize((*this)[i]);
628  return elementwiseNormalized;
629  }
630  }
631 
633  TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
634 
636  const Vec& operator+() const { return *this; }
640  const TNeg& operator-() const { return negate(); }
643  TNeg& operator-() { return updNegate(); }
647  const THerm& operator~() const { return transpose(); }
651  THerm& operator~() { return updTranspose(); }
652 
654  const TNeg& negate() const { return *reinterpret_cast<const TNeg*>(this); }
657  TNeg& updNegate() { return *reinterpret_cast< TNeg*>(this); }
658 
660  const THerm& transpose() const { return *reinterpret_cast<const THerm*>(this); }
663  THerm& updTranspose() { return *reinterpret_cast< THerm*>(this); }
664 
669  const TPosTrans& positionalTranspose() const
670  { return *reinterpret_cast<const TPosTrans*>(this); }
672  TPosTrans& updPositionalTranspose()
673  { return *reinterpret_cast<TPosTrans*>(this); }
674 
679  const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
682  TReal& real() { return *reinterpret_cast< TReal*>(this); }
683 
684  // Had to contort these next two routines to get them through VC++ 7.net
685 
690  const TImag& imag() const {
691  const int offs = ImagOffset;
692  const EImag* p = reinterpret_cast<const EImag*>(this);
693  return *reinterpret_cast<const TImag*>(p+offs);
694  }
697  TImag& imag() {
698  const int offs = ImagOffset;
699  EImag* p = reinterpret_cast<EImag*>(this);
700  return *reinterpret_cast<TImag*>(p+offs);
701  }
702 
706  const TWithoutNegator& castAwayNegatorIfAny() const
707  { return *reinterpret_cast<const TWithoutNegator*>(this); }
710  TWithoutNegator& updCastAwayNegatorIfAny()
711  { return *reinterpret_cast<TWithoutNegator*>(this); }
712 
713  // These are elementwise binary operators, (this op ee) by default but
714  // (ee op this) if 'FromLeft' appears in the name. The result is a packed
715  // Vec<M> but the element type may change. These are mostly used to
716  // implement global operators. We call these "scalar" operators but
717  // actually the "scalar" can be a composite type.
718 
719  //TODO: consider converting 'e' to Standard Numbers as precalculation and
720  // changing return type appropriately.
721  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Mul>
722  scalarMultiply(const EE& e) const {
723  Vec<M, typename CNT<E>::template Result<EE>::Mul> result;
724  for (int i=0; i<M; ++i) result[i] = (*this)[i] * e;
725  return result;
726  }
727  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Mul>
728  scalarMultiplyFromLeft(const EE& e) const {
729  Vec<M, typename CNT<EE>::template Result<E>::Mul> result;
730  for (int i=0; i<M; ++i) result[i] = e * (*this)[i];
731  return result;
732  }
733 
734  // TODO: should precalculate and store 1/e, while converting to Standard
735  // Numbers. Note that return type should change appropriately.
736  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Dvd>
737  scalarDivide(const EE& e) const {
738  Vec<M, typename CNT<E>::template Result<EE>::Dvd> result;
739  for (int i=0; i<M; ++i) result[i] = (*this)[i] / e;
740  return result;
741  }
742  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Dvd>
743  scalarDivideFromLeft(const EE& e) const {
744  Vec<M, typename CNT<EE>::template Result<E>::Dvd> result;
745  for (int i=0; i<M; ++i) result[i] = e / (*this)[i];
746  return result;
747  }
748 
749  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Add>
750  scalarAdd(const EE& e) const {
751  Vec<M, typename CNT<E>::template Result<EE>::Add> result;
752  for (int i=0; i<M; ++i) result[i] = (*this)[i] + e;
753  return result;
754  }
755  // Add is commutative, so no 'FromLeft'.
756 
757  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Sub>
758  scalarSubtract(const EE& e) const {
759  Vec<M, typename CNT<E>::template Result<EE>::Sub> result;
760  for (int i=0; i<M; ++i) result[i] = (*this)[i] - e;
761  return result;
762  }
763  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Sub>
764  scalarSubtractFromLeft(const EE& e) const {
765  Vec<M, typename CNT<EE>::template Result<E>::Sub> result;
766  for (int i=0; i<M; ++i) result[i] = e - (*this)[i];
767  return result;
768  }
769 
770  // Generic assignments for any element type not listed explicitly, including scalars.
771  // These are done repeatedly for each element and only work if the operation can
772  // be performed leaving the original element type.
773  template <class EE> Vec& operator =(const EE& e) {return scalarEq(e);}
774  template <class EE> Vec& operator+=(const EE& e) {return scalarPlusEq(e);}
775  template <class EE> Vec& operator-=(const EE& e) {return scalarMinusEq(e);}
776  template <class EE> Vec& operator*=(const EE& e) {return scalarTimesEq(e);}
777  template <class EE> Vec& operator/=(const EE& e) {return scalarDivideEq(e);}
778 
779  // Generalized element assignment & computed assignment methods. These will work
780  // for any assignment-compatible element, not just scalars.
781  template <class EE> Vec& scalarEq(const EE& ee)
782  { for(int i=0;i<M;++i) d[i*STRIDE] = ee; return *this; }
783  template <class EE> Vec& scalarPlusEq(const EE& ee)
784  { for(int i=0;i<M;++i) d[i*STRIDE] += ee; return *this; }
785  template <class EE> Vec& scalarMinusEq(const EE& ee)
786  { for(int i=0;i<M;++i) d[i*STRIDE] -= ee; return *this; }
787  template <class EE> Vec& scalarMinusEqFromLeft(const EE& ee)
788  { for(int i=0;i<M;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
789  template <class EE> Vec& scalarTimesEq(const EE& ee)
790  { for(int i=0;i<M;++i) d[i*STRIDE] *= ee; return *this; }
791  template <class EE> Vec& scalarTimesEqFromLeft(const EE& ee)
792  { for(int i=0;i<M;++i) d[i*STRIDE] = ee * d[i*STRIDE]; return *this; }
793  template <class EE> Vec& scalarDivideEq(const EE& ee)
794  { for(int i=0;i<M;++i) d[i*STRIDE] /= ee; return *this; }
795  template <class EE> Vec& scalarDivideEqFromLeft(const EE& ee)
796  { for(int i=0;i<M;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
797 
798  // Specialize for int to avoid warnings and ambiguities.
799  Vec& scalarEq(int ee) {return scalarEq(Precision(ee));}
800  Vec& scalarPlusEq(int ee) {return scalarPlusEq(Precision(ee));}
801  Vec& scalarMinusEq(int ee) {return scalarMinusEq(Precision(ee));}
802  Vec& scalarTimesEq(int ee) {return scalarTimesEq(Precision(ee));}
803  Vec& scalarDivideEq(int ee) {return scalarDivideEq(Precision(ee));}
804  Vec& scalarMinusEqFromLeft(int ee) {return scalarMinusEqFromLeft(Precision(ee));}
805  Vec& scalarTimesEqFromLeft(int ee) {return scalarTimesEqFromLeft(Precision(ee));}
806  Vec& scalarDivideEqFromLeft(int ee) {return scalarDivideEqFromLeft(Precision(ee));}
807 
810  void setToNaN() {
811  (*this) = CNT<ELT>::getNaN();
812  }
813 
815  void setToZero() {
816  (*this) = ELT(0);
817  }
818 
824  template <int MM>
825  const Vec<MM,ELT,STRIDE>& getSubVec(int i) const {
826  assert(0 <= i && i + MM <= M);
827  return Vec<MM,ELT,STRIDE>::getAs(&(*this)[i]);
828  }
834  template <int MM>
835  Vec<MM,ELT,STRIDE>& updSubVec(int i) {
836  assert(0 <= i && i + MM <= M);
837  return Vec<MM,ELT,STRIDE>::updAs(&(*this)[i]);
838  }
839 
840 
844  template <int MM>
845  static const Vec& getSubVec(const Vec<MM,ELT,STRIDE>& v, int i) {
846  assert(0 <= i && i + M <= MM);
847  return getAs(&v[i]);
848  }
852  template <int MM>
853  static Vec& updSubVec(Vec<MM,ELT,STRIDE>& v, int i) {
854  assert(0 <= i && i + M <= MM);
855  return updAs(&v[i]);
856  }
857 
861  Vec<M-1,ELT,1> drop1(int p) const {
862  assert(0 <= p && p < M);
863  Vec<M-1,ELT,1> out;
864  int nxt=0;
865  for (int i=0; i<M-1; ++i, ++nxt) {
866  if (nxt==p) ++nxt; // skip the loser
867  out[i] = (*this)[nxt];
868  }
869  return out;
870  }
871 
875  template <class EE> Vec<M+1,ELT,1> append1(const EE& v) const {
876  Vec<M+1,ELT,1> out;
877  Vec<M,ELT,1>::updAs(&out[0]) = (*this);
878  out[M] = v;
879  return out;
880  }
881 
882 
888  template <class EE> Vec<M+1,ELT,1> insert1(int p, const EE& v) const {
889  assert(0 <= p && p <= M);
890  if (p==M) return append1(v);
891  Vec<M+1,ELT,1> out;
892  int nxt=0;
893  for (int i=0; i<M; ++i, ++nxt) {
894  if (i==p) out[nxt++] = v;
895  out[nxt] = (*this)[i];
896  }
897  return out;
898  }
899 
902  static const Vec& getAs(const ELT* p)
903  { return *reinterpret_cast<const Vec*>(p); }
906  static Vec& updAs(ELT* p)
907  { return *reinterpret_cast<Vec*>(p); }
908 
909 
913  static Vec<M,ELT,1> getNaN() { return Vec<M,ELT,1>(CNT<ELT>::getNaN()); }
914 
916  bool isNaN() const {
917  for (int i=0; i<M; ++i)
918  if (CNT<ELT>::isNaN((*this)[i]))
919  return true;
920  return false;
921  }
922 
925  bool isInf() const {
926  bool seenInf = false;
927  for (int i=0; i<M; ++i) {
928  const ELT& e = (*this)[i];
929  if (!CNT<ELT>::isFinite(e)) {
930  if (!CNT<ELT>::isInf(e))
931  return false; // something bad was found
932  seenInf = true;
933  }
934  }
935  return seenInf;
936  }
937 
940  bool isFinite() const {
941  for (int i=0; i<M; ++i)
942  if (!CNT<ELT>::isFinite((*this)[i]))
943  return false;
944  return true;
945  }
946 
949  static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
950 
953  template <class E2, int RS2>
954  bool isNumericallyEqual(const Vec<M,E2,RS2>& v, double tol) const {
955  for (int i=0; i<M; ++i)
956  if (!CNT<ELT>::isNumericallyEqual((*this)[i], v[i], tol))
957  return false;
958  return true;
959  }
960 
964  template <class E2, int RS2>
965  bool isNumericallyEqual(const Vec<M,E2,RS2>& v) const {
966  const double tol = std::max(getDefaultTolerance(),v.getDefaultTolerance());
967  return isNumericallyEqual(v, tol);
968  }
969 
974  bool isNumericallyEqual
975  (const ELT& e,
976  double tol = getDefaultTolerance()) const
977  {
978  for (int i=0; i<M; ++i)
979  if (!CNT<ELT>::isNumericallyEqual((*this)[i], e, tol))
980  return false;
981  return true;
982  }
983 
984  // Functions to be used for Scripting in MATLAB and languages that do not support operator overloading
986  std::string toString() const {
987  std::stringstream stream;
988  stream << (*this);
989  return stream.str();
990  }
991 
993  void set(int i, const E& value)
994  { (*this)[i] = value; }
995 
997  const E& get(int i) const
998  { return operator[](i); }
999 
1000 private:
1001  // TODO: should be an array of scalars rather than elements to control
1002  // packing more carefully.
1003  ELT d[NActualElements]; // data
1004 };
1005 
1007 // Global operators involving two vectors. //
1008 // v+v, v-v, v==v, v!=v //
1010 
1011 // v3 = v1 + v2 where all v's have the same length M.
1012 template <int M, class E1, int S1, class E2, int S2> inline
1013 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Add
1014 operator+(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {
1015  return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
1016  ::AddOp::perform(l,r);
1017 }
1018 
1019 // v3 = v1 - v2, similar to +
1020 template <int M, class E1, int S1, class E2, int S2> inline
1021 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Sub
1022 operator-(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {
1023  return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
1024  ::SubOp::perform(l,r);
1025 }
1026 
1028 template <int M, class E1, int S1, class E2, int S2> inline bool
1029 operator==(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1030 { for (int i=0; i < M; ++i) if (l[i] != r[i]) return false;
1031  return true; }
1033 template <int M, class E1, int S1, class E2, int S2> inline bool
1034 operator!=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {return !(l==r);}
1035 
1037 template <int M, class E1, int S1, class E2> inline bool
1038 operator==(const Vec<M,E1,S1>& v, const E2& e)
1039 { for (int i=0; i < M; ++i) if (v[i] != e) return false;
1040  return true; }
1042 template <int M, class E1, int S1, class E2> inline bool
1043 operator!=(const Vec<M,E1,S1>& v, const E2& e) {return !(v==e);}
1044 
1046 template <int M, class E1, int S1, class E2, int S2> inline bool
1047 operator<(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1048 { for (int i=0; i < M; ++i) if (l[i] >= r[i]) return false;
1049  return true; }
1051 template <int M, class E1, int S1, class E2> inline bool
1052 operator<(const Vec<M,E1,S1>& v, const E2& e)
1053 { for (int i=0; i < M; ++i) if (v[i] >= e) return false;
1054  return true; }
1055 
1057 template <int M, class E1, int S1, class E2, int S2> inline bool
1058 operator>(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1059 { for (int i=0; i < M; ++i) if (l[i] <= r[i]) return false;
1060  return true; }
1062 template <int M, class E1, int S1, class E2> inline bool
1063 operator>(const Vec<M,E1,S1>& v, const E2& e)
1064 { for (int i=0; i < M; ++i) if (v[i] <= e) return false;
1065  return true; }
1066 
1069 template <int M, class E1, int S1, class E2, int S2> inline bool
1070 operator<=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1071 { for (int i=0; i < M; ++i) if (l[i] > r[i]) return false;
1072  return true; }
1075 template <int M, class E1, int S1, class E2> inline bool
1076 operator<=(const Vec<M,E1,S1>& v, const E2& e)
1077 { for (int i=0; i < M; ++i) if (v[i] > e) return false;
1078  return true; }
1079 
1082 template <int M, class E1, int S1, class E2, int S2> inline bool
1083 operator>=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1084 { for (int i=0; i < M; ++i) if (l[i] < r[i]) return false;
1085  return true; }
1088 template <int M, class E1, int S1, class E2> inline bool
1089 operator>=(const Vec<M,E1,S1>& v, const E2& e)
1090 { for (int i=0; i < M; ++i) if (v[i] < e) return false;
1091  return true; }
1092 
1094 // Global operators involving a vector and a scalar. //
1096 
1097 // I haven't been able to figure out a nice way to templatize for the
1098 // built-in reals without introducing a lot of unwanted type matches
1099 // as well. So we'll just grind them out explicitly here.
1100 
1101 // SCALAR MULTIPLY
1102 
1103 // v = v*real, real*v
1104 template <int M, class E, int S> inline
1105 typename Vec<M,E,S>::template Result<float>::Mul
1106 operator*(const Vec<M,E,S>& l, const float& r)
1107  { return Vec<M,E,S>::template Result<float>::MulOp::perform(l,r); }
1108 template <int M, class E, int S> inline
1109 typename Vec<M,E,S>::template Result<float>::Mul
1110 operator*(const float& l, const Vec<M,E,S>& r) {return r*l;}
1111 
1112 template <int M, class E, int S> inline
1113 typename Vec<M,E,S>::template Result<double>::Mul
1114 operator*(const Vec<M,E,S>& l, const double& r)
1115  { return Vec<M,E,S>::template Result<double>::MulOp::perform(l,r); }
1116 template <int M, class E, int S> inline
1117 typename Vec<M,E,S>::template Result<double>::Mul
1118 operator*(const double& l, const Vec<M,E,S>& r) {return r*l;}
1119 
1120 template <int M, class E, int S> inline
1121 typename Vec<M,E,S>::template Result<long double>::Mul
1122 operator*(const Vec<M,E,S>& l, const long double& r)
1123  { return Vec<M,E,S>::template Result<long double>::MulOp::perform(l,r); }
1124 template <int M, class E, int S> inline
1125 typename Vec<M,E,S>::template Result<long double>::Mul
1126 operator*(const long double& l, const Vec<M,E,S>& r) {return r*l;}
1127 
1128 // v = v*int, int*v -- just convert int to v's precision float
1129 template <int M, class E, int S> inline
1130 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
1131 operator*(const Vec<M,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
1132 template <int M, class E, int S> inline
1133 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
1134 operator*(int l, const Vec<M,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
1135 
1136 // Complex, conjugate, and negator are all easy to templatize.
1137 
1138 // v = v*complex, complex*v
1139 template <int M, class E, int S, class R> inline
1140 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1141 operator*(const Vec<M,E,S>& l, const std::complex<R>& r)
1142  { return Vec<M,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
1143 template <int M, class E, int S, class R> inline
1144 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1145 operator*(const std::complex<R>& l, const Vec<M,E,S>& r) {return r*l;}
1146 
1147 // v = v*conjugate, conjugate*v (convert conjugate->complex)
1148 template <int M, class E, int S, class R> inline
1149 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1150 operator*(const Vec<M,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
1151 template <int M, class E, int S, class R> inline
1152 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1153 operator*(const conjugate<R>& l, const Vec<M,E,S>& r) {return r*(std::complex<R>)l;}
1154 
1155 // v = v*negator, negator*v: convert negator to standard number
1156 template <int M, class E, int S, class R> inline
1157 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
1158 operator*(const Vec<M,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
1159 template <int M, class E, int S, class R> inline
1160 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
1161 operator*(const negator<R>& l, const Vec<M,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
1162 
1163 
1164 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
1165 // but when it is on the left it means scalar * pseudoInverse(vec), which is
1166 // a row.
1167 
1168 // v = v/real, real/v
1169 template <int M, class E, int S> inline
1170 typename Vec<M,E,S>::template Result<float>::Dvd
1171 operator/(const Vec<M,E,S>& l, const float& r)
1172  { return Vec<M,E,S>::template Result<float>::DvdOp::perform(l,r); }
1173 template <int M, class E, int S> inline
1174 typename CNT<float>::template Result<Vec<M,E,S> >::Dvd
1175 operator/(const float& l, const Vec<M,E,S>& r)
1176  { return CNT<float>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1177 
1178 template <int M, class E, int S> inline
1179 typename Vec<M,E,S>::template Result<double>::Dvd
1180 operator/(const Vec<M,E,S>& l, const double& r)
1181  { return Vec<M,E,S>::template Result<double>::DvdOp::perform(l,r); }
1182 template <int M, class E, int S> inline
1183 typename CNT<double>::template Result<Vec<M,E,S> >::Dvd
1184 operator/(const double& l, const Vec<M,E,S>& r)
1185  { return CNT<double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1186 
1187 template <int M, class E, int S> inline
1188 typename Vec<M,E,S>::template Result<long double>::Dvd
1189 operator/(const Vec<M,E,S>& l, const long double& r)
1190  { return Vec<M,E,S>::template Result<long double>::DvdOp::perform(l,r); }
1191 template <int M, class E, int S> inline
1192 typename CNT<long double>::template Result<Vec<M,E,S> >::Dvd
1193 operator/(const long double& l, const Vec<M,E,S>& r)
1194  { return CNT<long double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1195 
1196 // v = v/int, int/v -- just convert int to v's precision float
1197 template <int M, class E, int S> inline
1198 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Dvd
1199 operator/(const Vec<M,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
1200 template <int M, class E, int S> inline
1201 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Dvd
1202 operator/(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
1203 
1204 
1205 // Complex, conjugate, and negator are all easy to templatize.
1206 
1207 // v = v/complex, complex/v
1208 template <int M, class E, int S, class R> inline
1209 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
1210 operator/(const Vec<M,E,S>& l, const std::complex<R>& r)
1211  { return Vec<M,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
1212 template <int M, class E, int S, class R> inline
1213 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
1214 operator/(const std::complex<R>& l, const Vec<M,E,S>& r)
1215  { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1216 
1217 // v = v/conjugate, conjugate/v (convert conjugate->complex)
1218 template <int M, class E, int S, class R> inline
1219 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
1220 operator/(const Vec<M,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
1221 template <int M, class E, int S, class R> inline
1222 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
1223 operator/(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l/r;}
1224 
1225 // v = v/negator, negator/v: convert negator to number
1226 template <int M, class E, int S, class R> inline
1227 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
1228 operator/(const Vec<M,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
1229 template <int M, class E, int S, class R> inline
1230 typename CNT<R>::template Result<Vec<M,E,S> >::Dvd
1231 operator/(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
1232 
1233 
1234 // Add and subtract are odd as scalar ops. They behave as though the
1235 // scalar stands for a vector each of whose elements is that scalar,
1236 // and then a normal vector add or subtract is done.
1237 
1238 // SCALAR ADD
1239 
1240 // v = v+real, real+v
1241 template <int M, class E, int S> inline
1242 typename Vec<M,E,S>::template Result<float>::Add
1243 operator+(const Vec<M,E,S>& l, const float& r)
1244  { return Vec<M,E,S>::template Result<float>::AddOp::perform(l,r); }
1245 template <int M, class E, int S> inline
1246 typename Vec<M,E,S>::template Result<float>::Add
1247 operator+(const float& l, const Vec<M,E,S>& r) {return r+l;}
1248 
1249 template <int M, class E, int S> inline
1250 typename Vec<M,E,S>::template Result<double>::Add
1251 operator+(const Vec<M,E,S>& l, const double& r)
1252  { return Vec<M,E,S>::template Result<double>::AddOp::perform(l,r); }
1253 template <int M, class E, int S> inline
1254 typename Vec<M,E,S>::template Result<double>::Add
1255 operator+(const double& l, const Vec<M,E,S>& r) {return r+l;}
1256 
1257 template <int M, class E, int S> inline
1258 typename Vec<M,E,S>::template Result<long double>::Add
1259 operator+(const Vec<M,E,S>& l, const long double& r)
1260  { return Vec<M,E,S>::template Result<long double>::AddOp::perform(l,r); }
1261 template <int M, class E, int S> inline
1262 typename Vec<M,E,S>::template Result<long double>::Add
1263 operator+(const long double& l, const Vec<M,E,S>& r) {return r+l;}
1264 
1265 // v = v+int, int+v -- just convert int to v's precision float
1266 template <int M, class E, int S> inline
1267 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
1268 operator+(const Vec<M,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
1269 template <int M, class E, int S> inline
1270 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
1271 operator+(int l, const Vec<M,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
1272 
1273 // Complex, conjugate, and negator are all easy to templatize.
1274 
1275 // v = v+complex, complex+v
1276 template <int M, class E, int S, class R> inline
1277 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1278 operator+(const Vec<M,E,S>& l, const std::complex<R>& r)
1279  { return Vec<M,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
1280 template <int M, class E, int S, class R> inline
1281 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1282 operator+(const std::complex<R>& l, const Vec<M,E,S>& r) {return r+l;}
1283 
1284 // v = v+conjugate, conjugate+v (convert conjugate->complex)
1285 template <int M, class E, int S, class R> inline
1286 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1287 operator+(const Vec<M,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
1288 template <int M, class E, int S, class R> inline
1289 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1290 operator+(const conjugate<R>& l, const Vec<M,E,S>& r) {return r+(std::complex<R>)l;}
1291 
1292 // v = v+negator, negator+v: convert negator to standard number
1293 template <int M, class E, int S, class R> inline
1294 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
1295 operator+(const Vec<M,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
1296 template <int M, class E, int S, class R> inline
1297 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
1298 operator+(const negator<R>& l, const Vec<M,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
1299 
1300 // SCALAR SUBTRACT -- careful, not commutative.
1301 
1302 // v = v-real, real-v
1303 template <int M, class E, int S> inline
1304 typename Vec<M,E,S>::template Result<float>::Sub
1305 operator-(const Vec<M,E,S>& l, const float& r)
1306  { return Vec<M,E,S>::template Result<float>::SubOp::perform(l,r); }
1307 template <int M, class E, int S> inline
1308 typename CNT<float>::template Result<Vec<M,E,S> >::Sub
1309 operator-(const float& l, const Vec<M,E,S>& r)
1310  { return CNT<float>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1311 
1312 template <int M, class E, int S> inline
1313 typename Vec<M,E,S>::template Result<double>::Sub
1314 operator-(const Vec<M,E,S>& l, const double& r)
1315  { return Vec<M,E,S>::template Result<double>::SubOp::perform(l,r); }
1316 template <int M, class E, int S> inline
1317 typename CNT<double>::template Result<Vec<M,E,S> >::Sub
1318 operator-(const double& l, const Vec<M,E,S>& r)
1319  { return CNT<double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1320 
1321 template <int M, class E, int S> inline
1322 typename Vec<M,E,S>::template Result<long double>::Sub
1323 operator-(const Vec<M,E,S>& l, const long double& r)
1324  { return Vec<M,E,S>::template Result<long double>::SubOp::perform(l,r); }
1325 template <int M, class E, int S> inline
1326 typename CNT<long double>::template Result<Vec<M,E,S> >::Sub
1327 operator-(const long double& l, const Vec<M,E,S>& r)
1328  { return CNT<long double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1329 
1330 // v = v-int, int-v // just convert int to v's precision float
1331 template <int M, class E, int S> inline
1332 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Sub
1333 operator-(const Vec<M,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
1334 template <int M, class E, int S> inline
1335 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Sub
1336 operator-(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
1337 
1338 
1339 // Complex, conjugate, and negator are all easy to templatize.
1340 
1341 // v = v-complex, complex-v
1342 template <int M, class E, int S, class R> inline
1343 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
1344 operator-(const Vec<M,E,S>& l, const std::complex<R>& r)
1345  { return Vec<M,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
1346 template <int M, class E, int S, class R> inline
1347 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
1348 operator-(const std::complex<R>& l, const Vec<M,E,S>& r)
1349  { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1350 
1351 // v = v-conjugate, conjugate-v (convert conjugate->complex)
1352 template <int M, class E, int S, class R> inline
1353 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
1354 operator-(const Vec<M,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
1355 template <int M, class E, int S, class R> inline
1356 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
1357 operator-(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l-r;}
1358 
1359 // v = v-negator, negator-v: convert negator to standard number
1360 template <int M, class E, int S, class R> inline
1361 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Sub
1362 operator-(const Vec<M,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
1363 template <int M, class E, int S, class R> inline
1364 typename CNT<R>::template Result<Vec<M,E,S> >::Sub
1365 operator-(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
1366 
1367 // Vec I/O
1368 template <int M, class E, int S, class CHAR, class TRAITS> inline
1369 std::basic_ostream<CHAR,TRAITS>&
1370 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Vec<M,E,S>& v) {
1371  o << "~[" << v[0]; for(int i=1;i<M;++i) o<<','<<v[i]; o<<']'; return o;
1372 }
1373 
1376 template <int M, class E, int S, class CHAR, class TRAITS> inline
1377 std::basic_istream<CHAR,TRAITS>&
1378 operator>>(std::basic_istream<CHAR,TRAITS>& is, Vec<M,E,S>& v) {
1379  CHAR tilde;
1380  is >> tilde; if (is.fail()) return is;
1381  if (tilde != CHAR('~')) {
1382  tilde = CHAR(0);
1383  is.unget(); if (is.fail()) return is;
1384  }
1385 
1386  CHAR openBracket, closeBracket;
1387  is >> openBracket; if (is.fail()) return is;
1388  if (openBracket==CHAR('('))
1389  closeBracket = CHAR(')');
1390  else if (openBracket==CHAR('['))
1391  closeBracket = CHAR(']');
1392  else {
1393  closeBracket = CHAR(0);
1394  is.unget(); if (is.fail()) return is;
1395  }
1396 
1397  // If we saw a "~" but then we didn't see any brackets, that's an
1398  // error. Set the fail bit and return.
1399  if (tilde != CHAR(0) && closeBracket == CHAR(0)) {
1400  is.setstate( std::ios::failbit );
1401  return is;
1402  }
1403 
1404  for (int i=0; i < M; ++i) {
1405  is >> v[i];
1406  if (is.fail()) return is;
1407  if (i != M-1) {
1408  CHAR c; is >> c; if (is.fail()) return is;
1409  if (c != ',') is.unget();
1410  if (is.fail()) return is;
1411  }
1412  }
1413 
1414  // Get the closing bracket if there was an opening one. If we don't
1415  // see the expected character we'll set the fail bit in the istream.
1416  if (closeBracket != CHAR(0)) {
1417  CHAR closer; is >> closer; if (is.fail()) return is;
1418  if (closer != closeBracket) {
1419  is.unget(); if (is.fail()) return is;
1420  is.setstate( std::ios::failbit );
1421  }
1422  }
1423 
1424  return is;
1425 }
1426 
1427 } //namespace SimTK
1428 
1429 
1430 #endif //SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_