/* Copyright (c) 1997-2006
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
*/

#ifndef _POLYMAKE_GENERIC_VECTOR_H
#define _POLYMAKE_GENERIC_VECTOR_H "$Project: polymake $$Id: GenericVector.h 7402 2006-07-27 23:44:14Z gawrilow $"

#include <algorithm>

#ifndef _POLYMAKE_PLAUSIBLE_CHECKS_H
#  include <plausible_checks.h>
#endif
#ifndef _POLYMAKE_GENERIC_IO_H
#  include <GenericIO.h>
#endif
#ifndef _POLYMAKE_TRANSFORMED_CONTAINER_H
#  include <TransformedContainer.h>
#endif
#ifndef _POLYMAKE_CONVERTERS_H
#  include <converters.h>
#endif
#ifndef _POLYMAKE_INDEXED_SUBSET_H
#  include <IndexedSubset.h>
#endif
#ifndef _POLYMAKE_INDEXED_SUBSET_H
#  include <SelectedSubset.h>
#endif
#ifndef _POLYMAKE_CONTAINER_CHAIN_H
#  include <ContainerChain.h>
#endif
#ifndef _POLYMAKE_SPARSE_H
#  include <sparse.h>
#endif

namespace pm {

/* ----------------
 *  Generic vector
 * ---------------- */

template <typename E> class Vector;
template <typename E, typename Predicate=conv<E,bool> > class SparseVector;
template <typename VectorRef, typename Operation> class LazyVector1;
template <typename VectorRef1, typename VectorRef2, typename Operation> class LazyVector2;

template <typename _Vector, typename E=typename _Vector::element_type>
class GenericVector : public Generic<_Vector>, public operators::base {
protected:
   GenericVector() { }
   GenericVector(const GenericVector&) { }
public:
   typedef E element_type;
   typedef typename if_else<check_container_feature<_Vector,sparse>::value, SparseVector<E>, Vector<E> >::type persistent_type;
   typedef GenericVector generic_type;
   typedef typename if_else<check_container_feature<_Vector,pure_sparse>::value, pure_sparse,
	   typename if_else<check_container_feature<_Vector,sparse>::value, sparse, dense>::type>::type
      sparse_discr;
protected:
   bool is_shared() const { return false; }

   template <typename Container>
   void _assign(const Container& c, dense)
   {
      std::copy(ensure(c, (dense*)0).begin(), entire(*this));
   }

   template <typename Container>
   void _assign(const Container& c, pure_sparse)
   {
      assign_sparse(this->top(), ensure(c, (pure_sparse*)0).begin());
   }

   template <typename _Vector2>
   bool trivial_assignment(const GenericVector<_Vector2,E>&) const { return false; }

   bool trivial_assignment(const GenericVector& v) const { return this==&v; }

   template <typename _Vector2>
   void swap(GenericVector<_Vector2,E>& v, dense)
   {
      std::swap_ranges(entire(*this), v.top().begin());
   }

   template <typename _Vector2>
   void swap(GenericVector<_Vector2,E>& v, pure_sparse)
   {
      swap_sparse(this->top(), ensure(v.top(), (pure_sparse*)0));
   }

   template <typename _Vector2>
   void assign(const GenericVector<_Vector2>& v)
   {
      DIMENSION_CHECK(v.dim() != dim(), "operator= - vector dimension mismatch");
      _assign(v.top(), sparse_discr());
   }

   template <typename E2, int _size>
   void assign(const E2 (&a)[_size])
   {
      DIMENSION_CHECK(_size != dim(), "operator= - vector dimension mismatch");
      _assign(array2container(a), sparse_discr());
   }

   template <typename Operation>
   void _assign_op(const Operation& op_arg)
   {
      perform_assign(entire(*this), op_arg);
   }

   template <typename Container2, typename Operation>
   void _assign_op(const Container2& c2, const Operation& op_arg, dense, dense)
   {
      perform_assign(entire(*this), c2.begin(), op_arg);
   }

   template <typename Container2, typename Operation>
   void _assign_op(const Container2& c2, const Operation& op_arg, dense, sparse)
   {
      typedef typename Entire<Container2>::const_iterator iterator2;
      typedef binary_op_builder<Operation, typename _Vector::const_iterator, iterator2> opb;
      const typename opb::operation& op=opb::create(op_arg);
      int i_prev=0;
      typename _Vector::iterator dst=this->top().begin();
      for (iterator2 src2=entire(c2); !src2.at_end(); ++src2) {
	 int i=src2.index();
	 std::advance(dst, i-i_prev);
	 op.assign(*dst, *src2);
	 i_prev=i;
      }
   }

   template <typename Container2, typename Operation, typename discr2>
   void _assign_op(const Container2& c2, const Operation& op_arg, pure_sparse, discr2)
   {
      typedef typename ensure_features<Container2, pure_sparse>::const_iterator iterator2;
      typedef binary_op_builder<Operation, typename _Vector::const_iterator, iterator2> opb;
      const typename opb::operation& op=opb::create(op_arg);
      _Vector& me=this->top();
      typename _Vector::iterator dst=me.begin();
      iterator2 src2=ensure(c2, (pure_sparse*)0).begin();
      int state=(dst.at_end() ? 0 : zipper_first) + (src2.at_end() ? 0 : zipper_second);
      while (state >= zipper_both) {
	 const int idiff=dst.index()-src2.index();
	 if (idiff<0) {
	    ++dst;
	    if (dst.at_end()) state -= zipper_first;
	 } else if (idiff>0) {
	    me.insert(dst, src2.index(), op(operations::partial_right(), dst, *src2));
	    ++src2;
	    if (src2.at_end()) state -= zipper_second;
	 } else {
	    if (op.assign(*dst,*src2)) ++dst;
	    else me.erase(dst++);
	    if (dst.at_end()) state -= zipper_first;
	    ++src2;
	    if (src2.at_end()) state -= zipper_second;
	 }
      }
      if (state & zipper_second) {
	 do {
	    me.insert(dst, src2.index(), op(operations::partial_right(), dst, *src2));  ++src2;
	 } while (!src2.at_end());
      }
   }

