CarpetLib: Add placeholders for "const" and "pure" attributes to vect<T>

Add __attribute__((const)) and __attribute__((pure)) to vect<T>
functions where appropriate.  These attributes are not active, because
it is not clear that this really improves performance, and errors
could lead to wrong code.

darcs-hash:20060731152218-dae7b-20c523f60890b2ad611157764e0b993bfb185c8d.gz

1 files changed, 233 insertions, 59 deletions

diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index bc9f88cf6..f58f1e07b 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -12,6 +12,37 @@ using namespace std;
 
 
 
+#if 0
+
+// A pure function returns a value that depends only on the function
+// arguments and on global variables, and the function has no side
+// effects.
+#ifdef HAVE_CCTK_CXX_ATTRIBUTE_PURE
+#  define PURE __attribute__((pure))
+#else
+#  define PURE
+#endif
+
+// A const function returns a value that depends only on the function
+// arguments, and the function has no side effects.  This is even more
+// strict than pure functions.  Const functions cannot dereference
+// pointers or references (or this).
+#ifdef HAVE_CCTK_CXX_ATTRIBUTE_CONST
+#  define CONST __attribute__((const))
+#else
+#  define CONST
+#endif
+
+#else
+
+// Don't take any risks
+#  define PURE
+#  define CONST
+
+#endif
+
+
+
 // Forward definition
 template<typename T, int D> class vect;
 
@@ -39,21 +70,21 @@ public:
   // Constructors
   
   /** Explicit empty constructor.  */
-  explicit vect () { }
+  explicit vect () CONST { }
   
   /** Copy constructor.  */
