/Users/jab/Documents/Teaching/02585/GEL2_and_demos/GEL/src/CGLA/ArithVec.h

#ifndef __CGLA_ARITHVEC_H__
#define __CGLA_ARITHVEC_H__

#include <iostream>
#include "CGLA.h"

#include <numeric>

namespace CGLA 
{

  template <class T, class V, unsigned int N> 
    class ArithVec
    {

    protected:

      T data[N];

    protected:

      //----------------------------------------------------------------------
      // Constructors
      //----------------------------------------------------------------------

      ArithVec() 
        {
        }

      explicit ArithVec(T _a)
        {
          std::fill_n(data, N, _a);
        }

      ArithVec(T _a, T _b)
        {
          assert(N==2);
          data[0] = _a;
          data[1] = _b;
        }

      ArithVec(T _a, T _b, T _c)
        {
          assert(N==3);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
        }

      ArithVec(T _a, T _b, T _c, T _d)
        {
          assert(N==4);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
          data[3] = _d;
        }


    public:

      typedef T ScalarType;
        
      typedef V VectorType;

      static unsigned int get_dim() {return N;}
        
      void set(T _a, T _b)
        {
          assert(N==2);
          data[0] = _a;
          data[1] = _b;
        }

      void set(T _a, T _b, T _c)
        {
          assert(N==3);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
        }

      void set(T _a, T _b, T _c, T _d)
        {
          assert(N==4);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
          data[3] = _d;
        }

  const T& operator [] ( unsigned int i ) const
        {
          assert(i<N);
          return data[i];
        }

  T& operator [] ( unsigned int i ) 
        {
          assert(i<N);
          return data[i];
        }

  const T& operator () ( unsigned int i ) const
        {
          assert(i<N);
          return data[i];
        }

  T& operator () ( unsigned int i ) 
        {
          assert(i<N);
          return data[i];
        }


      T* get() {return &data[0];}

      const T* get() const {return &data[0];}

      //----------------------------------------------------------------------
      // Comparison operators
      //----------------------------------------------------------------------

      bool operator==(const V& v) const 
        {
          return std::inner_product(data, &data[N], v.get(), true,
                                    std::logical_and<bool>(), std::equal_to<T>());
        }

#define for_all_i(expr) for(unsigned int i=0;i<N;i++) {expr;}

      bool operator==(T k) const 
        { 
          for_all_i(if (data[i] != k) return false)
            return true;
        }
#undef for_all_i  

      bool operator!=(const V& v) const 
        {
          return std::inner_product(data, &data[N], v.get(), false,
                                    std::logical_or<bool>(), std::not_equal_to<T>());
        }

      bool operator!=(T k) const 
        { 
          return !(*this==k);
        }


      //----------------------------------------------------------------------
      // Comparison functions ... of geometric significance 
      //----------------------------------------------------------------------

      bool  all_l  (const V& v) const
        {
          return std::inner_product(data, &data[N], v.get(), true, 
                                    std::logical_and<bool>(), std::less<T>());
        }

      bool  all_le (const V& v) const
        {
          return std::inner_product(data, &data[N], v.get(), true, 
                                    std::logical_and<bool>(), std::less_equal<T>());
        }

      bool  all_g  (const V& v) const
        {
          return std::inner_product(data, &data[N], v.get(), true, 
                                    std::logical_and<bool>(), std::greater<T>());
        }

      bool  all_ge (const V& v) const
        {
          return std::inner_product(data, &data[N], v.get(), true, 
                                    std::logical_and<bool>(), std::greater_equal<T>());
        }


      //----------------------------------------------------------------------
      // Assignment operators
      //----------------------------------------------------------------------

      const V& operator *=(T k) 
        { 
          std::transform(data, &data[N], data, std::bind2nd(std::multiplies<T>(), k));
          return static_cast<const V&>(*this);
        }

      const V& operator /=(T k)
        { 
          std::transform(data, &data[N], data, std::bind2nd(std::divides<T>(), k));
          return static_cast<const V&>(*this);
        }

      const V& operator +=(T k) 
        {
          std::transform(data, &data[N], data, std::bind2nd(std::plus<T>(), k));
          return  static_cast<const V&>(*this);
        }

      const V& operator -=(T k) 
        { 
          std::transform(data, &data[N], data, std::bind2nd(std::minus<T>(), k));
          return  static_cast<const V&>(*this);
        }