   template <typename E2>
   void _fill(const E2& x, dense)
   {
      std::fill(entire(*this), x);
   }

   template <typename E2>
   void _fill(const E2& x, pure_sparse)
   {
      _Vector& me=this->top();
      typename _Vector::iterator e=me.begin();
      if (x) {
	 for (int i=0, n=dim(); i<n; ++i) {
	    if (e.at_end() || i<e.index()) {
	       me.insert(e,i,x);
	    } else {
	       *e=x;  ++e;
	    }
	 }
      } else {
	 while (!e.at_end())
	    me.erase(e++);
      }
   }

   // undefined for dense
   void remove0s(dense);

   void remove0s(pure_sparse)
   {
      _Vector& me=this->top();
      typename _Vector::iterator e=me.begin();
      while (!e.at_end())
	 if (*e) ++e;
	 else me.erase(e++);
   }

public:
   int dim() const { return get_dim(this->top()); }

   template <typename _Vector2>
   void swap(GenericVector<_Vector2,E>& v)
   {
      DIMENSION_CHECK(dim() != v.dim(), "swap - vector dimension mismatch");
      typedef typename if_else<check_container_feature<_Vector,sparse>::value ||
	                       check_container_feature<_Vector2,sparse>::value,
	                       pure_sparse, dense>::type
	 sparse_discr2;
      if (!trivial_assignment(v))
	 swap(v, sparse_discr2());
   }

   template <typename Operation>
   void assign_op(const Operation& op)
   {
      if (this->top().is_shared())
	 this->top() = LazyVector1<const _Vector&, Operation>(this->top(),op);
      else
	 _assign_op(op);
   }

   template <typename Container2, typename Operation>
   void assign_op(const Container2& c, const Operation& op)
   {
      typedef typename if_else<check_container_feature<typename Concrete<Container2>::type,unlimited>::value,
	                       dense, sparse_discr>::type
	 sparse_discr1;
      typedef typename if_else<check_container_feature<typename Concrete<Container2>::type,sparse>::value,
	                       sparse, dense>::type
	 sparse_discr2;
      if (this->top().is_shared())
	 this->top() = LazyVector2<const _Vector&, const typename Concrete<Container2>::type&, Operation>
	                          (this->top(), Concrete<Container2>::get(c), op);
      else
	 _assign_op(Concrete<Container2>::get(c), op, sparse_discr1(), sparse_discr2());
   }

   template <typename E2>
   void fill(const E2& x)
   {
      this->top()._fill(x, sparse_discr());
   }

   _Vector& operator= (const GenericVector& v)
   {
      if (!trivial_assignment(v))
	 this->top().assign(v);
      return this->top();
   }

   template <typename Vector2>
   _Vector& operator= (const GenericVector<Vector2>& v)
   {
      this->top().assign(v);
      return this->top();
   }

   template <typename E2, int _size>
   _Vector& operator= (const E2 (&a)[_size])
   {
      this->top().assign(a);
      return this->top();
   }

   _Vector& negate()
   {
      assign_op(BuildUnary<operations::neg>());
      return this->top();
   }

   friend
   _Vector& inv_sign(GenericVector& me)
   {
      return me.negate();
   }

   friend
   _Vector& inv_sign(const GenericVector& me)
   {
      return non_const(me.top()).negate();
   }

   template <typename Vector2>
   _Vector& operator+= (const GenericVector<Vector2>& v)
   {
      DIMENSION_CHECK((dim() != v.dim()), "operator+= - vector dimension mismatch");
      assign_op(v, BuildBinary<operations::add>());
      return this->top();
   }

   template <typename Vector2>
   _Vector& operator-= (const GenericVector<Vector2>& v)
   {
      DIMENSION_CHECK((dim() != v.dim()), "operator-= - vector dimension mismatch");
      assign_op(v, BuildBinary<operations::sub>());
      return this->top();
   }

   template <typename Right>
   _Vector& operator*= (const Right& r)
   {
      if (r)
	 assign_op(constant(r), BuildBinary<operations::mul>());
      else
	 fill(r);
      return this->top();
   }

   template <typename Right>
   _Vector& operator/= (const Right& r)
   {
      assign_op(constant(r), BuildBinary<operations::div>());
      return this->top();
   }

   /// Divides by the first element.
   _Vector& dehomogenize()
   {
      const E first=this->top().front();
      return this->top()/=first;
   }

   /// Compare with an implicit zero vector
   bool operator! () const
   {
      return entire(attach_selector(this->top(), conv<E,bool>())).at_end();
   }

   /// Look for a non-zero element
   operator bool() const
   {
      return ! operator!();
   }

   /// Remove all zero elements which might have been overseen in some previous operation.
   void remove0s()
   {
      remove0s((sparse_discr)0);
   }