-  vect (const vect& a) {
+  vect (const vect& a) PURE {
     for (int d=0; d<D; ++d) elt[d]=a.elt[d];
   }
   
   /** Constructor from a single element.  This constructor might be
       confusing, but it is very convenient.  */
-  vect (const T& x) {
+  vect (const T& x) PURE {
     for (int d=0; d<D; ++d) elt[d]=x;
   }
   
   /** Constructor for 2-element vectors from 2 elements.  */
-  vect (const T& x, const T& y) {
+  vect (const T& x, const T& y) PURE {
     assert (D==2);
     // Note: this statement may give "index out of range" warnings.
     // You can safely ignore these.
@@ -61,7 +92,7 @@ public:
   }
   
   /** Constructor for 3-element vectors from 3 elements.  */
-  vect (const T& x, const T& y, const T& z) {
+  vect (const T& x, const T& y, const T& z) PURE {
     assert (D==3);
     // Note: this statement may give "index out of range" warnings.
     // You can safely ignore these.
@@ -69,7 +100,7 @@ public:
   }
   
   /** Constructor for 4-element vectors from 4 elements.  */
-  vect (const T& x, const T& y, const T& z, const T& t) {
+  vect (const T& x, const T& y, const T& z, const T& t) PURE {
     assert (D==4);
     // Note: this statement may give "index out of range" warnings.
     // You can safely ignore these.
@@ -79,81 +110,81 @@ public:
 #if 0
   // This creates confusion
   /** Constructor from a pointer, i.e.\ a C array.  */
-  explicit vect (const T* const x) {
+  explicit vect (const T* const x) PURE {
     for (int d=0; d<D; ++d) elt[d]=x[d];
   }
 #endif
   
   /** Constructor from a vector with a different type.  */
   template<typename S>
-  /*explicit*/ vect (const vect<S,D>& a) {
+  /*explicit*/ vect (const vect<S,D>& a) /*PURE*/ {
     for (int d=0; d<D; ++d) elt[d]=(T)a[d];
   }
   
   /** Create a new 0-element vector with a specific type.  */
-  static vect make () {
+  static vect make () CONST {
     assert (D==0);
     return vect();
   }
   
   /** Create a new 1-element vector with a specific type.  */
-  static vect make (const T& x) {
+  static vect make (const T& x) PURE {
     assert (D==1);
     return vect(x);
   }
   
   /** Create a new 2-element vector with a specific type.  */
-  static vect make (const T& x, const T& y) {
+  static vect make (const T& x, const T& y) PURE {
     assert (D==2);
     return vect(x, y);
   }
   
   /** Create a new 3-element vector with a specific type.  */
-  static vect make (const T& x, const T& y, const T& z) {
+  static vect make (const T& x, const T& y, const T& z) PURE {
     assert (D==3);
     return vect(x, y, z);
   }
   
   /** Create a new 4-element vector with a specific type.  */
-  static vect make (const T& x, const T& y, const T& z, const T& t) {
+  static vect make (const T& x, const T& y, const T& z, const T& t) PURE {
     assert (D==4);
     return vect(x, y, z, t);
   }
   
   /** Treat a constant pointer as a reference to a constant vector.  */
-  static const vect& ref (const T* const x) {
+  static const vect& ref (const T* const x) PURE {
     return *(const vect*)x;
   }
   
   /** Treat a pointer as a reference to a vector.  */
-  static vect& ref (T* const x) {
+  static vect& ref (T* const x) PURE {
     return *(vect*)x;
   }
   
   /** Create a vector with one element set to 1 and all other elements
       set to zero.  */
-  static vect dir (const int d) {
+  static vect dir (const int d) CONST {
     vect r=(T)0;
     r[d]=1;
     return r;
   }
   
   /** Create a vector with e[i] = i.  */
-  static vect seq () {
+  static vect seq () CONST {
     vect r;
     for (int d=0; d<D; ++d) r[d]=d;
     return r;
   }
   
   /** Create a vector with e[i] = n + i.  */
-  static vect seq (const int n) {
+  static vect seq (const int n) CONST {
     vect r;
     for (int d=0; d<D; ++d) r[d]=n+d;
     return r;
   }
   
   /** Create a vector with e[i] = n + s * i.  */
-  static vect seq (const int n, const int s) {
+  static vect seq (const int n, const int s) CONST {
     vect r;
     for (int d=0; d<D; ++d) r[d]=n+s*d;
     return r;
@@ -164,13 +195,13 @@ public:
   /** Return a non-writable element of a vector.  */
   // (Don't return a reference; *this might be a temporary)
   // Do return a reference, so that a vector can be accessed as array
-  const T& operator[] (const int d) const {
+  const T& operator[] (const int d) const PURE {
     assert(d>=0 && d<D);
     return elt[d];
   }
   
   /** Return a writable element of a vector as reference.  */
-  T& operator[] (const int d) {
+  T& operator[] (const int d) PURE {
     assert(d>=0 && d<D);
     return elt[d];
   }
@@ -178,7 +209,7 @@ public:
 #if 0
   // This creates confusion
   /** Return a pointer to a vector.  */
-  operator const T* () const {
+  operator const T* () const PURE {
     return this;
   }
 #endif
@@ -186,7 +217,7 @@ public:
   /** Return a combination of the vector elements e[a[i]].  The
       element combination is selected by another vector.  */
   template<typename TT, int DD>
-  vect<T,DD> operator[] (const vect<TT,DD>& a) const {
+  vect<T,DD> operator[] (const vect<TT,DD>& a) const /*PURE*/ {
     vect<T,DD> r;
     // (*this)[] performs index checking
     for (int d=0; d<DD; ++d) r[d] = (*this)[a[d]];
@@ -277,38 +308,38 @@ public:
   // Non-modifying operators
   
   /** Return a new vector where one element has been replaced.  */
-  vect replace (const int d, const T& x) const {
+  vect replace (const int d, const T& x) const PURE {
     assert (d>=0 && d<D);
     vect r;
     for (int dd=0; dd<D; ++dd) r[dd]=dd==d?x:elt[dd];
     return r;
   }
   
-  vect reverse () const {
+  vect reverse () const PURE {
     vect r;
     for (int d=0; d<D; ++d) r[d]=elt[D-1-d];
     return r;
   }
   
-  vect operator+ () const {
+  vect operator+ () const PURE {
     vect r;
     for (int d=0; d<D; ++d) r[d]=+elt[d];
     return r;
   }
   
-  vect operator- () const {
+  vect operator- () const PURE {
     vect r;
     for (int d=0; d<D; ++d) r[d]=-elt[d];
     return r;
   }
   
-  vect<bool,D> operator! () const {
+  vect<bool,D> operator! () const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=!elt[d];
     return r;
   }
   
-  vect operator~ () const {
+  vect operator~ () const PURE {
     vect r;
     for (int d=0; d<D; ++d) r[d]=~elt[d];
     return r;
@@ -316,157 +347,157 @@ public:
   
 #if 0
   // These are now defined as template functions
-  vect operator+ (const T& x) const {
+  vect operator+ (const T& x) const PURE {
     vect r(*this);
     r+=x;
     return r;
   }
   
-  vect operator- (const T& x) const {
+  vect operator- (const T& x) const PURE {
     vect r(*this);
     r-=x;
     return r;
   }
   
-  vect operator* (const T& x) const {
+  vect operator* (const T& x) const PURE {
     vect r(*this);
     r*=x;
     return r;
   }
   
-  vect operator/ (const T& x) const {
+  vect operator/ (const T& x) const PURE {
     vect r(*this);
     r/=x;
     return r;
   }
   
-  vect operator% (const T& x) const {
+  vect operator% (const T& x) const PURE {
     vect r(*this);
     r%=x;
     return r;
   }
   
-  vect operator& (const T& x) const {
+  vect operator& (const T& x) const PURE {
     vect r(*this);
     r&=x;
     return r;
   }
   
-  vect operator| (const T& x) const {
+  vect operator| (const T& x) const PURE {
     vect r(*this);
     r|=x;
     return r;
   }
   
-  vect operator^ (const T& x) const {
+  vect operator^ (const T& x) const PURE {
     vect r(*this);
     r^=x;
     return r;
   }
   
-  vect<bool,D> operator&& (const T& x) const {
+  vect<bool,D> operator&& (const T& x) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]&&x;
     return r;
   }
   
-  vect<bool,D> operator|| (const T& x) const {
+  vect<bool,D> operator|| (const T& x) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]||x;
     return r;
   }
   
-  vect operator+ (const vect& a) const {
+  vect operator+ (const vect& a) const PURE {
     vect r(*this);
     r+=a;
     return r;
   }
   
-  vect operator- (const vect& a) const {
+  vect operator- (const vect& a) const PURE {
     vect r(*this);
     r-=a;
     return r;
   }
   
-  vect operator* (const vect& a) const {
+  vect operator* (const vect& a) const PURE {
     vect r(*this);
     r*=a;
     return r;
   }
   
-  vect operator/ (const vect& a) const {
+  vect operator/ (const vect& a) const PURE {
     vect r(*this);
     r/=a;
     return r;
   }
   
-  vect operator% (const vect& a) const {
+  vect operator% (const vect& a) const PURE {
     vect r(*this);
     r%=a;
     return r;
   }
   
-  vect operator& (const vect& a) const {
+  vect operator& (const vect& a) const PURE {
     vect r(*this);
     r&=a;
     return r;
   }
   
-  vect operator| (const vect& a) const {
+  vect operator| (const vect& a) const PURE {
     vect r(*this);
     r|=a;
     return r;
   }
   
-  vect operator^ (const vect& a) const {
+  vect operator^ (const vect& a) const PURE {
     vect r(*this);
     r^=a;
     return r;
   }
   
-  vect<bool,D> operator&& (const vect& a) const {
+  vect<bool,D> operator&& (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]&&a[d];
     return r;
   }
   
-  vect<bool,D> operator|| (const vect& a) const {
+  vect<bool,D> operator|| (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]||a[d];
     return r;
   }
   
-  vect<bool,D> operator== (const vect& a) const {
+  vect<bool,D> operator== (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]==a[d];
     return r;
   }
   
-  vect<bool,D> operator!= (const vect& a) const {
+  vect<bool,D> operator!= (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]!=a[d];
     return r;
   }
   
-  vect<bool,D> operator< (const vect& a) const {
+  vect<bool,D> operator< (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]<a[d];
     return r;
   }
   
-  vect<bool,D> operator<= (const vect& a) const {
+  vect<bool,D> operator<= (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]<=a[d];
     return r;
   }
   
-  vect<bool,D> operator> (const vect& a) const {
+  vect<bool,D> operator> (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]>a[d];
     return r;
   }
   
-  vect<bool,D> operator>= (const vect& a) const {
+  vect<bool,D> operator>= (const vect& a) const PURE {
     vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]>=a[d];
     return r;
@@ -477,7 +508,7 @@ public:
       elements set to either a[i] or b[i], depending on whether
       (*this)[i] is true or not.  */
   template<typename TT>
-  vect<TT,D> ifthen (const vect<TT,D>& a, const vect<TT,D>& b) const {
+  vect<TT,D> ifthen (const vect<TT,D>& a, const vect<TT,D>& b) const /*PURE*/ {
     vect<TT,D> r;
     for (int d=0; d<D; ++d) r[d]=elt[d]?a[d]:b[d];
     return r;
@@ -490,10 +521,10 @@ public:
     vect &vec;
     int d;
   public:
-    iter (vect &a): vec(a), d(0) { }
+    iter (vect &a) PURE: vec(a), d(0) { }
     iter& operator++ () { assert(d<D); ++d; return *this; }
-    bool operator bool () { return d==D; }
-    T& operator* { return vec[d]; }
+    bool operator bool () const PURE { return d==D; }
+    T& operator* () PURE { return vec[d]; }
   };
 #endif
   
@@ -511,6 +542,9 @@ public:
     true or not.  */
 template<typename S,typename T,int D>
 inline vect<T,D> either (const vect<S,D>& a,
+                         const vect<T,D>& b, const vect<T,D>& c) PURE;
+template<typename S,typename T,int D>
+inline vect<T,D> either (const vect<S,D>& a,
                          const vect<T,D>& b, const vect<T,D>& c) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]?b[d]:c[d];
@@ -519,6 +553,8 @@ inline vect<T,D> either (const vect<S,D>& a,
 
 /** Transpose a vector of a vector */
 template<typename T, int D, int DD>
+inline vect<vect<T,D>,DD> xpose (vect<vect<T,DD>,D> const & a) PURE;
+template<typename T, int D, int DD>
 inline vect<vect<T,D>,DD> xpose (vect<vect<T,DD>,D> const & a) {
   vect<vect<T,D>,DD> r;
   for (int dd=0; dd<DD; ++dd) for (int d=0; d<D; ++d) r[dd][d] = a[d][dd];
@@ -527,6 +563,8 @@ inline vect<vect<T,D>,DD> xpose (vect<vect<T,DD>,D> const & a) {
 
 /** Return the element-wise absolute value.  */
 template<typename T,int D>
+inline vect<T,D> abs (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> abs (const vect<T,D>& a) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=abs(a[d]);
@@ -535,6 +573,8 @@ inline vect<T,D> abs (const vect<T,D>& a) {
 
 /** Return the element-wise ceiling.  */
 template<typename T,int D>
+inline vect<T,D> ceil (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> ceil (const vect<T,D>& a) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=ceil(a[d]);
@@ -543,6 +583,8 @@ inline vect<T,D> ceil (const vect<T,D>& a) {
 
 /** Return the element-wise floor.  */
 template<typename T,int D>
+inline vect<T,D> floor (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> floor (const vect<T,D>& a) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=floor(a[d]);
@@ -551,6 +593,8 @@ inline vect<T,D> floor (const vect<T,D>& a) {
 
 /** Return the element-wise maximum of two vectors.  */
 template<typename T,int D>
+inline vect<T,D> max (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> max (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=max(a[d],b[d]);
@@ -559,6 +603,8 @@ inline vect<T,D> max (const vect<T,D>& a, const vect<T,D>& b) {
 
 /** Return the element-wise minimum of two vectors.  */
 template<typename T,int D>
+inline vect<T,D> min (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> min (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=min(a[d],b[d]);
@@ -567,6 +613,8 @@ inline vect<T,D> min (const vect<T,D>& a, const vect<T,D>& b) {
 
 /** Return the element-wise power of two vectors.  */
 template<typename T,typename U,int D>
+inline vect<T,D> pow (const vect<T,D>& a, const vect<U,D>& b) PURE;
+template<typename T,typename U,int D>
 inline vect<T,D> pow (const vect<T,D>& a, const vect<U,D>& b) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=pow(a[d],b[d]);
@@ -575,6 +623,8 @@ inline vect<T,D> pow (const vect<T,D>& a, const vect<U,D>& b) {
 
 /** Return the element-wise integer power of two vectors.  */
 template<typename T,int D>
+inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=ipow(a[d],b[d]);
@@ -582,12 +632,16 @@ inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator+ (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator+ (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r+=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator+ (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator+ (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r+=a;
@@ -595,12 +649,16 @@ inline vect<T,D> operator+ (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator- (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator- (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r-=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator- (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator- (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r-=a;
@@ -608,12 +666,16 @@ inline vect<T,D> operator- (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator* (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator* (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r*=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator* (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator* (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r*=a;
@@ -621,12 +683,16 @@ inline vect<T,D> operator* (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator/ (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator/ (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r/=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator/ (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator/ (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r/=a;
@@ -634,12 +700,16 @@ inline vect<T,D> operator/ (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator% (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator% (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r%=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator% (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator% (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r%=a;
@@ -647,12 +717,16 @@ inline vect<T,D> operator% (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator| (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator| (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r|=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator| (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator| (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r|=a;
@@ -660,12 +734,16 @@ inline vect<T,D> operator| (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator& (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator& (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r&=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator& (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator& (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r&=a;
@@ -673,12 +751,16 @@ inline vect<T,D> operator& (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator^ (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator^ (const vect<T,D>& a, const T& x) {
   vect<T,D> r(a);
   r^=x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator^ (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator^ (const T& x, const vect<T,D>& a) {
   vect<T,D> r(x);
   r^=a;
@@ -686,12 +768,16 @@ inline vect<T,D> operator^ (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator&& (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator&& (const vect<T,D>& a, const T& x) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]&&x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator&& (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator&& (const T& x, const vect<T,D>& a) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=x&&a[d];
@@ -699,12 +785,16 @@ inline vect<T,D> operator&& (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator|| (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<T,D> operator|| (const vect<T,D>& a, const T& x) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]||x;
   return r;
 }
 template<typename T,int D>
+inline vect<T,D> operator|| (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<T,D> operator|| (const T& x, const vect<T,D>& a) {
   vect<T,D> r;
   for (int d=0; d<D; ++d) r[d]=x||a[d];
@@ -712,12 +802,16 @@ inline vect<T,D> operator|| (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator== (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator== (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]==x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator== (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator== (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x==a[d];
@@ -725,12 +819,16 @@ inline vect<bool,D> operator== (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator!= (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator!= (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]!=x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator!= (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator!= (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x!=a[d];
@@ -738,12 +836,16 @@ inline vect<bool,D> operator!= (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator< (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator< (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]<x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator< (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator< (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x<a[d];
@@ -751,12 +853,16 @@ inline vect<bool,D> operator< (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator<= (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator<= (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]<=x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator<= (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator<= (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x<=a[d];
@@ -764,12 +870,16 @@ inline vect<bool,D> operator<= (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator> (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator> (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]>x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator> (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator> (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x>a[d];
@@ -777,12 +887,16 @@ inline vect<bool,D> operator> (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator>= (const vect<T,D>& a, const T& x) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator>= (const vect<T,D>& a, const T& x) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]>=x;
   return r;
 }
 template<typename T,int D>
+inline vect<bool,D> operator>= (const T& x, const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator>= (const T& x, const vect<T,D>& a) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=x>=a[d];
@@ -790,6 +904,8 @@ inline vect<bool,D> operator>= (const T& x, const vect<T,D>& a) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator+ (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator+ (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r+=b;
@@ -797,6 +913,8 @@ inline vect<T,D> operator+ (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator- (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator- (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r-=b;
@@ -804,6 +922,8 @@ inline vect<T,D> operator- (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator* (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator* (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r*=b;
@@ -811,6 +931,8 @@ inline vect<T,D> operator* (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator/ (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator/ (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r/=b;
@@ -818,6 +940,8 @@ inline vect<T,D> operator/ (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator% (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator% (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r%=b;
@@ -825,6 +949,8 @@ inline vect<T,D> operator% (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator& (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator& (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r&=b;
@@ -832,6 +958,8 @@ inline vect<T,D> operator& (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator| (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator| (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r|=b;
@@ -839,6 +967,8 @@ inline vect<T,D> operator| (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<T,D> operator^ (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<T,D> operator^ (const vect<T,D>& a, const vect<T,D>& b) {
   vect<T,D> r(a);
   r^=b;
@@ -846,6 +976,8 @@ inline vect<T,D> operator^ (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator&& (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator&& (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]&&b[d];
@@ -853,6 +985,8 @@ inline vect<bool,D> operator&& (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator|| (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator|| (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]||b[d];
@@ -860,6 +994,8 @@ inline vect<bool,D> operator|| (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator== (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator== (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]==b[d];
@@ -867,6 +1003,8 @@ inline vect<bool,D> operator== (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator!= (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator!= (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]!=b[d];
@@ -874,6 +1012,8 @@ inline vect<bool,D> operator!= (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator< (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator< (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]<b[d];
@@ -881,6 +1021,8 @@ inline vect<bool,D> operator< (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator<= (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator<= (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]<=b[d];
@@ -888,6 +1030,8 @@ inline vect<bool,D> operator<= (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator> (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator> (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]>b[d];
@@ -895,6 +1039,8 @@ inline vect<bool,D> operator> (const vect<T,D>& a, const vect<T,D>& b) {
 }
 
 template<typename T,int D>
+inline vect<bool,D> operator>= (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline vect<bool,D> operator>= (const vect<T,D>& a, const vect<T,D>& b) {
   vect<bool,D> r;
   for (int d=0; d<D; ++d) r[d]=a[d]>=b[d];
@@ -907,6 +1053,8 @@ inline vect<bool,D> operator>= (const vect<T,D>& a, const vect<T,D>& b) {
 
 /** Return true iff any of the elements are true (boolean sum).  */
 template<int D>
+inline bool any (const vect<bool,D>& a) PURE;
+template<int D>
 inline bool any (const vect<bool,D>& a) {
   bool r(false);
   for (int d=0; d<D; ++d) r|=a[d];
@@ -915,6 +1063,8 @@ inline bool any (const vect<bool,D>& a) {
 
 /** Return true iff all of the elements are true (boolean product).  */
 template<int D>
+inline bool all (const vect<bool,D>& a) PURE;
+template<int D>
 inline bool all (const vect<bool,D>& a) {
   bool r(true);
   for (int d=0; d<D; ++d) r&=a[d];
@@ -923,12 +1073,16 @@ inline bool all (const vect<bool,D>& a) {
 
 /** Count the number of elements in the vector.  */
 template<typename T,int D>
+inline int count (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline int count (const vect<T,D>& a) {
   return D;
 }
 
 /** Return the dot product of two vectors.  */
 template<typename T,int D>
+inline T dot (const vect<T,D>& a, const vect<T,D>& b) PURE;
+template<typename T,int D>
 inline T dot (const vect<T,D>& a, const vect<T,D>& b) {
   T r(0);
   for (int d=0; d<D; ++d) r+=a[d]*b[d];
@@ -937,12 +1091,16 @@ inline T dot (const vect<T,D>& a, const vect<T,D>& b) {
 
 /** Return the Euklidean length.  */
 template<typename T,int D>
+inline T hypot (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline T hypot (const vect<T,D>& a) {
   return sqrt(dot(a,a));
 }
 
 /** Return the maximum element.  */
 template<typename T,int D>
+inline T maxval (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline T maxval (const vect<T,D>& a) {
   assert (D>0);
   T r(a[0]);
@@ -952,6 +1110,8 @@ inline T maxval (const vect<T,D>& a) {
 
 /** Return the minimum element.  */
 template<typename T,int D>
+inline T minval (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline T minval (const vect<T,D>& a) {
   assert (D>0);
   T r(a[0]);
@@ -961,6 +1121,8 @@ inline T minval (const vect<T,D>& a) {
 
 /** Return the index of the first maximum element.  */
 template<typename T,int D>
+inline int maxloc (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline int maxloc (const vect<T,D>& a) {
   assert (D>0);
   int r(0);
@@ -970,6 +1132,8 @@ inline int maxloc (const vect<T,D>& a) {
 
 /** Return the index of the first minimum element.  */
 template<typename T,int D>
+inline int minloc (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline int minloc (const vect<T,D>& a) {
   assert (D>0);
   int r(0);
@@ -979,6 +1143,8 @@ inline int minloc (const vect<T,D>& a) {
 
 /** Return the product of the elements.  */
 template<typename T,int D>
+inline T prod (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline T prod (const vect<T,D>& a) {
   T r(1);
   for (int d=0; d<D; ++d) r*=a[d];
@@ -987,12 +1153,16 @@ inline T prod (const vect<T,D>& a) {
 
 /** Return the size (number of elements) of the vector.  */
 template<typename T,int D>
+inline int size (const vect<T,D>& a) CONST;
+template<typename T,int D>
 inline int size (const vect<T,D>& a) {
   return D;
 }
 
 /** Return the sum of the elements.  */
 template<typename T,int D>
+inline T sum (const vect<T,D>& a) PURE;
+template<typename T,int D>
 inline T sum (const vect<T,D>& a) {
   T r(0);
   for (int d=0; d<D; ++d) r+=a[d];
@@ -1098,17 +1268,21 @@ inline ostream& operator<< (ostream& os, const vect<T,D>& a) {
 
 /** Constructor for 2-element vectors from 2 elements.  */
 template<typename T>
+inline vect<T,2>::vect<T,2> (const T& x, const T& y) PURE;
+template<typename T>
 inline vect<T,2>::vect<T,2> (const T& x, const T& y) {
   elt[0]=x; elt[1]=y;
 }
 
 /** Constructor for 3-element vectors from 3 elements.  */
+vect (const T& x, const T& y, const T& z) PURE;
 vect (const T& x, const T& y, const T& z) {
   assert (D==3);
   elt[0]=x; elt[1]=y; elt[2]=z;
 }
 
 /** Constructor for 4-element vectors from 4 elements.  */
+vect (const T& x, const T& y, const T& z, const T& t) PURE;
 vect (const T& x, const T& y, const T& z, const T& t) {
   assert (D==4);
   elt[0]=x; elt[1]=y; elt[2]=z; elt[3]=t;