Vector.h

Go to the documentation of this file.
00001 
00007 #ifndef VECTOR_H
00008 #define VECTOR_H
00009 
00010 #include <math.h>
00011 #include <stdio.h>
00012 #include "common.h"
00013 
00014 class Vector;
00015 
00016 class FloatVector {
00017  public:
00018   float x,y,z;
00019   inline FloatVector(void) { ; }
00020   inline FloatVector(const Vector &v);
00021 };
00022 
00023 #ifdef ARCH_POWERPC
00024 
00025 #include "builtins.h"
00026 
00027 inline double namd_rsqrt(double x)
00028 {
00029   double r0, r1, r2;
00030 
00031   /*------------------------------------------*/
00032   /* use reciprocal sqrt estimate instruction */
00033   /*------------------------------------------*/
00034   r0 = __frsqrte(x);
00035 
00036   /*----------------------*/
00037   /* 1st newton iteration */
00038   /*----------------------*/
00039   r1 = r0 + 0.5*r0*(1.0 - x*r0*r0);
00040 
00041   /*----------------------*/
00042   /* 2nd newton iteration */
00043   /*----------------------*/
00044   r2 = r1 + 0.5*r1*(1.0 - x*r1*r1);
00045 
00046   return r2;
00047 }
00048 
00049 inline double namd_reciprocal (double x) {
00050   register double rx;
00051   
00052   rx = __fres(x);             // 0th estimate (13 bits)
00053   rx = rx + rx*(1.0 - x*rx);  // 1st Newton iteration (26 bits)
00054   rx = rx + rx*(1.0 - x*rx);  // 2nd Newton iteration (52 bits = full mantissa for a double)
00055 
00056   return rx;
00057 }
00058 
00059 #else
00060 #define namd_rsqrt(x)        (1.0 / sqrt (x))
00061 #define namd_reciprocal(x)   (1.0 / x)
00062 #endif
00063 
00064 class Vector {
00065    public:
00066      BigReal x,y,z;
00067      
00068      inline Vector(void) : x(-99999), y(-99999), z(-99999) { ; }
00069 //     inline Vector(void) { ; }
00070      
00071      inline Vector( BigReal newx, BigReal newy, BigReal newz)
00072        : x(newx), y(newy), z(newz) { ; }
00073 
00074      inline Vector( BigReal newv )  // allow Vector v = 0; etc.
00075        : x(newv), y(newv), z(newv) { ; }
00076 
00077      inline Vector(const FloatVector &v) : x(v.x), y(v.y), z(v.z) { ; }
00078 
00079      inline BigReal &operator[](int i) {
00080        return i==0 ? x
00081              :i==1 ? y
00082              :i==2 ? z
00083              :(NAMD_die("Vector reference out of bounds."), x);
00084 
00085      }
00086 
00087      //  v1 = const;
00088      inline Vector& operator=(const BigReal &v2) {
00089        x = v2;
00090        y = v2;
00091        z = v2;
00092        return *this;
00093      }
00094 
00095      //  v1 += v2;
00096      inline void operator+=(const Vector &v2) {
00097        x += v2.x;
00098        y += v2.y;
00099        z += v2.z;
00100      }
00101 
00102      // v1 -= v2;
00103      inline void operator-=(const Vector &v2) {
00104        x -= v2.x;
00105        y -= v2.y;
00106        z -= v2.z;
00107      }
00108 
00109      // v1 *= const
00110      inline void operator*=(const BigReal &v2) {
00111        x *= v2;
00112        y *= v2;
00113        z *= v2;
00114      }
00115 
00116      // v1 /= const
00117      inline void operator/=(const BigReal& v2) {
00118        BigReal v2_recip = namd_reciprocal(v2);
00119        x *= v2_recip;
00120        y *= v2_recip;
00121        z *= v2_recip;
00122      }
00123 
00124      inline friend int operator == (const Vector& v1, const Vector& v2) {
00125        return v1.x == v2.x && v1.y == v2.y && v1.z == v2.z;
00126      }
00127      inline friend int operator != (const Vector& v1, const Vector& v2) {
00128        // return !(v1.x == v2.x && v1.y == v2.y && v1.z == v2.z);
00129        return v1.x != v2.x || v1.y != v2.y || v1.z != v2.z;
00130      }
00131 
00132      // addition of two vectors
00133      inline friend Vector operator+(const Vector& v1, const Vector& v2) {
00134        return Vector( v1.x+v2.x, v1.y+v2.y, v1.z+v2.z);
00135      }
00136 
00137      // negation
00138      inline friend Vector operator-(const Vector &v1) {
00139        return Vector( -v1.x, -v1.y, -v1.z);
00140      }
00141 
00142      // subtraction
00143      inline friend Vector operator-(const Vector &v1, const Vector &v2) {
00144        return Vector( v1.x-v2.x, v1.y-v2.y, v1.z-v2.z);
00145      }
00146      // inner ("dot") product
00147      inline friend BigReal operator*(const Vector &v1, const Vector &v2) {
00148        return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
00149      }
00150      // scalar product
00151      inline friend Vector operator*(const BigReal &f, const Vector &v1) {
00152        return Vector(f*v1.x, f*v1.y, f*v1.z);
00153      }
00154      // scalar product
00155      inline friend Vector operator*(const Vector &v1, const BigReal &f) {
00156        return Vector(f*v1.x, f*v1.y, f*v1.z);
00157      }
00158      // division by a scalar
00159      inline friend Vector operator/(const Vector &v1, const BigReal &f) {
00160 //       if (!f)
00161 //         NAMD_die("Division by 0 on a vector operation.");
00162        return Vector(v1.x/f, v1.y/f, v1.z/f);
00163      }
00164      
00165      // return the norm
00166      inline BigReal length(void) const {
00167        return sqrt(x*x+y*y+z*z);
00168      }
00169      
00170      inline BigReal length2(void) const {
00171        return (x*x + y*y + z*z);
00172      }
00173 
00174      inline BigReal rlength (void) {
00175        return namd_rsqrt (x*x + y*y + z*z);
00176      }
00177 
00178      // return the unit vector in the same direction
00179      inline Vector unit(void) const {
00180        return Vector(x, y, z)/length();
00181      }
00182      
00183      
00184      // one cross product  v3 = cross(v1, v2)
00185      inline friend Vector cross(const Vector &v1, const Vector &v2) {
00186        return Vector( v1.y*v2.z-v2.y*v1.z,
00187                      // -v1.x*v2.z+v2.x*v1.z,
00188                       v2.x*v1.z-v1.x*v2.z,
00189                       v1.x*v2.y-v2.x*v1.y  );
00190      }
00191 
00192      // multiplying a cross product by a scalar is very common
00193      // one cross product  v3 = k*cross(v1, v2)
00194      inline friend Vector cross(const Real &k, const Vector &v1, const Vector &v2) {
00195        return Vector( k*(v1.y*v2.z-v2.y*v1.z),
00196                       // k*(-v1.x*v2.z+v2.x*v1.z),
00197                       k*(v2.x*v1.z-v1.x*v2.z),
00198                       k*(v1.x*v2.y-v2.x*v1.y) );
00199      }
00200 
00201      inline friend Vector cross(const BigReal &k, const Vector &v1, const Vector &v2) {
00202        return Vector( k*(v1.y*v2.z-v2.y*v1.z),
00203                       // k*(-v1.x*v2.z+v2.x*v1.z),
00204                       k*(v2.x*v1.z-v1.x*v2.z),
00205                       k*(v1.x*v2.y-v2.x*v1.y) );
00206      }
00207 
00208      // A = A x B  -- why do you want this function, anyway?
00209      void cross(const Vector &v2) {
00210        BigReal xx =  y*v2.z-v2.y*z;
00211        // BigReal yy = -x*v2.z+v2.x*z;
00212        BigReal yy = v2.x*z-x*v2.z;
00213        z =  x*v2.y-v2.x*y;
00214        y=yy;
00215        x=xx;
00216      }
00217 
00218      // returns (*this) * V2
00219      BigReal dot(const Vector &v2) {
00220        return x*v2.x + y*v2.y + z*v2.z;
00221      }
00222 
00223      // set the vector based on a string.  If bad, return FALSE
00224      // the string can be in the form "x y z" or "x, y, z"
00225      Bool set(const char *s) {
00226         double a[3];    // read into doubles, since I don't know what
00227         char tmp[100];  // a "BigReal" is in real life
00228         // cheap way to get commas, etc.  a poor regex
00229        int i=sscanf(s, "%lf%99[ \t,]%lf%99[ \t,]%lf%99s",
00230                     a, tmp, a+1, tmp, a+2, tmp);
00231        if (i != 5) return FALSE;
00232        const char *t = s;       // now count commas (for "1,,,,2,  , 3")
00233        int flg = 0;                 // and check for "1 2,,3"
00234        i = 0;
00235        for (;*t;t++) {
00236           if (*t == ',') { 
00237              if (flg == 0) {   // expecting non-whitespace
00238                 return FALSE;  //    so error
00239              }
00240              flg = 0;          // now expect non-whitespace
00241              i++;              // and increment comma counter
00242           }
00243           else if (*t != ' ' && *t != '\t') {  // got non-whitespace
00244              flg = 1;          // so the next can be whitespace or commas
00245           }
00246        }
00247        if (i == 0 || i == 2) {  // allow "1 2 3" or "1, 2,3" forms
00248           x = a[0]; y = a[1]; z = a[2];
00249           return TRUE;
00250        }
00251        return FALSE;
00252      }
00253 };
00254 
00255 class zVector : public Vector {
00256   public:
00257      inline zVector(void) : Vector(0,0,0) { ; }
00258      inline zVector(const Vector &v) : Vector(v) { ; }
00259 };
00260 
00261 
00262 class AlignVector : public Vector {
00263  public:
00264   BigReal pad;
00265   inline AlignVector(void) : Vector(0,0,0) { pad = 0.0; }
00266   inline AlignVector(const Vector &v) : Vector(v) { pad = 0.0; }
00267 
00268   inline AlignVector( BigReal newx, BigReal newy, BigReal newz)
00269     : Vector (newx, newy, newz) { pad = 0.0; }
00270   
00271   inline AlignVector( BigReal newv )  // allow Vector v = 0; etc.
00272     : Vector (newv) { pad = 0.0; }
00273   
00274   inline AlignVector(const FloatVector &v) : Vector (v) { pad = 0.0; }
00275 };
00276 
00277 
00278 inline FloatVector::FloatVector(const Vector &v) : x(v.x), y(v.y), z(v.z) { ; }
00279 
00280 //#define TEST_VECTOR_CLASS
00281 #ifdef TEST_VECTOR_CLASS
00282 main()
00283 {
00284   Vector v1(1.1,2.2, 3.3);
00285   Vector v2(-1, 55, 32.1);
00286   Vector v3(v1+2*v2);
00287   Vector v4;
00288   std::cout << v1 << "  " << v2 << "  " << v3 << "  " << v4 << '\n';
00289   std::cout << v1*v2 << "  "  << v3-v1-2*v2 <<"  "<< v2 * v3 <<"  "<< v3*v2 <<'\n';
00290   v4 = v3*5 - v2/4;
00291   std::cout << v4 << "  " << v3*5.0 - v2/4.0 << '\n';
00292   std::cout << v4[0] << "  "  << v4[1] << "  " << v4[2] << '\n';
00293 //  std::cout.flush();
00294 //  std::cout << v4[3];
00295   std::cout << cross(v1, v2) << '\n';
00296   std::cout << v1 << '\n';  
00297   v1 += v2;
00298   std::cout << v1 << '\n';
00299   v1 -= v2;
00300   std::cout << v1 << '\n';
00301   {
00302     Vector v1(1.0, 2.0, 3.0);  // some more examples, but I was too lazy to
00303     Vector v2 = v1.unit();     // fix the names
00304     std::cout << v2 << '\n';
00305     std::cout << v2.dot(v1) << '\n';
00306     std::cout << v1.length() << '\n';
00307     v1 *= -1;
00308     v1 += v2*14;
00309     std::cout << v1 << '\n';
00310   }
00311 }
00312 #endif
00313 
00314 #endif
00315 

Generated on Thu Nov 23 01:17:15 2017 for NAMD by  doxygen 1.4.7