   template <typename IndexSet>
   IndexedSlice<_Vector&, const typename Concrete<IndexSet>::type&>
   slice(const IndexSet& indices)
   {
      INDEX_CHECK(!complement_helper<const typename Concrete<IndexSet>::type&>::dimension_check(Concrete<IndexSet>::get(indices), dim()),
		  "slice - indices out of range");
      return IndexedSlice<_Vector&, const typename Concrete<IndexSet>::type&>(this->top(), Concrete<IndexSet>::get(indices));
   }

   template <typename IndexSet>
   const IndexedSlice<const _Vector&, const typename Concrete<IndexSet>::type&>
   slice(const IndexSet& indices) const
   {
      INDEX_CHECK(!complement_helper<const typename Concrete<IndexSet>::type&>::dimension_check(Concrete<IndexSet>::get(indices), dim()),
		  "slice - indices out of range");
      return IndexedSlice<const _Vector&, const typename Concrete<IndexSet>::type&>(this->top(), Concrete<IndexSet>::get(indices));
   }

   IndexedSlice<_Vector&, sequence>
   slice(int sstart, int ssize=0)
   {
      INDEX_CHECK((sstart<0 || sstart+ssize>dim()), "slice - indices out of range");
      if (!ssize) ssize=dim()-sstart;
      return IndexedSlice<_Vector&, sequence>(this->top(), sequence(sstart,ssize));
   }

   const IndexedSlice<const _Vector&, sequence>
   slice(int sstart, int ssize=0) const
   {
      INDEX_CHECK((sstart<0 || sstart+ssize>dim()), "slice - indices out of range");
      if (!ssize) ssize=dim()-sstart;
      return IndexedSlice<const _Vector&, sequence>(this->top(), sequence(sstart,ssize));
   }

   template <typename Result>
   struct rebind_generic {
      typedef GenericVector<Result, E> type;
   };
};

struct is_vector;

template <typename Vector, typename E>
struct spec_object_traits< GenericVector<Vector,E> >
   : spec_or_model_traits<Vector,is_container> {
   typedef is_vector generic_tag;
   static const bool allow_sparse=true;
};

/* ---------------------------------------
 *  FixedVector
 * a built-in array decorated as a vector
 * --------------------------------------- */

template <typename E, int _size>
class FixedVector
   : public fixed_array<E,_size>,
     public GenericVector<FixedVector<E,_size>, E> {
public:
   FixedVector(int=0)
   {
      if (is_pod<E>::value) clear();
   }

   FixedVector(int, const E& init)
   {
      std::fill(entire(*this), init);
   }

   FixedVector(const E (&init)[_size])
   {
      std::copy(entire(*this), init+0);
   }

   template <typename E2>
   explicit FixedVector(const E2 (&init)[_size])
   {
      std::copy(entire(*this), attach_converter<E>(array2container(init)).begin());
   }

   template <typename Iterator>
   FixedVector(int, Iterator src)
   {
      std::copy(src, entire(*this));
   }

   template <typename Vector>
   FixedVector(const GenericVector<Vector,E>& v)
   {
      DIMENSION_CHECK(v.dim() != this->size(), "FixedVector - dimension mismatch in constructor");
      this->_assign(v.top(), dense());
   }

   template <typename Vector, typename E2>
   explicit FixedVector(const GenericVector<Vector,E2>& v)
   {
      DIMENSION_CHECK(v.dim() != this->size(), "FixedVector - dimension mismatch in constructor");
      this->_assign(attach_converter<E>(v.top()), dense());
   }

   FixedVector& operator= (const FixedVector& v)
   {
      return FixedVector::generic_type::operator=(v);
   }
   using FixedVector::generic_type::operator=;

   void clear()
   {
      operations::clear<E> zero;
      std::fill(entire(*this), zero());
   }
protected:
   friend class GenericVector<FixedVector>;
};

template <typename E, int _size>
struct spec_object_traits< FixedVector<E,_size> >
   : spec_object_traits<is_container> {
   static const int is_resizeable=0;
};

/* --------------------------------------------
 *  LazyVector1
 * lazy evaluation of an unary vector operator
 * -------------------------------------------- */

template <typename VectorRef, typename Operation>
class LazyVector1
   : public TransformedContainer<VectorRef, Operation>,
     public GenericVector< LazyVector1<VectorRef,Operation>,
                           typename object_traits<typename TransformedContainer<VectorRef,Operation>::value_type>::persistent_type > {
   typedef TransformedContainer<VectorRef, Operation> _super;
public:
   LazyVector1(typename _super::arg_type src_arg, const Operation& op_arg=Operation())
      : _super(src_arg, op_arg) { }

   int dim() const { return get_dim(this->get_container()); }
};

template <typename VectorRef, typename Operation>
struct spec_object_traits< LazyVector1<VectorRef,Operation> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true, is_temporary=true, is_always_const=true;
};

template <typename VectorRef, typename Operation, typename Feature>
struct check_container_feature<LazyVector1<VectorRef,Operation>, Feature>
   : check_container_feature<TransformedContainer<VectorRef,Operation>, Feature> { };

/* --------------------------------------------
 *  LazyVector2
 * lazy evaluation of a binary vector operator
 * -------------------------------------------- */