      const V& operator *=(const V& v) 
        { 
          std::transform(data, &data[N], v.get(), data, std::multiplies<T>());
          return  static_cast<const V&>(*this);
        }

      const V& operator /=(const V& v)
        {
          std::transform(data, &data[N], v.get(), data, std::divides<T>());
          return  static_cast<const V&>(*this);
        }

      const V& operator +=(const V& v) 
        {
          std::transform(data, &data[N], v.get(), data, std::plus<T>());
          return  static_cast<const V&>(*this);
        }
                
      const V& operator -=(const V& v) 
        { 
          std::transform(data, &data[N], v.get(), data, std::minus<T>());
          return  static_cast<const V&>(*this);
        }


      //----------------------------------------------------------------------
      // Unary operators on vectors
      //----------------------------------------------------------------------

      const V operator - () const
        {
          V v_new;
          std::transform(data, &data[N], v_new.get(), std::negate<T>());
          return v_new;
        }

      //----------------------------------------------------------------------
      // Binary operators on vectors
      //----------------------------------------------------------------------

      const V operator * (const V& v1) const
        {
          V v_new;
          std::transform(data, &data[N], v1.get(), v_new.get(), std::multiplies<T>());
          return v_new;
        }

      const V operator + (const V& v1) const
        {
          V v_new;
          std::transform(data, &data[N], v1.get(), v_new.get(), std::plus<T>());
          return v_new;
        }

      const V operator - (const V& v1) const
        {
          V v_new;
          std::transform(data, &data[N], v1.get(), v_new.get(), std::minus<T>());
          return v_new;
        }

      const V operator / (const V& v1) const
        {
          V v_new;
          std::transform(data, &data[N], v1.get(), v_new.get(), std::divides<T>());
          return v_new;
        }

      //----------------------------------------------------------------------
      // Binary operators on vector and scalar
      //----------------------------------------------------------------------

      const V operator * (T k) const
        {
          V v_new;
          std::transform(data, &data[N], v_new.get(), std::bind2nd(std::multiplies<T>(), k));
          return v_new;
        }
  

      const V operator / (T k) const
        {
          V v_new;
          std::transform(data, &data[N], v_new.get(), std::bind2nd(std::divides<T>(), k));
          return v_new;      
        }


      const T min_coord() const 
        {
          return *std::min_element(data, &data[N]);
        }

      const T max_coord() const
        {
          return *std::max_element(data, &data[N]);
        }

    };

  template <class T, class V, unsigned int N> 
    inline std::ostream& operator<<(std::ostream&os, const ArithVec<T,V,N>& v)
    {
      os << "[ ";
      for(unsigned int i=0;i<N;i++) os << v[i] << " ";
      os << "]";
      return os;
    }

  template <class T,class V, unsigned int N>
    inline std::istream& operator>>(std::istream&is, ArithVec<T,V,N>& v)
    {
      for(unsigned int i=0;i<N;i++) is>>v[i];
      return is;
    }


  template <class T,class V, unsigned int N>
    inline T dot(const ArithVec<T,V,N>& v0, const ArithVec<T,V,N>& v1)
    {
      return std::inner_product(v0.get(), v0.get() + N, v1.get(), T(0));
    }

  template <class T,class V, unsigned int N>
    inline T sqr_length(const ArithVec<T,V,N>& v)
    {
      return dot(v,v);
    }

  template<class T, class V, unsigned int N>
    inline const V operator * (double k, const ArithVec<T,V,N>& v) 
    {
      return v * k;
    }

  template<class T, class V, unsigned int N>
    inline const V operator * (float k, const ArithVec<T,V,N>& v) 
    {
      return v * k;
    }

  template<class T, class V, unsigned int N>
    inline const V operator * (int k, const ArithVec<T,V,N>& v) 
    {
      return v * k;
    }

  template <class T,class V, unsigned int N>
    inline V v_min(const ArithVec<T,V,N>& v0, const ArithVec<T,V,N>& v1)
    {
      V v;
      std::transform(v0.get(), v0.get() + N, v1.get(), v.get(), std::ptr_fun(s_min<T>));
      return v;
    }

  template <class T,class V, unsigned int N>
    inline V v_max(const ArithVec<T,V,N>& v0, const ArithVec<T,V,N>& v1)
    {
      V v;
      std::transform(v0.get(), v0.get() + N, v1.get(), v.get(), std::ptr_fun(s_max<T>));
      return v;
    }


}

#endif