/**********************************************************

This software is part of J.-S. Caux's ABACUS++ library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  ABACUS_Vect.h

Purpose:  Declares vector class.

***********************************************************/

#ifndef ABACUS_VECT_H
#define ABACUS_VECT_H

namespace ABACUS {

  template <class T>
    class Vect {
  private:
    int dim;
    T* V;
  public:
    Vect();
    explicit Vect (int N);
    Vect (const T& a, int N);  // initialize all N elements are a
    Vect (const T* a, int N);  // initialize to array
    Vect (const Vect& rhs);   // Copy constructor
    Vect& operator= (const Vect& rhs);  // assignment
    Vect& operator= (const T& a);      // assign a to all elements
    inline T& operator[] (const int i);
    inline const T& operator[] (const int i) const;
    Vect& operator+= (const Vect& rhs);
    Vect& operator-= (const Vect& rhs);
    bool operator== (const Vect& rhs);  // checks equality of size and of all elements
    bool operator!= (const Vect& rhs);  // checks inequality
    bool Append (const T& rhs);  // appends rhs to the vector
    bool Append (const Vect& rhs);  // appends rhs to the vector
    bool Increase_Size (int nr_to_add); // resizes the array to accommodate nr_to_add more entries
    bool Increase_Size (int nr_to_add, T setval); // resizes the array to accommodate nr_to_add more entries
    inline int size() const;
    inline double norm() const;  // returns norm of vector
    inline T max() const;        // returns maximal value
    inline T min() const;        // returns maximal value
    inline T sum() const;        // returns sum of all elements
    inline bool includes(T check) const;  // whether check == one of the elements or not
    void QuickSort (int l, int r);
    void QuickSort (Vect<int>& index, int l, int r);
    void QuickSort ();
    void QuickSort (Vect<int>& index);
    ~Vect();
  };

  template <class T>
    Vect<T>::Vect() : dim(0), V(0) {}

  template <class T>
    Vect<T>::Vect (int N) : dim(N), V(new T[N]) {}

  template <class T>
    Vect<T>::Vect (const T& a, int N) : dim(N), V(new T[N])
    {
      for (int i = 0; i < N; ++i) V[i] = a;
    }

  template <class T>
    Vect<T>::Vect (const T* a, int N) : dim(N), V(new T[N])
    {
      for (int i = 0; i < N; ++i) V[i] = *a++;
    }

  template <class T>
    Vect<T>::Vect (const Vect<T>& rhs) : dim(rhs.dim), V(new T[dim])
    {
      for (int i = 0; i < dim; ++i) V[i] = rhs[i];
    }

  template <class T>
    Vect<T>& Vect<T>::operator= (const Vect<T>& rhs)
    {
      if (this != &rhs) {
	if (dim != rhs.dim) {
	  if (V != 0) delete[] V;
	  dim = rhs.dim;
	  V = new T[dim];
	}
	for (int i = 0; i < dim; ++i) V[i] = rhs[i];
      }
      return *this;
    }

  template <class T>
    Vect<T>& Vect<T>::operator= (const T& a)
    {
      for (int i = 0; i < dim; ++i) V[i] = a;
      return *this;
    }

  template <class T>
    inline T& Vect<T>::operator[] (const int i)
    {
      return V[i];
    }

  template <class T>
    inline const T& Vect<T>::operator[] (const int i) const
    {
      return V[i];
    }

  template <class T>
    Vect<T>& Vect<T>::operator+= (const Vect<T>& rhs)
    {
      for (int i = 0; i < dim; ++i) V[i] += rhs[i];
      return *this;
    }

  template <class T>
    Vect<T>& Vect<T>::operator-= (const Vect<T>& rhs)
    {
      for (int i = 0; i < dim; ++i) V[i] -= rhs[i];
      return *this;
    }

  template <class T>
    bool Vect<T>::operator== (const Vect<T>& rhs)
    {
      bool answer = ((*this).size() == rhs.size());
      if (answer) {
	for (int i = 0; i < dim; ++i) answer = (answer && (V[i] == rhs[i]));
      }
      return answer;
    }

  template <class T>
    bool Vect<T>::operator!= (const Vect<T>& rhs)
    {
      return(!((*this) == rhs));
    }

  template <class T>
    bool Vect<T>::Append (const Vect<T>& rhs)  // appends rhs to the vector
    {
      T* newvect = new T[dim + rhs.size()];
      for (int i = 0; i < dim; ++i) newvect[i] = V[i];
      for (int i = 0; i < rhs.size(); ++i) newvect[i+ dim] = rhs[i];

      dim += rhs.size();
      delete[] V;
      V = new T[dim];
      for (int i = 0; i < dim; ++i) V[i] = newvect[i];

      delete[] newvect;

      return(true);
    }

  template <class T>
    bool Vect<T>::Append (const T& rhs)  // appends rhs to the vector
    {
      T* newvect = new T[dim + 1];
      for (int i = 0; i < dim; ++i) newvect[i] = V[i];
      newvect[dim] = rhs;

      dim += 1;
      delete[] V;
      V = new T[dim];
      for (int i = 0; i < dim; ++i) V[i] = newvect[i];

      delete[] newvect;

      return(true);
    }