template <typename VectorRef1, typename VectorRef2, typename Operation>
class LazyVector2
   : public TransformedContainerPair<VectorRef1, VectorRef2, Operation>,
     public GenericVector< LazyVector2<VectorRef1,VectorRef2,Operation>,
                           typename object_traits<typename TransformedContainerPair<VectorRef1,VectorRef2,Operation>::value_type>::persistent_type > {
   typedef TransformedContainerPair<VectorRef1, VectorRef2, Operation> _super;
protected:
   int dim(True) const { return get_dim(this->get_container2()); }
   int dim(False) const { return get_dim(this->get_container1()); }
public:
   LazyVector2(typename _super::first_arg_type src1_arg, typename _super::second_arg_type src2_arg, const Operation& op_arg=Operation())
      : _super(src1_arg, src2_arg, op_arg) { }

   int dim() const
   {
      return dim(bool2type< check_container_ref_feature<VectorRef1,unlimited>::value >());
   }
};

template <typename VectorRef1, typename VectorRef2, typename Operation>
struct spec_object_traits< LazyVector2<VectorRef1, VectorRef2, Operation> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true, is_temporary=true, is_always_const=true;
};

template <typename VectorRef1, typename VectorRef2, typename Operation, typename Feature>
struct check_container_feature<LazyVector2<VectorRef1, VectorRef2, Operation>, Feature>
   : check_container_feature<TransformedContainerPair<VectorRef1, VectorRef2, Operation>, Feature> { };

/// explicit conversion of vector elements to another type
template <typename TargetType, typename Vector> inline
const LazyVector1<const Vector&, conv<typename Vector::element_type, TargetType> >
convert_to(const GenericVector<Vector>& V)
{
   return V.top();
}

/* --------------------------------------------
 *  VectorTensorProduct
 * lazy evaluation of a vector tensor product
 * -------------------------------------------- */
template <typename VectorRef1, typename VectorRef2, typename Operation=BuildBinary<operations::mul> >
class VectorTensorProduct
   : public ContainerProduct<VectorRef1, VectorRef2, Operation>,
     public GenericVector< VectorTensorProduct<VectorRef1,VectorRef2,Operation>,
			   typename object_traits<typename ContainerProduct<VectorRef1, VectorRef2, Operation>::value_type>::persistent_type  > {
   typedef ContainerProduct<VectorRef1, VectorRef2, Operation> _super;
public:
   VectorTensorProduct(typename _super::first_arg_type src1_arg, typename _super::second_arg_type src2_arg,
		       const Operation& op_arg=Operation())
      : _super(src1_arg, src2_arg, op_arg) { }

   using _super::dim;
};

template <typename VectorRef1, typename VectorRef2, typename Operation>
struct spec_object_traits< VectorTensorProduct<VectorRef1,VectorRef2,Operation> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true;
};

template <typename VectorRef1, typename VectorRef2, typename Operation, typename Feature>
struct check_container_feature<VectorTensorProduct<VectorRef1,VectorRef2,Operation>, Feature>
   : check_container_feature<ContainerProduct<VectorRef1, VectorRef2, Operation>, Feature> { };

/* -------------
 *  VectorChain
 * ------------- */

template <typename VectorRef1, typename VectorRef2>
class VectorChain
   : public ContainerChain<VectorRef1, VectorRef2>,
     public GenericVector< VectorChain<VectorRef1,VectorRef2>,
			   typename identical<typename deref<VectorRef1>::type::element_type,
					      typename deref<VectorRef2>::type::element_type>::type > {
   typedef ContainerChain<VectorRef1, VectorRef2> _super;
public:
   VectorChain(typename _super::first_arg_type src1_arg, typename _super::second_arg_type src2_arg)
      : _super(src1_arg,src2_arg) { }

   VectorChain& operator= (const VectorChain& v)
   {
      return VectorChain::generic_type::operator=(v);
   }
   using VectorChain::generic_type::operator=;
   using _super::dim;
};

template <typename VectorRef1, typename VectorRef2>
struct spec_object_traits< VectorChain<VectorRef1, VectorRef2> >
   : spec_object_traits< ContainerChain<VectorRef1, VectorRef2> > { };

template <typename VectorRef1, typename VectorRef2, typename Feature>
struct check_container_feature< VectorChain<VectorRef1, VectorRef2>, Feature>
   : check_container_feature< ContainerChain<VectorRef1, VectorRef2>, Feature> { };

/* ---------------------
 *  SingleElementVector
 * --------------------- */

template <typename E>
class SingleElementVector
   : public single_elem_container<E>,
     public GenericVector<SingleElementVector<E>, typename deref<E>::type> {
public:
   SingleElementVector(typename function_argument<E>::type arg)
      : single_elem_container<E>(arg) { }
};

template <typename E>
class SingleElementVector<E&>
   : public single_elem_container<E&>,
     public GenericVector<SingleElementVector<E&>, E> {
protected:
   ~SingleElementVector();
public:
   using SingleElementVector::generic_type::operator=;
};

template <typename E>
class SingleElementVector<const E&>
   : public single_elem_container<const E&>,
     public GenericVector<SingleElementVector<const E&>, E> {
protected:
   ~SingleElementVector();
};

template <typename E, typename Predicate=conv<E,bool> >
class SingleElementSparseVector
   : public single_elem_container<E, true>,
     public GenericVector<SingleElementSparseVector<E,Predicate>, E> {
public:
   SingleElementSparseVector(typename function_argument<E>::type arg)
      : single_elem_container<E, true>(arg) { }
   using single_elem_container<E, true>::dim;
};

