/* 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_PERMUTATIONS_H
#define _POLYMAKE_PERMUTATIONS_H "$Project: polymake $$Id: permutations.h,v 1.10 2006/03/09 16:21:36 gawrilow Exp $"

#include <vector>
#include <GenericMatrix.h>
#include <algorithm>
#include <Integer.h>
#include <Map.h>
#include <Array.h>
#include <sstream>

namespace pm {

template <typename Permutation>
int permutation_sign (const Permutation& v)
{
   const int l=v.size();
   if (l<=1) return 1;

   std::vector<int> w(v.size());
   std::copy(v.begin(), entire(w));

   int sign=1;
   for (int i=0; i<l; ) {
      if (w[i]!=i) {
	 int k=w[i];
	 w[i]=w[k];
	 w[k]=k;
	 sign = -sign;
      } else {
	 ++i;
      }
   }
   return sign;
}

template <typename Permutation, typename InvPermutation> inline
void inverse_permutation (const Permutation& perm, InvPermutation& inv_perm)
{
   inv_perm.resize(perm.size());
   int pos=0;
   for (typename Entire<Permutation>::const_iterator i=entire(perm);  !i.at_end();  ++i, ++pos)
      inv_perm[*i]=pos;
}

template <typename Input1, typename Input2, typename Output, typename Comparator>
void find_permutation (Input1 src1, Input2 src2, Output dst, const Comparator& comparator)
{
   typedef Map<typename iterator_traits<Input1>::value_type, int, Comparator> map_type;
   map_type index_map(comparator);

   for (int i=0; !src1.at_end(); ++src1, ++i)
      index_map[*src1]=i;

   for (; !src2.at_end(); ++src2, ++dst) {
      typename map_type::iterator where=index_map.find(*src2);
      if (where.at_end()) {
	 std::string reason;
	 if (index_map.empty()) {
	    reason="not a permutation: first sequence is shorter";
	 } else {
	    std::ostringstream os;
	    wrap(os) << "not a permutation: no match for <" << *src2 << ">";
	    reason=os.str();
	 }
	 throw no_match(reason);
      }
      *dst=where->second;
      index_map.erase(where);
   }

   if (!index_map.empty())
      throw no_match("not a permutation: second sequence is shorter");
}

template <typename Container1, typename Container2, typename Comparator> inline
Array<int> find_permutation (const Container1& c1, const Container2& c2, const Comparator& comparator)
{
   Array<int> perm(c1.size());
   find_permutation(entire(c1), entire(c2), perm.begin(), comparator);
   return perm;
}

template <typename Container1, typename Container2> inline
Array<int> find_permutation (const Container1& c1, const Container2& c2)
{
   return find_permutation(c1,c2,polymake::operations::cmp());
}


enum permutation_sequence {
   permutations_heap,
   permutations_adjacent,
   permutations_lex
};

template <permutation_sequence kind=permutations_heap>
class AllPermutations;

template <permutation_sequence kind>
class permutation_iterator;

class permutation_iterator_base {
public:
   typedef forward_iterator_tag iterator_category;
   typedef std::vector<int> value_type;
   typedef const value_type* pointer;
   typedef const value_type& reference;
   typedef ptrdiff_t difference_type;
protected:
   value_type perm;

   permutation_iterator_base() { }
   permutation_iterator_base(int n) : perm(n) { reset(n); }
   void reset(int n) { std::copy(entire(sequence(0,n)), perm.begin()); }

public:
   reference operator* () const { return perm; }
   pointer operator-> () const { return &perm; }
};

/// Implementation of the Heap's algorithm by R. Sedgewick
template <>
class permutation_iterator<permutations_heap>
   : public permutation_iterator_base {
   friend class AllPermutations<permutations_heap>;
public:
   typedef permutation_iterator iterator;
   typedef iterator const_iterator;

protected:
   std::vector<int> cnt;
   int n, pos;

   permutation_iterator(int n_arg, bool)
      : permutation_iterator_base(n_arg),
	cnt(n_arg,0), n(n_arg), pos(n_arg) { }

public:
   permutation_iterator()
      : n(0), pos(0) { }

   explicit permutation_iterator(int n_arg)
      : permutation_iterator_base(n_arg),
	cnt(size_t(n_arg),0),
	n(n_arg), pos(n_arg>1) { }

   iterator& operator++ ()
   {
      do {
	 if (cnt[pos]<pos) {
	    std::swap(perm[pos], perm[pos%2 * cnt[pos]]);
	    ++cnt[pos];
	    pos=1;
	    break;
	 }
	 cnt[pos]=0;
      } while (++pos < n);
      return *this;
   }
   const iterator operator++ (int) { iterator copy(*this); operator++(); return copy; }

   bool at_end() const { return pos>=n; }

   bool operator== (const iterator& it) const
   {
      return n==it.n && pos==it.pos && std::equal(entire(cnt), it.cnt.begin());
   }
   bool operator!= (const iterator& it) const
   {
      return !operator==(it);
   }

   void reset()
   {
      permutation_iterator_base::reset(n);
      pos=1;
      std::fill(entire(cnt), 0);
   }
};

template <>
class permutation_iterator<permutations_adjacent>
   : public permutation_iterator_base {
   friend class AllPermutations<permutations_adjacent>;
public:
   typedef permutation_iterator iterator;
   typedef iterator const_iterator;

protected:
   std::vector<int> move;
   int n, next;

   void incr()
   {
      if (next==0) {
	 next=-1; return;
      }
      int lim=n, offset=0;
      std::vector<int>::iterator pos=move.begin();
      for (;;) {
	 int m=*pos;
	 if (m>0) {
	    if (m<lim) {
	       next=m+offset;
	       ++*pos;
	       break;
	    }
	 } else {
	    if (m<-1) {
	       next=offset-m-1;
	       ++*pos;
	       break;
	    }
	    ++offset;
	 }
	 if (--lim > 1) {
	    *pos=-m;
	    ++pos;
	 } else {
	    next=0;
	    break;
	 }
      }
   }

public:
   permutation_iterator() : n(0) { }

   explicit permutation_iterator(int n_arg)
      : permutation_iterator_base(n_arg),
	move(size_t(n_arg),1), n(n_arg), next(n-1)
   {
      if (next>0) {
	 next=1; move[0]=2;
      }
   }

   iterator& operator++ ()
   {
      std::swap(perm[next-1], perm[next]);
      incr();
      return *this;
   }
   const iterator operator++ (int) { iterator copy(*this); operator++(); return copy; }

   int next_swap() const { return next; }

   bool at_end() const { return next<0; }

   bool operator== (const iterator& it) const
   {
      return next==it.next && n==it.n && (!next || std::equal(entire(move), it.move.begin()));
   }
   bool operator!= (const iterator& it) const
   {
      return !operator==(it);
   }

   void reset()
   {
      permutation_iterator_base::reset(n);
      std::fill(entire(move),1);
      next=n-1;
      if (next>1) {
	 next=1; move[0]=2;
      }
   }
};

// permutations in lexicographic order via STL
template <>
class permutation_iterator<permutations_lex>
   : public permutation_iterator_base {
   friend class AllPermutations<permutations_lex>;
public:
   typedef permutation_iterator iterator;
   typedef iterator const_iterator;

protected:
   bool _at_end;

public:
   permutation_iterator() { }

   explicit permutation_iterator(int n, bool end_arg=false)
      : permutation_iterator_base(n), _at_end(end_arg || n==0) { }

   iterator& operator++ ()
   {
      _at_end=!std::next_permutation(perm.begin(), perm.end());
      return *this;
   }
   const iterator operator++ (int) { iterator copy(*this); operator++(); return copy; }

   bool at_end() const { return _at_end; }

   bool operator== (const iterator& it) const
   {
      return _at_end ? it._at_end : !it._at_end && perm.size()==it.perm.size() && std::equal(entire(perm), it.perm.begin());
   }
   bool operator!= (const iterator& it) const
   {
      return !operator==(it);
   }

   void reset()
   {
      permutation_iterator_base::reset(perm.size());
      _at_end=perm.empty();
   }
};

template <permutation_sequence kind>
struct check_iterator_feature<permutation_iterator<kind>, end_sensitive> : True { };
template <permutation_sequence kind>
struct check_iterator_feature<permutation_iterator<kind>, resettable> : True { };

template <permutation_sequence kind>
class AllPermutations {
protected:
   int n;
public:
   explicit AllPermutations(int n_arg=0) : n(n_arg) { }

   typedef permutation_iterator<kind> iterator;
   typedef iterator const_iterator;
   typedef typename iterator::value_type value_type;
   typedef typename iterator::reference reference;
   typedef reference const_reference;

   iterator begin() const { return iterator(n); }
   iterator end() const { return iterator(n,true); }

   value_type front() const
   {
      return value_type(sequence(0,n).begin(), sequence(0,n).end());
   }

   Integer size() const { return n ? Integer::fac(n) : Integer(0); }
   bool empty() const { return n==0; }
};

template <permutation_sequence kind>
struct spec_object_traits< AllPermutations<kind> >
   : spec_object_traits<is_container> {
   static const bool is_always_const=true;
};

inline
AllPermutations<> all_permutations(int n)
{
   return AllPermutations<>(n);
}

template <permutation_sequence kind>
AllPermutations<kind> all_permutations(int n)
{
   return AllPermutations<kind>(n);
}


template <typename PermutationRef, typename Element=int>
class PermutationMatrix
: public GenericMatrix< PermutationMatrix<PermutationRef,Element>, Element > {
protected:
   alias<PermutationRef> perm;
   mutable std::vector<int> inv_perm;
public:
   typedef Element value_type;
   typedef const Element& reference;
   typedef reference const_reference;

   typedef typename alias<PermutationRef>::arg_type arg_type;
   PermutationMatrix(arg_type perm_arg) : perm(perm_arg) { }

   typename alias<PermutationRef>::const_reference get_perm() const { return perm; }

   const std::vector<int>& get_inv_perm() const
   {
      if (inv_perm.empty() && !get_perm().empty())
	 inverse_permutation(get_perm(), inv_perm);
      return inv_perm;
   }

   const Element& one() const
   {
      static const Element One(1);
      return One;
   }
};

template <typename PermutationRef, typename Element>
struct spec_object_traits< PermutationMatrix<PermutationRef,Element> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename PermutationRef, typename Element>
struct check_container_feature< PermutationMatrix<PermutationRef,Element>, pure_sparse> : True { };

template <typename PermutationRef, typename Element>
class matrix_random_access_methods< PermutationMatrix<PermutationRef,Element> > {
   typedef PermutationMatrix<PermutationRef,Element> master;
public:
   const Element& operator() (int i, int j) const
   {
      const master& me=*static_cast<const master*>(this);
      INDEX_CHECK((i<0 || i>=me.rows()), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.rows()), "operator() - col index out of range");
      operations::clear<Element> zero;
      return me.get_perm()[i]==j ? me.one() : zero();
   }
};

template <typename PermutationRef, typename Element>
class Rows< PermutationMatrix<PermutationRef,Element> >
   : public modified_container_pair_impl< Rows< PermutationMatrix<PermutationRef,Element> >,
	       cons< Container1< PermutationRef >,
	       cons< Container2< unlimited_constant_container<const Element&> >,
	       cons< Accessor< UnitVector_factory<> >,
		     Hidden< PermutationMatrix<PermutationRef,Element> > > > >
            > {
   typedef modified_container_pair_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container1& get_container1() const
   {
      return this->hidden().get_perm();
   }
   const typename _super::container2& get_container2() const
   {
      return constant(this->hidden().one());
   }
   typename _super::accessor get_accessor() const
   {
      return this->size();
   }
};

template <typename PermutationRef, typename Element>
class Cols< PermutationMatrix<PermutationRef,Element> >
   : public modified_container_pair_impl< Cols< PermutationMatrix<PermutationRef,Element> >,
	       cons< Container1< std::vector<int> >,
	       cons< Container2< unlimited_constant_container<const Element&> >,
	       cons< Accessor< UnitVector_factory<> >,
		     Hidden< PermutationMatrix<PermutationRef,Element> > > > >
            > {
   typedef modified_container_pair_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::container1& get_container1() const
   {
      return this->hidden().get_inv_perm();
   }
   const typename _super::container2& get_container2() const
   {
      return constant(this->hidden().one());
   }
   typename _super::accessor get_accessor() const
   {
      return this->size();
   }
};

template <typename Element, typename Permutation> inline
PermutationMatrix<const Permutation&, Element>
permutation_matrix(const Permutation& perm)
{
   return perm;
}

} // end namespace pm

namespace polymake {
   using pm::permutation_sign;
   using pm::inverse_permutation;
   using pm::find_permutation;
   using pm::permutation_sequence;
   using pm::permutations_heap;
   using pm::permutations_adjacent;
   using pm::permutations_lex;
   using pm::AllPermutations;
   using pm::all_permutations;
   using pm::PermutationMatrix;
   using pm::permutation_matrix;
}

#endif // _POLYMAKE_PERMUTATIONS_H

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