  template <class T>
    bool Vect<T>::Increase_Size (int nr_to_add) // resizes the array to accommodate nr_to_add more entries
    {
      int resized_dim = dim + nr_to_add;

      T* resized_vect = new T[resized_dim];
      for (int i = 0; i < dim; ++i) resized_vect[i] = V[i];
      for (int i = dim; i < resized_dim; ++i) resized_vect[i] = T(0);

      dim = resized_dim;
      delete[] V;
      V = new T[dim];
      for (int i = 0; i < dim; ++i) V[i] = resized_vect[i];

      delete[] resized_vect;

      return(true);
    }

  template <class T>
    bool Vect<T>::Increase_Size (int nr_to_add, T setval) // resizes the array to accommodate nr_to_add more entries
    {
      int resized_dim = dim + nr_to_add;

      T* resized_vect = new T[resized_dim];
      for (int i = 0; i < dim; ++i) resized_vect[i] = V[i];
      for (int i = dim; i < resized_dim; ++i) resized_vect[i] = setval;

      dim = resized_dim;
      delete[] V;
      V = new T[dim];
      for (int i = 0; i < dim; ++i) V[i] = resized_vect[i];

      delete[] resized_vect;

      return(true);
    }


  template <class T>
    inline int Vect<T>::size() const
    {
      return dim;
    }

  template <class T>
    inline double Vect<T>::norm () const
    {
      double normsq = 0.0;
      for (int i = 0; i < dim; ++i) normsq += abs(V[i]) * abs(V[i]);
      return sqrt(normsq);
    }

  template <>
    inline double Vect<double>::norm () const
    {
      double normsq = 0.0;
      for (int i = 0; i < dim; ++i) normsq += V[i] * V[i];
      return(sqrt(normsq));
    }

  template <>
    inline double Vect<std::complex<double> >::norm () const
    {
      double normsq = 0.0;
      for (int i = 0; i < dim; ++i) normsq += std::norm(V[i]);
      return(sqrt(normsq));
    }

  template <class T>
    inline T Vect<T>::max() const
    {
      T maxval = V[0];
      for (int i = 0; i < dim; ++i) if (V[i] > maxval) maxval = V[i];
      return maxval;
    }

  template <class T>
    inline T Vect<T>::min() const
    {
      T minval = V[0];
      for (int i = 0; i < dim; ++i) if (V[i] < minval) minval = V[i];
      return minval;
    }

  template <class T>
    inline T Vect<T>::sum() const
    {
      T total = T(0);
      for (int i = 0; i < dim; ++i) total += V[i];
      return total;
    }

  template <class T>
    inline bool Vect<T>::includes (T check) const
    {
      int index = 0;
      while (index < dim && V[index] != check) index++;

      return(index < dim);
    }

  template <class T>
    void Vect<T>::QuickSort (int l, int r)
    {
      int i = l, j = r;
      T pivot = V[l + (r-l)/2];

      while (i <= j) {
	while (V[i] < pivot) i++;
	while (V[j] > pivot) j--;
	if (i <= j) {
	  std::swap(V[i],V[j]);
	  i++;
	  j--;
	}
      };

      if (l < j) (*this).QuickSort(l, j);
      if (i < r) (*this).QuickSort(i, r);
    }

  template <class T>
    void Vect<T>::QuickSort ()
    {
      if ((*this).size() > 1) (*this).QuickSort (0, (*this).size() - 1);
    }

  template <class T>
    void Vect<T>::QuickSort (Vect<int>& index, int l, int r)
    {
      int i = l, j = r;
      T pivot = V[l + (r-l)/2];

      while (i <= j) {
	while (V[i] < pivot) i++;
	while (V[j] > pivot) j--;
	if (i <= j) {
	  std::swap(V[i],V[j]);
	  std::swap(index[i],index[j]);
	  i++;
	  j--;
	}
      };

      if (l < j) (*this).QuickSort(index, l, j);
      if (i < r) (*this).QuickSort(index, i, r);
    }

  template <class T>
    void Vect<T>::QuickSort (Vect<int>& index)
    {
      if (index.size() != (*this).size()) ABACUSerror("Wrong dim for index in Vect QuickSort.");
      (*this).QuickSort (index, 0, (*this).size() - 1);
    }

  template <class T>
    inline std::ostream& operator<< (std::ostream& s, const Vect<T>& vector)
    {
      for (int i = 0; i < vector.size() - 1; ++i) s << vector[i] << " ";
      if (vector.size() >= 1) s << vector[vector.size() - 1];

      return (s);
    }

  template <class T>
    Vect<T>::~Vect<T>()
    {
      if (V != 0) delete[] V;
    }


  // TYPEDEFS
  typedef ABACUS::Vect<int> Vect_INT;
  typedef ABACUS::Vect<double> Vect_DP;
  typedef ABACUS::Vect<std::complex<double> > Vect_CX;


} // namespace ABACUS

#endif