template <typename E, typename Predicate>
class SingleElementSparseVector<E&, Predicate>
   : public single_elem_container<E&, true>,
     public GenericVector<SingleElementSparseVector<E&,Predicate>, E> {
protected:
   ~SingleElementSparseVector();
public:
   using single_elem_container<E&, true>::dim;
};

template <typename E, typename Predicate>
class SingleElementSparseVector<const E&, Predicate>
   : public single_elem_container<const E&, true>,
     public GenericVector<SingleElementSparseVector<const E&,Predicate>, E> {
protected:
   ~SingleElementSparseVector();
public:
   using single_elem_container<const E&, true>::dim;
};

template <typename E, typename Predicate>
struct check_container_feature<SingleElementSparseVector<E, Predicate>, pure_sparse> : True { };

template <typename ElementRef, typename IndexRef>
struct SingleElementSparseVector_factory {
   typedef ElementRef first_argument_type;
   typedef int second_argument_type;
   typedef SingleElementSparseVector<ElementRef> vector_type;
   typedef typename inherit_ref<vector_type, ElementRef>::type result_type;

   result_type operator() (typename function_argument<ElementRef>::type x, int) const
   {
      return static_cast<result_type>(vector_type::create(x));
   }
   // never called, but must be defined
   template <typename Iterator2>
   result_type operator() (operations::partial_left, typename function_argument<ElementRef>::type x, const Iterator2&) const
   {
      return static_cast<result_type>(vector_type::create(x));
   }

   template <typename Iterator1>
   result_type operator() (operations::partial_right, const Iterator1&, int) const
   {
      return static_cast<result_type>(vector_type::create_empty());
   }
};

template <typename E> inline
SingleElementVector<E&>&
scalar2vector(E& x)
{
   return reinterpret_cast<SingleElementVector<E&>&>(x);
}

template <typename E> inline
const SingleElementVector<const E&>&
scalar2vector(const E& x)
{
   return reinterpret_cast<const SingleElementVector<const E&>&>(x);
}

template <typename E> inline
const SingleElementSparseVector<const E&>&
scalar2sparse_vector(const E& x)
{
   conv<E, bool> non_zero;
   if (non_zero(x))
      return reinterpret_cast<const SingleElementSparseVector<const E&>&>(x);
   else
      return SingleElementSparseVector<const E&>::create_empty();
}

template <typename E, typename Predicate> inline
const SingleElementSparseVector<const E&, Predicate>&
scalar2sparse_vector(const E& x, Predicate non_zero)
{
   if (non_zero(x))
      return reinterpret_cast<const SingleElementSparseVector<const E&, Predicate>&>(x);
   else
      return SingleElementSparseVector<const E&, Predicate>::create_empty();
}

/* -------------------
 *  SameElementVector
 * ------------------- */

template <typename ElemRef>
class SameElementVector
   : public repeated_item_container<ElemRef>,
     public GenericVector<SameElementVector<ElemRef>,
			  typename object_traits<typename deref<ElemRef>::type>::persistent_type> {
   typedef repeated_item_container<ElemRef> _base;
public:
   SameElementVector(typename _base::arg_type value_arg, int dim_arg)
      : _base(value_arg,dim_arg) { }
};

template <typename ElemRef>
struct spec_object_traits< SameElementVector<ElemRef> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename E> inline
const SameElementVector<const E&>
same_element_vector(const E& x, int dim)
{
   return SameElementVector<const E&>(x, dim);
}

template <typename E> inline
const SameElementVector<E>
ones_vector(int dim)
{
   return SameElementVector<E>(E(1), dim);
}

template <typename E> inline
const SameElementVector<E>
zero_vector(int dim)
{
   return SameElementVector<E>(E(0), dim);
}

template <typename E, bool _inverse=false>
struct apparent_data_accessor {
   typedef void argument_type;
   typedef const E& result_type;

   apparent_data_accessor(const E* data_arg=0) : _data(data_arg) { }

   template <typename Iterator>
   result_type operator() (const Iterator&) const
   {
      return *_data;
   }

   // these methods are for construct_dense
   typedef void first_argument_type;
   typedef void second_argument_type;

   template <typename Iterator1, typename Iterator2>
   result_type operator() (const Iterator1&, const Iterator2&) const
   {
      if (_inverse) return zero();
      return *_data;
   }
   template <typename Iterator1, typename Iterator2>
   result_type operator() (operations::partial_left, const Iterator1&, const Iterator2&) const
   {
      if (_inverse) return zero();
      return *_data;
   }
   template <typename Iterator1, typename Iterator2>
   result_type operator() (operations::partial_right, const Iterator1&, const Iterator2&) const
   {
      if (_inverse) return *_data;
      return zero();
   }
protected:
   const E* _data;
   operations::clear<E> zero;
};

template <typename SetRef, typename ElemRef>
class SameElementSparseVector
   : public modified_container_impl< SameElementSparseVector<SetRef,ElemRef>,
				     list( Container< typename Set_with_dim_helper<SetRef>::container >,
					   Operation< pair<apparent_data_accessor<typename deref<ElemRef>::type, complement_helper<SetRef>::value>,
					                   operations::identity<int> > > ) >,
     public GenericVector< SameElementSparseVector<SetRef,ElemRef>,
			   typename object_traits<typename deref<ElemRef>::type>::persistent_type > {
   typedef modified_container_impl<SameElementSparseVector> _super;
protected:
   typedef Set_with_dim_helper<SetRef> helper;
   typename helper::alias_type set;
   alias<ElemRef> apparent_elem;

public:
   typedef typename helper::alias_type::arg_type first_arg_type;
   typedef typename alias<ElemRef>::arg_type second_arg_type;

   SameElementSparseVector(first_arg_type set_arg, second_arg_type data_arg, int dim_arg=0)
      : set(helper::create(set_arg,dim_arg)), apparent_elem(data_arg) { }

   const typename _super::container& get_container() const
   {
      return set;
   }
   typename alias<ElemRef>::const_reference get_apparent_element() const
   {
      return apparent_elem;
   }
   typename _super::operation get_operation() const
   {
      return typename _super::operation(&get_apparent_element(), operations::identity<int>());
   }
   using _super::dim;
};

