text: // Copyright (c)2008-2011, Preferred Infrastructure Inc.
  // 
  // All rights reserved.
  // 
  // Redistribution and use in source and binary forms, with or without
  // modification, are permitted provided that the following conditions are met:
  // 
  //     * Redistributions of source code must retain the above copyright
  //       notice, this list of conditions and the following disclaimer.
  // 
  //     * Redistributions in binary form must reproduce the above
  //       copyright notice, this list of conditions and the following
  //       disclaimer in the documentation and/or other materials provided
  //       with the distribution.
  // 
  //     * Neither the name of Preferred Infrastructure nor the names of other
  //       contributors may be used to endorse or promote products derived
  //       from this software without specific prior written permission.
  // 
  // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  
  ///generated by vector_h.rb
  
  #ifndef INCLUDE_GUARD_PFI_MATH_VECTOR_H_
  #define INCLUDE_GUARD_PFI_MATH_VECTOR_H_
  #include <cassert>
  #include <cmath>
  #include <iostream>
  
  namespace pfi{
  namespace math{
  namespace vector{
  namespace component_by_name{    
    ///linear-algebraic vector of type t, dimension 2.
    template<class t> class vector2{
    public:
      ///member
      t  x,y;
      ///zero constructor
      vector2():x(0),y(0){}
      ///constructor from components
      vector2(t x,t y):x(x),y(y){}
      ///constructor from a pointer
      vector2(t *v):x(v[0]),y(v[1]){}
    
      ///default normalization;
      vector2 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=x;v[1]=y;;
      }
  
      ///component access operator.
      t& operator[](int index){
        assert(index >= 0 && index < 2);
        if (index == 0)
          return x;
        else
          return y;
      }
    };
  
    ///vector addition
    template<class t>
      vector2<t> operator+(const vector2<t> & a, const vector2<t> & b){
      return vector2<t>(a.x + b.x,a.y + b.y);
    }
  
    ///vector subtraction
    template<class t>
      vector2<t> operator-(const vector2<t> & a, const vector2<t> & b){
      return vector2<t>(a.x - b.x,a.y - b.y);
    }
  
    ///vector negate
    template<class t>
      vector2<t> operator-(const vector2<t> & a){
      return vector2<t>(-a.x,-a.y);
    }
  
    ///scalar multiplication
    template<class t>
      vector2<t> operator*(const t &m, const vector2<t> & a){
      return vector2<t>(m*a.x,m*a.y);
    }
    template<class t>
      vector2<t> operator*(const vector2<t> & a, const t &m){
      return vector2<t>(m*a.x,m*a.y);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector2<t> operator/(const vector2<t> & a, const t &m){
      return vector2<t>(a.x/m,a.y/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector2<t> &a, const vector2<t> & b){
      return a.x*b.x + a.y*b.y;
    }
  
    ///outer product; outer product uses % operator.
    template<class t>
      t operator%(const vector2<t> & a, const vector2<t> & b){
      return a.x*b.y - a.y*b.x;
    }
  
    /// +=
    template<class t>
      vector2<t>& operator+=(vector2<t> &a, const vector2<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector2<t>& operator-=(vector2<t> &a, const vector2<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector2<t>& operator*=(vector2<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector2<t>& operator/=(vector2<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector2<t> & a, const vector2<t> & b){
      return a.x==b.x && a.y==b.y;
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector2<t> & a, const vector2<t> & b){
      if(a.x != b.x) return a.x < b.x;
      return a.y < b.y;
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector2<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector2<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector2<t> normalize(const vector2<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector2<t> & a, const vector2<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector2<t> vector2<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector2<t>& a){
      t tmp0,tmp1; 
      cin  >> tmp0 >> tmp1;
      a = vector2<t>(tmp0,tmp1);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector2<t> &a){
      return cout << a.x<< " " << a.y;
    }
  
  }//component_by_name
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 2.
    template<class t> class vector2{
    public:
      ///member
      t  r[2];
      ///zero constructor
      vector2(){
        r[0]=0; r[1]=0; 
      }
      ///constructor from components
      vector2(t r0,t r1){
        r[0]=r0; r[1]=r1; 
      }
      ///constructor from a pointer
      vector2(t *v){
        r[0]=v[0]; r[1]=v[1]; 
      }
    
      ///default normalization;
      vector2 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector2<t> operator+(const vector2<t> & a, const vector2<t> & b){
      return vector2<t>(a.r[0] + b.r[0],a.r[1] + b.r[1]);
    }
  
    ///vector subtraction
    template<class t>
      vector2<t> operator-(const vector2<t> & a, const vector2<t> & b){
      return vector2<t>(a.r[0] - b.r[0],a.r[1] - b.r[1]);
    }
  
    ///vector negate
    template<class t>
      vector2<t> operator-(const vector2<t> & a){
      return vector2<t>(-a.r[0],-a.r[1]);
    }
  
    ///scalar multiplication
    template<class t>
      vector2<t> operator*(const t &m, const vector2<t> & a){
      return vector2<t>(m*a.r[0],m*a.r[1]);
    }
    template<class t>
      vector2<t> operator*(const vector2<t> & a, const t &m){
      return vector2<t>(m*a.r[0],m*a.r[1]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector2<t> operator/(const vector2<t> & a, const t &m){
      return vector2<t>(a.r[0]/m,a.r[1]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector2<t> &a, const vector2<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1];
    }
  
    ///outer product; outer product uses % operator.
    template<class t>
      t operator%(const vector2<t> & a, const vector2<t> & b){
      return a.x*b.y - a.y*b.x;
    }
  
    /// +=
    template<class t>
      vector2<t>& operator+=(vector2<t> &a, const vector2<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector2<t>& operator-=(vector2<t> &a, const vector2<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector2<t>& operator*=(vector2<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector2<t>& operator/=(vector2<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector2<t> & a, const vector2<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector2<t> & a, const vector2<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      return a.r[1] < b.r[1];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector2<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector2<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector2<t> normalize(const vector2<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector2<t> & a, const vector2<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector2<t> vector2<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector2<t>& a){
      t tmp0,tmp1; 
      cin  >> tmp0 >> tmp1;
      a = vector2<t>(tmp0,tmp1);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector2<t> &a){
      return cout << a.r[0]<< " " << a.r[1];
    }
  
  }//component_by_array
  
  
  namespace component_by_name{    
    ///linear-algebraic vector of type t, dimension 3.
    template<class t> class vector3{
    public:
      ///member
      t  x,y,z;
      ///zero constructor
      vector3():x(0),y(0),z(0){}
      ///constructor from components
      vector3(t x,t y,t z):x(x),y(y),z(z){}
      ///constructor from a pointer
      vector3(t *v):x(v[0]),y(v[1]),z(v[2]){}
    
      ///default normalization;
      vector3 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=x;v[1]=y;v[2]=z;;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return x;
          case 1: return y;
          case 2: return z;
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector3<t> operator+(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.x + b.x,a.y + b.y,a.z + b.z);
    }
  
    ///vector subtraction
    template<class t>
      vector3<t> operator-(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.x - b.x,a.y - b.y,a.z - b.z);
    }
  
    ///vector negate
    template<class t>
      vector3<t> operator-(const vector3<t> & a){
      return vector3<t>(-a.x,-a.y,-a.z);
    }
  
    ///scalar multiplication
    template<class t>
      vector3<t> operator*(const t &m, const vector3<t> & a){
      return vector3<t>(m*a.x,m*a.y,m*a.z);
    }
    template<class t>
      vector3<t> operator*(const vector3<t> & a, const t &m){
      return vector3<t>(m*a.x,m*a.y,m*a.z);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector3<t> operator/(const vector3<t> & a, const t &m){
      return vector3<t>(a.x/m,a.y/m,a.z/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector3<t> &a, const vector3<t> & b){
      return a.x*b.x + a.y*b.y + a.z*b.z;
    }
  
    ///outer product; outer product uses % operator.
    template<class t>
      vector3<t> operator%(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.y*b.z - a.z*b.y,a.z*b.x - a.x*b.z,a.x*b.y - a.y*b.x);
    }
  
    /// +=
    template<class t>
      vector3<t>& operator+=(vector3<t> &a, const vector3<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector3<t>& operator-=(vector3<t> &a, const vector3<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector3<t>& operator*=(vector3<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector3<t>& operator/=(vector3<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector3<t> & a, const vector3<t> & b){
      return a.x==b.x && a.y==b.y && a.z==b.z;
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector3<t> & a, const vector3<t> & b){
      if(a.x != b.x) return a.x < b.x;
      if(a.y != b.y) return a.y < b.y;
      return a.z < b.z;
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector3<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector3<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector3<t> normalize(const vector3<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector3<t> & a, const vector3<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector3<t> vector3<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector3<t>& a){
      t tmp0,tmp1,tmp2; 
      cin  >> tmp0 >> tmp1 >> tmp2;
      a = vector3<t>(tmp0,tmp1,tmp2);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector3<t> &a){
      return cout << a.x<< " " << a.y<< " " << a.z;
    }
  
  }//component_by_name
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 3.
    template<class t> class vector3{
    public:
      ///member
      t  r[3];
      ///zero constructor
      vector3(){
        r[0]=0; r[1]=0; r[2]=0; 
      }
      ///constructor from components
      vector3(t r0,t r1,t r2){
        r[0]=r0; r[1]=r1; r[2]=r2; 
      }
      ///constructor from a pointer
      vector3(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; 
      }
    
      ///default normalization;
      vector3 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector3<t> operator+(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2]);
    }
  
    ///vector subtraction
    template<class t>
      vector3<t> operator-(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2]);
    }
  
    ///vector negate
    template<class t>
      vector3<t> operator-(const vector3<t> & a){
      return vector3<t>(-a.r[0],-a.r[1],-a.r[2]);
    }
  
    ///scalar multiplication
    template<class t>
      vector3<t> operator*(const t &m, const vector3<t> & a){
      return vector3<t>(m*a.r[0],m*a.r[1],m*a.r[2]);
    }
    template<class t>
      vector3<t> operator*(const vector3<t> & a, const t &m){
      return vector3<t>(m*a.r[0],m*a.r[1],m*a.r[2]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector3<t> operator/(const vector3<t> & a, const t &m){
      return vector3<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector3<t> &a, const vector3<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2];
    }
  
    ///outer product; outer product uses % operator.
    template<class t>
      vector3<t> operator%(const vector3<t> & a, const vector3<t> & b){
      return vector3<t>(a.y*b.z - a.z*b.y,a.z*b.x - a.x*b.z,a.x*b.y - a.y*b.x);
    }
  
    /// +=
    template<class t>
      vector3<t>& operator+=(vector3<t> &a, const vector3<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector3<t>& operator-=(vector3<t> &a, const vector3<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector3<t>& operator*=(vector3<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector3<t>& operator/=(vector3<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector3<t> & a, const vector3<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector3<t> & a, const vector3<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      return a.r[2] < b.r[2];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector3<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector3<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector3<t> normalize(const vector3<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector3<t> & a, const vector3<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector3<t> vector3<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector3<t>& a){
      t tmp0,tmp1,tmp2; 
      cin  >> tmp0 >> tmp1 >> tmp2;
      a = vector3<t>(tmp0,tmp1,tmp2);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector3<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2];
    }
  
  }//component_by_array
  
  
  namespace component_by_name{    
    ///linear-algebraic vector of type t, dimension 4.
    template<class t> class vector4{
    public:
      ///member
      t  x,y,z,w;
      ///zero constructor
      vector4():x(0),y(0),z(0),w(0){}
      ///constructor from components
      vector4(t x,t y,t z,t w):x(x),y(y),z(z),w(w){}
      ///constructor from a pointer
      vector4(t *v):x(v[0]),y(v[1]),z(v[2]),w(v[3]){}
    
      ///default normalization;
      vector4 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=x;v[1]=y;v[2]=z;v[3]=w;;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return x;
          case 1: return y;
          case 2: return z;
          case 3: return w;
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector4<t> operator+(const vector4<t> & a, const vector4<t> & b){
      return vector4<t>(a.x + b.x,a.y + b.y,a.z + b.z,a.w + b.w);
    }
  
    ///vector subtraction
    template<class t>
      vector4<t> operator-(const vector4<t> & a, const vector4<t> & b){
      return vector4<t>(a.x - b.x,a.y - b.y,a.z - b.z,a.w - b.w);
    }
  
    ///vector negate
    template<class t>
      vector4<t> operator-(const vector4<t> & a){
      return vector4<t>(-a.x,-a.y,-a.z,-a.w);
    }
  
    ///scalar multiplication
    template<class t>
      vector4<t> operator*(const t &m, const vector4<t> & a){
      return vector4<t>(m*a.x,m*a.y,m*a.z,m*a.w);
    }
    template<class t>
      vector4<t> operator*(const vector4<t> & a, const t &m){
      return vector4<t>(m*a.x,m*a.y,m*a.z,m*a.w);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector4<t> operator/(const vector4<t> & a, const t &m){
      return vector4<t>(a.x/m,a.y/m,a.z/m,a.w/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector4<t> &a, const vector4<t> & b){
      return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
    }
  
  
    /// +=
    template<class t>
      vector4<t>& operator+=(vector4<t> &a, const vector4<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector4<t>& operator-=(vector4<t> &a, const vector4<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector4<t>& operator*=(vector4<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector4<t>& operator/=(vector4<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector4<t> & a, const vector4<t> & b){
      return a.x==b.x && a.y==b.y && a.z==b.z && a.w==b.w;
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector4<t> & a, const vector4<t> & b){
      if(a.x != b.x) return a.x < b.x;
      if(a.y != b.y) return a.y < b.y;
      if(a.z != b.z) return a.z < b.z;
      return a.w < b.w;
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector4<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector4<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector4<t> normalize(const vector4<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector4<t> & a, const vector4<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector4<t> vector4<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector4<t>& a){
      t tmp0,tmp1,tmp2,tmp3; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3;
      a = vector4<t>(tmp0,tmp1,tmp2,tmp3);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector4<t> &a){
      return cout << a.x<< " " << a.y<< " " << a.z<< " " << a.w;
    }
  
  }//component_by_name
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 4.
    template<class t> class vector4{
    public:
      ///member
      t  r[4];
      ///zero constructor
      vector4(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; 
      }
      ///constructor from components
      vector4(t r0,t r1,t r2,t r3){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; 
      }
      ///constructor from a pointer
      vector4(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; 
      }
    
      ///default normalization;
      vector4 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector4<t> operator+(const vector4<t> & a, const vector4<t> & b){
      return vector4<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3]);
    }
  
    ///vector subtraction
    template<class t>
      vector4<t> operator-(const vector4<t> & a, const vector4<t> & b){
      return vector4<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3]);
    }
  
    ///vector negate
    template<class t>
      vector4<t> operator-(const vector4<t> & a){
      return vector4<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3]);
    }
  
    ///scalar multiplication
    template<class t>
      vector4<t> operator*(const t &m, const vector4<t> & a){
      return vector4<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3]);
    }
    template<class t>
      vector4<t> operator*(const vector4<t> & a, const t &m){
      return vector4<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector4<t> operator/(const vector4<t> & a, const t &m){
      return vector4<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector4<t> &a, const vector4<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3];
    }
  
  
    /// +=
    template<class t>
      vector4<t>& operator+=(vector4<t> &a, const vector4<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector4<t>& operator-=(vector4<t> &a, const vector4<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector4<t>& operator*=(vector4<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector4<t>& operator/=(vector4<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector4<t> & a, const vector4<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector4<t> & a, const vector4<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      return a.r[3] < b.r[3];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector4<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector4<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector4<t> normalize(const vector4<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector4<t> & a, const vector4<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector4<t> vector4<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector4<t>& a){
      t tmp0,tmp1,tmp2,tmp3; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3;
      a = vector4<t>(tmp0,tmp1,tmp2,tmp3);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector4<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 5.
    template<class t> class vector5{
    public:
      ///member
      t  r[5];
      ///zero constructor
      vector5(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; 
      }
      ///constructor from components
      vector5(t r0,t r1,t r2,t r3,t r4){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; 
      }
      ///constructor from a pointer
      vector5(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; 
      }
    
      ///default normalization;
      vector5 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector5<t> operator+(const vector5<t> & a, const vector5<t> & b){
      return vector5<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4]);
    }
  
    ///vector subtraction
    template<class t>
      vector5<t> operator-(const vector5<t> & a, const vector5<t> & b){
      return vector5<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4]);
    }
  
    ///vector negate
    template<class t>
      vector5<t> operator-(const vector5<t> & a){
      return vector5<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4]);
    }
  
    ///scalar multiplication
    template<class t>
      vector5<t> operator*(const t &m, const vector5<t> & a){
      return vector5<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4]);
    }
    template<class t>
      vector5<t> operator*(const vector5<t> & a, const t &m){
      return vector5<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector5<t> operator/(const vector5<t> & a, const t &m){
      return vector5<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector5<t> &a, const vector5<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4];
    }
  
  
    /// +=
    template<class t>
      vector5<t>& operator+=(vector5<t> &a, const vector5<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector5<t>& operator-=(vector5<t> &a, const vector5<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector5<t>& operator*=(vector5<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector5<t>& operator/=(vector5<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector5<t> & a, const vector5<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector5<t> & a, const vector5<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      return a.r[4] < b.r[4];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector5<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector5<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector5<t> normalize(const vector5<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector5<t> & a, const vector5<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector5<t> vector5<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector5<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4;
      a = vector5<t>(tmp0,tmp1,tmp2,tmp3,tmp4);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector5<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 6.
    template<class t> class vector6{
    public:
      ///member
      t  r[6];
      ///zero constructor
      vector6(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; r[5]=0; 
      }
      ///constructor from components
      vector6(t r0,t r1,t r2,t r3,t r4,t r5){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; r[5]=r5; 
      }
      ///constructor from a pointer
      vector6(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; r[5]=v[5]; 
      }
    
      ///default normalization;
      vector6 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];v[5]=r[5];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
          case 5: return r[5];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector6<t> operator+(const vector6<t> & a, const vector6<t> & b){
      return vector6<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4],a.r[5] + b.r[5]);
    }
  
    ///vector subtraction
    template<class t>
      vector6<t> operator-(const vector6<t> & a, const vector6<t> & b){
      return vector6<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4],a.r[5] - b.r[5]);
    }
  
    ///vector negate
    template<class t>
      vector6<t> operator-(const vector6<t> & a){
      return vector6<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4],-a.r[5]);
    }
  
    ///scalar multiplication
    template<class t>
      vector6<t> operator*(const t &m, const vector6<t> & a){
      return vector6<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5]);
    }
    template<class t>
      vector6<t> operator*(const vector6<t> & a, const t &m){
      return vector6<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector6<t> operator/(const vector6<t> & a, const t &m){
      return vector6<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m,a.r[5]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector6<t> &a, const vector6<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4] + a.r[5]*b.r[5];
    }
  
  
    /// +=
    template<class t>
      vector6<t>& operator+=(vector6<t> &a, const vector6<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector6<t>& operator-=(vector6<t> &a, const vector6<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector6<t>& operator*=(vector6<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector6<t>& operator/=(vector6<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector6<t> & a, const vector6<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4] && a.r[5]==b.r[5];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector6<t> & a, const vector6<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      if(a.r[4] != b.r[4]) return a.r[4] < b.r[4];
      return a.r[5] < b.r[5];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector6<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector6<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector6<t> normalize(const vector6<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector6<t> & a, const vector6<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector6<t> vector6<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector6<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4,tmp5; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5;
      a = vector6<t>(tmp0,tmp1,tmp2,tmp3,tmp4,tmp5);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector6<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4]<< " " << a.r[5];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 7.
    template<class t> class vector7{
    public:
      ///member
      t  r[7];
      ///zero constructor
      vector7(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; r[5]=0; r[6]=0; 
      }
      ///constructor from components
      vector7(t r0,t r1,t r2,t r3,t r4,t r5,t r6){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; r[5]=r5; r[6]=r6; 
      }
      ///constructor from a pointer
      vector7(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; r[5]=v[5]; r[6]=v[6]; 
      }
    
      ///default normalization;
      vector7 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];v[5]=r[5];v[6]=r[6];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
          case 5: return r[5];
          case 6: return r[6];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector7<t> operator+(const vector7<t> & a, const vector7<t> & b){
      return vector7<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4],a.r[5] + b.r[5],a.r[6] + b.r[6]);
    }
  
    ///vector subtraction
    template<class t>
      vector7<t> operator-(const vector7<t> & a, const vector7<t> & b){
      return vector7<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4],a.r[5] - b.r[5],a.r[6] - b.r[6]);
    }
  
    ///vector negate
    template<class t>
      vector7<t> operator-(const vector7<t> & a){
      return vector7<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4],-a.r[5],-a.r[6]);
    }
  
    ///scalar multiplication
    template<class t>
      vector7<t> operator*(const t &m, const vector7<t> & a){
      return vector7<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6]);
    }
    template<class t>
      vector7<t> operator*(const vector7<t> & a, const t &m){
      return vector7<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector7<t> operator/(const vector7<t> & a, const t &m){
      return vector7<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m,a.r[5]/m,a.r[6]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector7<t> &a, const vector7<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4] + a.r[5]*b.r[5] + a.r[6]*b.r[6];
    }
  
  
    /// +=
    template<class t>
      vector7<t>& operator+=(vector7<t> &a, const vector7<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector7<t>& operator-=(vector7<t> &a, const vector7<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector7<t>& operator*=(vector7<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector7<t>& operator/=(vector7<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector7<t> & a, const vector7<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4] && a.r[5]==b.r[5] && a.r[6]==b.r[6];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector7<t> & a, const vector7<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      if(a.r[4] != b.r[4]) return a.r[4] < b.r[4];
      if(a.r[5] != b.r[5]) return a.r[5] < b.r[5];
      return a.r[6] < b.r[6];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector7<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector7<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector7<t> normalize(const vector7<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector7<t> & a, const vector7<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector7<t> vector7<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector7<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5 >> tmp6;
      a = vector7<t>(tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector7<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4]<< " " << a.r[5]<< " " << a.r[6];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 8.
    template<class t> class vector8{
    public:
      ///member
      t  r[8];
      ///zero constructor
      vector8(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; r[5]=0; r[6]=0; r[7]=0; 
      }
      ///constructor from components
      vector8(t r0,t r1,t r2,t r3,t r4,t r5,t r6,t r7){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; r[5]=r5; r[6]=r6; r[7]=r7; 
      }
      ///constructor from a pointer
      vector8(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; r[5]=v[5]; r[6]=v[6]; r[7]=v[7]; 
      }
    
      ///default normalization;
      vector8 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];v[5]=r[5];v[6]=r[6];v[7]=r[7];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
          case 5: return r[5];
          case 6: return r[6];
          case 7: return r[7];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector8<t> operator+(const vector8<t> & a, const vector8<t> & b){
      return vector8<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4],a.r[5] + b.r[5],a.r[6] + b.r[6],a.r[7] + b.r[7]);
    }
  
    ///vector subtraction
    template<class t>
      vector8<t> operator-(const vector8<t> & a, const vector8<t> & b){
      return vector8<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4],a.r[5] - b.r[5],a.r[6] - b.r[6],a.r[7] - b.r[7]);
    }
  
    ///vector negate
    template<class t>
      vector8<t> operator-(const vector8<t> & a){
      return vector8<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4],-a.r[5],-a.r[6],-a.r[7]);
    }
  
    ///scalar multiplication
    template<class t>
      vector8<t> operator*(const t &m, const vector8<t> & a){
      return vector8<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7]);
    }
    template<class t>
      vector8<t> operator*(const vector8<t> & a, const t &m){
      return vector8<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector8<t> operator/(const vector8<t> & a, const t &m){
      return vector8<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m,a.r[5]/m,a.r[6]/m,a.r[7]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector8<t> &a, const vector8<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4] + a.r[5]*b.r[5] + a.r[6]*b.r[6] + a.r[7]*b.r[7];
    }
  
  
    /// +=
    template<class t>
      vector8<t>& operator+=(vector8<t> &a, const vector8<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector8<t>& operator-=(vector8<t> &a, const vector8<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector8<t>& operator*=(vector8<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector8<t>& operator/=(vector8<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector8<t> & a, const vector8<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4] && a.r[5]==b.r[5] && a.r[6]==b.r[6] && a.r[7]==b.r[7];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector8<t> & a, const vector8<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      if(a.r[4] != b.r[4]) return a.r[4] < b.r[4];
      if(a.r[5] != b.r[5]) return a.r[5] < b.r[5];
      if(a.r[6] != b.r[6]) return a.r[6] < b.r[6];
      return a.r[7] < b.r[7];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector8<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector8<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector8<t> normalize(const vector8<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector8<t> & a, const vector8<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector8<t> vector8<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector8<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5 >> tmp6 >> tmp7;
      a = vector8<t>(tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector8<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4]<< " " << a.r[5]<< " " << a.r[6]<< " " << a.r[7];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 9.
    template<class t> class vector9{
    public:
      ///member
      t  r[9];
      ///zero constructor
      vector9(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; r[5]=0; r[6]=0; r[7]=0; r[8]=0; 
      }
      ///constructor from components
      vector9(t r0,t r1,t r2,t r3,t r4,t r5,t r6,t r7,t r8){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; r[5]=r5; r[6]=r6; r[7]=r7; r[8]=r8; 
      }
      ///constructor from a pointer
      vector9(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; r[5]=v[5]; r[6]=v[6]; r[7]=v[7]; r[8]=v[8]; 
      }
    
      ///default normalization;
      vector9 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];v[5]=r[5];v[6]=r[6];v[7]=r[7];v[8]=r[8];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
          case 5: return r[5];
          case 6: return r[6];
          case 7: return r[7];
          case 8: return r[8];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector9<t> operator+(const vector9<t> & a, const vector9<t> & b){
      return vector9<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4],a.r[5] + b.r[5],a.r[6] + b.r[6],a.r[7] + b.r[7],a.r[8] + b.r[8]);
    }
  
    ///vector subtraction
    template<class t>
      vector9<t> operator-(const vector9<t> & a, const vector9<t> & b){
      return vector9<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4],a.r[5] - b.r[5],a.r[6] - b.r[6],a.r[7] - b.r[7],a.r[8] - b.r[8]);
    }
  
    ///vector negate
    template<class t>
      vector9<t> operator-(const vector9<t> & a){
      return vector9<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4],-a.r[5],-a.r[6],-a.r[7],-a.r[8]);
    }
  
    ///scalar multiplication
    template<class t>
      vector9<t> operator*(const t &m, const vector9<t> & a){
      return vector9<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7],m*a.r[8]);
    }
    template<class t>
      vector9<t> operator*(const vector9<t> & a, const t &m){
      return vector9<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7],m*a.r[8]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector9<t> operator/(const vector9<t> & a, const t &m){
      return vector9<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m,a.r[5]/m,a.r[6]/m,a.r[7]/m,a.r[8]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector9<t> &a, const vector9<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4] + a.r[5]*b.r[5] + a.r[6]*b.r[6] + a.r[7]*b.r[7] + a.r[8]*b.r[8];
    }
  
  
    /// +=
    template<class t>
      vector9<t>& operator+=(vector9<t> &a, const vector9<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector9<t>& operator-=(vector9<t> &a, const vector9<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector9<t>& operator*=(vector9<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector9<t>& operator/=(vector9<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector9<t> & a, const vector9<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4] && a.r[5]==b.r[5] && a.r[6]==b.r[6] && a.r[7]==b.r[7] && a.r[8]==b.r[8];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector9<t> & a, const vector9<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      if(a.r[4] != b.r[4]) return a.r[4] < b.r[4];
      if(a.r[5] != b.r[5]) return a.r[5] < b.r[5];
      if(a.r[6] != b.r[6]) return a.r[6] < b.r[6];
      if(a.r[7] != b.r[7]) return a.r[7] < b.r[7];
      return a.r[8] < b.r[8];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector9<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector9<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector9<t> normalize(const vector9<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector9<t> & a, const vector9<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector9<t> vector9<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector9<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5 >> tmp6 >> tmp7 >> tmp8;
      a = vector9<t>(tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector9<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4]<< " " << a.r[5]<< " " << a.r[6]<< " " << a.r[7]<< " " << a.r[8];
    }
  
  }//component_by_array
  
  
  namespace component_by_array{    
    ///linear-algebraic vector of type t, dimension 10.
    template<class t> class vector10{
    public:
      ///member
      t  r[10];
      ///zero constructor
      vector10(){
        r[0]=0; r[1]=0; r[2]=0; r[3]=0; r[4]=0; r[5]=0; r[6]=0; r[7]=0; r[8]=0; r[9]=0; 
      }
      ///constructor from components
      vector10(t r0,t r1,t r2,t r3,t r4,t r5,t r6,t r7,t r8,t r9){
        r[0]=r0; r[1]=r1; r[2]=r2; r[3]=r3; r[4]=r4; r[5]=r5; r[6]=r6; r[7]=r7; r[8]=r8; r[9]=r9; 
      }
      ///constructor from a pointer
      vector10(t *v){
        r[0]=v[0]; r[1]=v[1]; r[2]=v[2]; r[3]=v[3]; r[4]=v[4]; r[5]=v[5]; r[6]=v[6]; r[7]=v[7]; r[8]=v[8]; r[9]=v[9]; 
      }
    
      ///default normalization;
      vector10 normalize(t a = 1);
      void writeToArray(t *v) const {
        v[0]=r[0];v[1]=r[1];v[2]=r[2];v[3]=r[3];v[4]=r[4];v[5]=r[5];v[6]=r[6];v[7]=r[7];v[8]=r[8];v[9]=r[9];;
      }
  
      ///component access operator.
      t& operator[](int index){
        switch(index){
          case 0: return r[0];
          case 1: return r[1];
          case 2: return r[2];
          case 3: return r[3];
          case 4: return r[4];
          case 5: return r[5];
          case 6: return r[6];
          case 7: return r[7];
          case 8: return r[8];
          case 9: return r[9];
        }
      }
    };
  
    ///vector addition
    template<class t>
      vector10<t> operator+(const vector10<t> & a, const vector10<t> & b){
      return vector10<t>(a.r[0] + b.r[0],a.r[1] + b.r[1],a.r[2] + b.r[2],a.r[3] + b.r[3],a.r[4] + b.r[4],a.r[5] + b.r[5],a.r[6] + b.r[6],a.r[7] + b.r[7],a.r[8] + b.r[8],a.r[9] + b.r[9]);
    }
  
    ///vector subtraction
    template<class t>
      vector10<t> operator-(const vector10<t> & a, const vector10<t> & b){
      return vector10<t>(a.r[0] - b.r[0],a.r[1] - b.r[1],a.r[2] - b.r[2],a.r[3] - b.r[3],a.r[4] - b.r[4],a.r[5] - b.r[5],a.r[6] - b.r[6],a.r[7] - b.r[7],a.r[8] - b.r[8],a.r[9] - b.r[9]);
    }
  
    ///vector negate
    template<class t>
      vector10<t> operator-(const vector10<t> & a){
      return vector10<t>(-a.r[0],-a.r[1],-a.r[2],-a.r[3],-a.r[4],-a.r[5],-a.r[6],-a.r[7],-a.r[8],-a.r[9]);
    }
  
    ///scalar multiplication
    template<class t>
      vector10<t> operator*(const t &m, const vector10<t> & a){
      return vector10<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7],m*a.r[8],m*a.r[9]);
    }
    template<class t>
      vector10<t> operator*(const vector10<t> & a, const t &m){
      return vector10<t>(m*a.r[0],m*a.r[1],m*a.r[2],m*a.r[3],m*a.r[4],m*a.r[5],m*a.r[6],m*a.r[7],m*a.r[8],m*a.r[9]);
    }
  
    ///scalar multiplication by reciprocal of \a m
    template<class t>
      vector10<t> operator/(const vector10<t> & a, const t &m){
      return vector10<t>(a.r[0]/m,a.r[1]/m,a.r[2]/m,a.r[3]/m,a.r[4]/m,a.r[5]/m,a.r[6]/m,a.r[7]/m,a.r[8]/m,a.r[9]/m);
    }
  
    ///inner product
    template<class t>
      t operator*(const vector10<t> &a, const vector10<t> & b){
      return a.r[0]*b.r[0] + a.r[1]*b.r[1] + a.r[2]*b.r[2] + a.r[3]*b.r[3] + a.r[4]*b.r[4] + a.r[5]*b.r[5] + a.r[6]*b.r[6] + a.r[7]*b.r[7] + a.r[8]*b.r[8] + a.r[9]*b.r[9];
    }
  
  
    /// +=
    template<class t>
      vector10<t>& operator+=(vector10<t> &a, const vector10<t> &b){
      return a=a+b;
    }
    /// -=
    template<class t>
      vector10<t>& operator-=(vector10<t> &a, const vector10<t> &b){
      return a=a-b;
    }
    /// *=
    template<class t>
      vector10<t>& operator*=(vector10<t> &a, const t &b){
      return a=b*a;
    }
    /// /=
    template<class t>
      vector10<t>& operator/=(vector10<t> &a, const t &b){
      return a=a/b;
    }
  
  
  
    ///equality between two vectors
    template<class t>
      bool operator==(const vector10<t> & a, const vector10<t> & b){
      return a.r[0]==b.r[0] && a.r[1]==b.r[1] && a.r[2]==b.r[2] && a.r[3]==b.r[3] && a.r[4]==b.r[4] && a.r[5]==b.r[5] && a.r[6]==b.r[6] && a.r[7]==b.r[7] && a.r[8]==b.r[8] && a.r[9]==b.r[9];
    }
  
    ///dictionary order of vectors
    template<class t>
      bool operator<(const vector10<t> & a, const vector10<t> & b){
      if(a.r[0] != b.r[0]) return a.r[0] < b.r[0];
      if(a.r[1] != b.r[1]) return a.r[1] < b.r[1];
      if(a.r[2] != b.r[2]) return a.r[2] < b.r[2];
      if(a.r[3] != b.r[3]) return a.r[3] < b.r[3];
      if(a.r[4] != b.r[4]) return a.r[4] < b.r[4];
      if(a.r[5] != b.r[5]) return a.r[5] < b.r[5];
      if(a.r[6] != b.r[6]) return a.r[6] < b.r[6];
      if(a.r[7] != b.r[7]) return a.r[7] < b.r[7];
      if(a.r[8] != b.r[8]) return a.r[8] < b.r[8];
      return a.r[9] < b.r[9];
    }
  
    ///square of absolute value
    template<class t>
      t norm(const vector10<t> & a){
      return a*a;
    }
  
    ///absolute value
    template<class t> 
      t abs(const vector10<t> & a){
      return std::sqrt(norm(a));
    }
  
    ///normalization (return new vector, of length \a a and parallel to \a *this)
    ///default value of \a a is 1
    template<class t> 
      vector10<t> normalize(const vector10<t> &a, t b=1){
      return (b/abs(a)) * (a);
    }
  
    ///projection factor of vector \a a in direction of \a b
    template<class t>
      t operator/(const vector10<t> & a, const vector10<t> & b){
      return (a*b) / (b*b);
    }
  
  
    /// normalization also provided as class method
    template<class t>
      vector10<t> vector10<t>::normalize(t a){
      return (a/abs(*this))*(*this);
    }
  
  
    ///read from istream
    template<class t> 
      std::istream& operator>>(std::istream& cin, vector10<t>& a){
      t tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9; 
      cin  >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5 >> tmp6 >> tmp7 >> tmp8 >> tmp9;
      a = vector10<t>(tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9);
      return cin;
    }
  
    ///write to ostream
    template<class t> 
      std::ostream& operator<<(std::ostream& cout, const vector10<t> &a){
      return cout << a.r[0]<< " " << a.r[1]<< " " << a.r[2]<< " " << a.r[3]<< " " << a.r[4]<< " " << a.r[5]<< " " << a.r[6]<< " " << a.r[7]<< " " << a.r[8]<< " " << a.r[9];
    }
  
  }//component_by_array
  
  
  }//vector
  }//math
  }//pfi
  #endif // #ifndef INCLUDE_GUARD_PFI_MATH_VECTOR_H_
  