#ifndef VEC_H
#define VEC_H

#include <string>
#include <sstream>
using namespace std;

#include <typedefs.h>


/*!
\brief Primitive vector.

Support arithmetic operations on a data type
 with the [] operator.

This class does not explicity manage memory.
 No array bounds error checking or anything,
 nothing to interfere with the cache. Let
 the emphasis be on the client data structures to provide 
 such functionality if they choose to.

\par Example
\verbatim
  unsigned int const n=3;
  double a1[n] = { 2.0, 7.0, 0.0 };
  double a2[n] = { 1.0, -3.0, 2.0 };
  double a3[n] = { 0.0, 0.0, 0.0 };

  vec<double*,double> v1(&a1[0],n);
  vec<double*,double> v2(&a2[0],n);
  vec<double*,double> v3(&a3[0],n);

  // Calculate equation
  // a3 = 2(a1-a2)
  v3 = a1;
  v3 -= a2;
  v3 *= 2.0;
\endverbatim  
*/
template < typename A, typename T >
class vec
{
public:

  /** The vector. */
  A ai;

  /** The number of dimensions. */
  uintc N;

  /** Pass in the vector and it number of dimensions. */
  vec(A ai_, uintc N_)
    : ai(ai_), N(N_) {}

  /** Scalar assigned to each component. */
  void operator = ( T const x );

  /** Copy vector to ai. */
  template< typename Z >
  void operator = ( Z const & zi );

  /** Copy vector to ai. */
  void operator = (vec<A,T> const & vb);

  /** Matching component multiplication */
  template< typename Z >
  void operator *= ( Z const & zi );

  /** Matching component addition. */
  template< typename Z >
  void operator += ( Z const & zi );

  /** Matching component subraction. */
  template< typename Z >
  void operator -= ( Z const & zi );

  /** Scalar addition. */
  void operator += ( T const x );

  /** Scalar multiplication. */
  void operator *= ( T const x );

  /** Element access. */
  T const operator [] (uintc i) const
    { return ai[i]; }

  /** Element access. */
  T & operator [] (uintc i)
    { return ai[i]; }

  
  /** Set the value in the ith component, all others zero. */
  void identity(uintc k, T const x=1);

  /** Difference X=X-rightshift(X) */
  void difference();
 
  /** X=X+rightshif(X) */
  void sum();

  /** Generate the binomial coefficients for N degree polynomial. */
  void binomial(uintc N);

  /** Write the vector out as a string separated with spaces. */
  operator string () const;

};


//---------------------------------------------------------
//  Implementation

template < typename A, typename T >
vec<A,T>::operator string() const
{
  if (N==0)
    return "";

  string s;
  { stringstream ss;  ss << ai[0]; s+=ss.str(); }
  
  for (uint i=1; i<N; ++i)
  {
    s += " ";
    { stringstream ss;  ss << ai[i]; s+=ss.str(); }
  }
  
  return s;
}

template < typename A, typename T >
void vec<A,T>::binomial(uintc N)
{
  identity(0);
  for (uint i=0; i<N; ++i)
    sum();
}

template < typename A, typename T >
void vec<A,T>::identity(uintc k, T const x)
{
  for (uint i=0; i<N; ++i)
    ai[i] = 0;
  ai[k]=x;
}

template < typename A, typename T >
void vec<A,T>::difference()
{
  for (uint i=N-1; i>0; --i)
    ai[i] = ai[i] - ai[i-1];
}

template < typename A, typename T >
void vec<A,T>::sum()
{
  for (uint i=N-1; i>0; --i)
    ai[i] = ai[i] + ai[i-1];
}

template < typename A, typename T >
template< typename Z >
void vec<A,T>::operator = ( Z const & zi )
{
  for (uint i=0; i<N; ++i)
    ai[i] = zi[i];
}

template < typename A, typename T >
template< typename Z >
void vec<A,T>::operator *= ( Z const & zi )
{
  for (uint i=0; i<N; ++i)
    ai[i] *= zi[i];
}

template < typename A, typename T >
template< typename Z >
void vec<A,T>::operator += ( Z const & zi )
{
  for (uint i=0; i<N; ++i)
    ai[i] += zi[i];
}

template < typename A, typename T >
template< typename Z >
void vec<A,T>::operator -= ( Z const & zi )
{
  for (uint i=0; i<N; ++i)
    ai[i] -= zi[i];
}

template < typename A, typename T >
void vec<A,T>::operator = ( T x )
{
  for (uint i=0; i<N; ++i)
    ai[i] = x; 
}

template < typename A, typename T >
void vec<A,T>::operator *= ( T x )
{
  for (uint i=0; i<N; ++i)
    ai[i] *= x; 
}

template < typename A, typename T >
void vec<A,T>::operator += ( T x )
{
  for (uint i=0; i<N; ++i)
    ai[i] += x; 
}

template < typename A, typename T >
void vec<A,T>::operator = ( vec<A,T> const & bi )
{
  for (uint i=0; i<N; ++i)
    ai[i] = bi[i];
}


#endif