template <typename SetRef, typename ElemRef>
struct check_container_feature<SameElementSparseVector<SetRef,ElemRef>, pure_sparse> : True { };

template <typename SetRef, typename ElemRef>
class construct_dense< SameElementSparseVector<SetRef,ElemRef> >
   : public modified_container_pair_impl< construct_dense< SameElementSparseVector<SetRef,ElemRef> >,
					  list( Container1< typename complement_helper<SetRef>::base_type >,
						Container2< sequence >,
						IteratorCoupler< zipping_coupler< operations::cmp<>, set_union_zipper> >,
						Operation< pair< typename SameElementSparseVector<SetRef,ElemRef>::operation::first_type,
						                 BuildBinaryIt<operations::zipper> > >,
						MasqueradedTop ) > {
   typedef modified_container_pair_impl<construct_dense> _super;
public:
   const typename _super::container1& get_container1() const
   {
      return Set_with_dim_helper<SetRef>::get_set(this->hidden().get_container());
   }

   sequence get_container2() const
   {
      return sequence(0, size());
   }

   typename _super::operation get_operation() const
   {
      return typename _super::operation(&(this->hidden().get_apparent_element()), BuildBinaryIt<operations::zipper>());
   }

   int size() const
   {
      return this->hidden().dim();
   }
};

template <typename E, typename Vector> inline
const SameElementSparseVector<const Vector&, E>
same_element_sparse_vector(const GenericSet<Vector,int>& s, int dim)
{
   return SameElementSparseVector<const Vector&, E>(s.top(), E(1), dim);
}

template <typename E, typename Vector> inline
const SameElementSparseVector<const Vector&, const E&>
same_element_sparse_vector(const GenericSet<Vector,int>& s, const E& x, int dim)
{
   return SameElementSparseVector<const Vector&, const E&>(s.top(), x, dim);
}

template <typename E, typename Vector> inline
const SameElementSparseVector<const Complement<Vector>&, E>
same_element_sparse_vector(const Complement<Vector,int>& s, int dim)
{
   return SameElementSparseVector<const Complement<Vector>& ,E>(s, E(1), dim);
}

template <typename E, typename Vector> inline
const SameElementSparseVector<const Complement<Vector>&, const E&>
same_element_sparse_vector(const Complement<Vector,int>& s, const E& x, int dim)
{
   return SameElementSparseVector<const Complement<Vector>&, const E&>(s, x, dim);
}

/* ------------
 *  UnitVector
 * ------------ */

template <typename ElemRef=void>
struct UnitVector_factory {
protected:
   int dim;
public:
   typedef int first_argument_type;
   typedef ElemRef second_argument_type;
   typedef SameElementSparseVector<masquerade<SingleElementSet, int>, ElemRef> vector_type;
   typedef const vector_type result_type;

   UnitVector_factory(int dim_arg=0) : dim(dim_arg) { }

   result_type operator() (first_argument_type pos,
			   typename function_argument<ElemRef>::type elem) const
   {
      return vector_type(pos,elem,dim);
   }
};

template <>
struct UnitVector_factory<void> : operations::incomplete {
protected:
   int dim;
public:
   UnitVector_factory(int dim_arg=0) : dim(dim_arg) { }
   operator int () const { return dim; }
};

template <typename Iterator1, typename Iterator2, typename Reference1, typename Reference2>
struct binary_op_builder< UnitVector_factory<void>, Iterator1, Iterator2, Reference1, Reference2 > {
   typedef UnitVector_factory<Reference2> operation;

   static operation create(int dim_arg) { return operation(dim_arg); }
};

template <typename E> inline
const SameElementSparseVector<masquerade<SingleElementSet, int>, const E&>
unit_vector(int dim, int i, const E& x)
{
   DIMENSION_CHECK(dim<=0, "unit_vector - dimension should be positive");
   INDEX_CHECK((i<0 || i>=dim), "unit_vector - index of non-zero element out of range");
   return SameElementSparseVector<masquerade<SingleElementSet, int>, const E&>(i,x,dim);
}

template <typename E> inline
const SameElementSparseVector<masquerade<SingleElementSet, int>, E>
unit_vector(int dim, int i)
{
   DIMENSION_CHECK(dim<=0, "unit_vector - dimension should be positive");
   INDEX_CHECK((i<0 || i>=dim), "unit_vector - index of non-zero element out of range");
   return SameElementSparseVector<masquerade<SingleElementSet, int>, E>(i,E(1),dim);
}

/* ----------------
 *  ExpandedVector
 * ---------------- */

template <typename VectorRef>
class ExpandedVector
   : public TransformedContainer<masquerade_add_features<VectorRef, sparse_compatible>,
				 pair<nothing, operations::fix2<int, operations::composed12< BuildUnaryIt<operations::index2element>, void, BuildBinary<operations::add> > > > >,
     public GenericVector< ExpandedVector<VectorRef>, typename deref<VectorRef>::type::element_type > {
   typedef TransformedContainer<masquerade_add_features<VectorRef, sparse_compatible>,
				pair<nothing, operations::fix2<int, operations::composed12< BuildUnaryIt<operations::index2element>, void, BuildBinary<operations::add> > > > >
      _super;
   int _dim;
public:
   ExpandedVector(typename _super::arg_type src_arg, int offset, int dim_arg)
      : _super(src_arg, offset), _dim(dim_arg) { }

   int dim() const { return _dim; }
};

template <typename VectorRef>
struct check_container_feature<ExpandedVector<VectorRef>, sparse> : True { };

template <typename VectorRef>
struct check_container_feature<ExpandedVector<VectorRef>, pure_sparse>
   : check_container_ref_feature<VectorRef, pure_sparse> { };

template <typename VectorRef=void>
struct ExpandedVector_factory {
protected:
   int offset, dim;
public:
   typedef VectorRef argument_type;
   typedef ExpandedVector<VectorRef> vector_type;
   typedef const vector_type result_type;

   ExpandedVector_factory(int offset_arg=0, int dim_arg=0)
      : offset(offset_arg), dim(dim_arg) { }

   template <typename Ref2>
   ExpandedVector_factory(const ExpandedVector_factory<Ref2>& f)
      : offset(f.offset), dim(f.dim) { }

   result_type operator() (typename function_argument<VectorRef>::type vec) const
   {
      return vector_type(vec,offset,dim);
   }

   template <typename> friend class ExpandedVector_factory;
};

template <>
struct ExpandedVector_factory<void> : operations::incomplete {
protected:
   int offset, dim;
public:
   ExpandedVector_factory(int offset_arg=0, int dim_arg=0)
      : offset(offset_arg), dim(dim_arg) { }
   template <typename> friend class ExpandedVector_factory;
};

template <typename Iterator, typename Reference>
struct unary_op_builder< ExpandedVector_factory<void>, Iterator, Reference > {
   typedef ExpandedVector_factory<Reference> operation;

   template <typename Ref2>
   static operation create(const ExpandedVector_factory<Ref2>& f) { return f; }
};

/* ---------------------------------------
 *  functor objects for vector operations
 * --------------------------------------- */

namespace operations {

template <typename OpRef>
struct neg_impl<OpRef, is_vector> {
   typedef OpRef argument_type;
   typedef const LazyVector1<typename attrib<OpRef>::plus_const, BuildUnary<neg> > result_type;

   result_type operator() (typename function_argument<OpRef>::const_type x) const { return x; }

   typename function_argument<OpRef>::lvalue
   assign(typename function_argument<OpRef>::lvalue x) const
   {
      return static_cast<typename function_argument<OpRef>::lvalue> (inv_sign(x));
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct add_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<typename attrib<LeftRef>::plus_const, typename attrib<RightRef>::plus_const,
			     BuildBinary<add> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() != r.dim()), "operator+ - vector dimension mismatch");
      return result_type(l,r);
   }

   template <typename Iterator2>
   typename function_argument<LeftRef>::const_type
   operator() (partial_left, typename function_argument<LeftRef>::const_type l, const Iterator2&) const
   {
      return l;
   }

   template <typename Iterator1>
   typename function_argument<RightRef>::const_type
   operator() (partial_right, const Iterator1&, typename function_argument<RightRef>::const_type r) const
   {
      return r;
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l+=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct sub_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<typename attrib<LeftRef>::plus_const, typename attrib<RightRef>::plus_const,
			     BuildBinary<sub> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() != r.dim()), "operator- - vector dimension mismatch");
      return result_type(l,r);
   }

   template <typename Iterator2>
   typename function_argument<LeftRef>::const_type
   operator() (partial_left, typename function_argument<LeftRef>::const_type l, const Iterator2&) const
   {
      return l;
   }
   template <typename Iterator1>
   typename neg<RightRef>::result_type
   operator() (partial_right, const Iterator1&, typename function_argument<RightRef>::const_type r) const
   {
      neg<RightRef> n;
      return n(r);
   }
   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l-=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_vector, is_scalar> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<typename attrib<LeftRef>::plus_const,
			     pm::masquerade_keep_ref<unlimited_constant_container, typename prevent_lazy<RightRef>::type>,
			     BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const {
      return result_type(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l*=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_scalar, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<pm::masquerade_keep_ref<unlimited_constant_container, typename prevent_lazy<LeftRef>::type>,
			     typename attrib<RightRef>::plus_const,
			     BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef typename add<typename deref<LeftRef>::type::value_type,
			typename deref<RightRef>::type::value_type>::result_type result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() != r.dim()), "operator* - vector dimension mismatch");
      return accumulate(attach_operation(l,r,BuildBinary<mul>()), BuildBinary<add>());
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct tensor_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const VectorTensorProduct<typename attrib<LeftRef>::plus_const,
				     typename attrib<RightRef>::plus_const, BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_vector, is_scalar> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<typename attrib<LeftRef>::plus_const,
			     pm::masquerade_keep_ref<unlimited_constant_container, typename prevent_lazy<RightRef>::type>,
			     BuildBinary<div> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l/=r);
   }
};

template <typename OpRef>
struct square_impl<OpRef, is_vector> {
   typedef OpRef argument_type;
   typedef typename mul_impl<OpRef, OpRef, false, cons<is_vector, is_vector> >::result_type
      result_type;

   result_type operator() (typename function_argument<OpRef>::const_type x) const
   {
      return accumulate(attach_operation(x, BuildUnary<square>()), BuildBinary<add>());
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct concat_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef VectorChain<LeftRef, RightRef> chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      return chain(l,r);
   }
   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l|=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct concat_impl<LeftRef, RightRef, check, cons<is_vector, is_scalar> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   // provisions for the type-heterogeneous case, where the scalar value needs to be converted first
   typedef typename deref<LeftRef>::type::element_type element_type;
   typedef VectorChain<LeftRef, pm::masquerade_keep_ref<SingleElementVector,
							typename if_else<identical<element_type, typename deref<RightRef>::type>::value,
									 RightRef, const element_type>::type> >
      chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      return chain(l,r);
   }
   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l|=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct concat_impl<LeftRef, RightRef, check, cons<is_scalar, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   // provisions for the type-heterogeneous case, where the scalar value needs to be converted first
   typedef typename deref<RightRef>::type::element_type element_type;
   typedef VectorChain<pm::masquerade_keep_ref<SingleElementVector,
					       typename if_else<identical<element_type, typename deref<LeftRef>::type>::value, 
								LeftRef, const element_type>::type>,
		       RightRef>
      chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      return chain(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> >
   : concat_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > { };

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_vector, is_scalar> >
   : concat_impl<LeftRef, RightRef, check, cons<is_vector, is_scalar> > { };

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_scalar, is_vector> >
   : concat_impl<LeftRef, RightRef, check, cons<is_scalar, is_vector> > { };

} // end namespace operations
namespace operators {

template <typename Vector1, typename Vector2> inline
typename operations::add_impl<const Vector1&, const Vector2&, true>::result_type
operator+ (const GenericVector<Vector1>& l, const GenericVector<Vector2>& r)
{
   operations::add_impl<const Vector1&, const Vector2&, true> op;
   return op(l.top(), r.top());
}

template <typename Vector1, typename Vector2> inline
typename operations::sub_impl<const Vector1&, const Vector2&, true>::result_type
operator- (const GenericVector<Vector1>& l, const GenericVector<Vector2>& r)
{
   operations::sub_impl<const Vector1&, const Vector2&, true> op;
   return op(l.top(), r.top());
}

template <typename Vector1, typename Vector2> inline
bool operator== (const GenericVector<Vector1>& l, const GenericVector<Vector2>& r)
{
   DIMENSION_CHECK((l.dim() != r.dim()), "operator== - vector dimension mismatch");
   operations::eq<const Vector1&, const Vector2&> op;
   return op(l.top(), r.top());
}

template <typename Vector1, typename Vector2> inline
bool operator< (const GenericVector<Vector1>& l, const GenericVector<Vector2>& r)
{
   DIMENSION_CHECK((l.dim() != r.dim()), "operator< - vector dimension mismatch");
   operations::lt<const Vector1&, const Vector2&> op;
   return op(l.top(), r.top());
}

} // end namespace operators
} // end namespace pm

namespace polymake {
   using pm::GenericVector;
   using pm::FixedVector;
   using pm::scalar2vector;
   using pm::same_element_vector;
   using pm::same_element_sparse_vector;
   using pm::ones_vector;
   using pm::zero_vector;
   using pm::unit_vector;
}

namespace std {
   template <typename Vector1, typename Vector2, typename E> inline
   void swap(pm::GenericVector<Vector1,E>& v1, pm::GenericVector<Vector2,E>& v2)
   {
      v1.top().swap(v2.top());
   }

   template <typename Vector1, typename Vector2, typename E> inline
   void swap(const pm::GenericVector<Vector1,E>& v1, const pm::GenericVector<Vector2,E>& v2)
   {
      pm::non_const(v1.top()).swap(pm::non_const(v2.top()));
   }

   // due to silly overloading rules
   template <typename VectorRef1, typename VectorRef2> inline
   void swap(pm::VectorChain<VectorRef1,VectorRef2>& v1, pm::VectorChain<VectorRef1,VectorRef2>& v2)
   {
      v1.swap(v2);
   }

   template <typename E, int _size> inline
   void swap(pm::FixedVector<E,_size>& v1, pm::FixedVector<E,_size>& v2)
   {
      v1.swap(v2);
   }
}

namespace std_ext {
   template <typename T> struct hash;

   template <typename E>
   struct hash< pm::Vector<E> > {
      template <typename _Vector>
      size_t operator() (const pm::GenericVector<_Vector,E>& v) const
      {
	 hash<E> element_hasher;
	 size_t a=1;
	 for (typename pm::ensure_features<_Vector, pm::cons<pm::end_sensitive, pm::sparse_compatible> >::const_iterator e=ensure(v.top(), (pm::cons<pm::end_sensitive, pm::sparse_compatible>*)0).begin(); !e.at_end(); ++e)
	    a=a*element_hasher(*e)+e.index();
	 return a;
      }
   };

   template <typename E, typename Predicate>
   struct hash< pm::SparseVector<E, Predicate> > : hash< pm::Vector<E> > { };
   template <typename E, int _size>
   struct hash< pm::FixedVector<E,_size> > : hash< pm::Vector<E> > { };
}

#endif // _POLYMAKE_GENERIC_VECTOR_H

// Local Variables:
// mode:C++
// End:
