/* Copyright (c) 1997-2004
   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.3 2004/09/20 15:11:51 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());
}


class AllPermutations;

/// Implementation of the Heap's algorithm by R. Sedgewick
class Permutation_iterator {
   friend class AllPermutations;
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;
   typedef Permutation_iterator iterator;
   typedef iterator const_iterator;

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

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

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

   explicit Permutation_iterator(int n_arg)
      : perm(sequence(0,n_arg).begin(), sequence(0,n_arg).end()),
	cnt(size_t(n_arg),0),
	n(n_arg), pos(n_arg>1) { }

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

   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 operator== (const iterator& it) const
   {
      return n==it.n && std::equal(entire(cnt), it.cnt.begin());
   }
   bool operator!= (const iterator& it) const
   {
      return !operator==(it);
   }

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

   void reset()
   {
      pos=1;
      std::copy(sequence(0,n).begin(), entire(perm));
      std::fill(entire(cnt), 0);
   }
};

template <>
struct check_iterator_feature<Permutation_iterator, end_sensitive> : True { };
template <>
struct check_iterator_feature<Permutation_iterator, resettable> : True { };

class AllPermutations {
protected:
   int n;
public:
   explicit AllPermutations(int n_arg=0) : n(n_arg) { }

   typedef Permutation_iterator iterator;
   typedef iterator const_iterator;
   typedef iterator::value_type value_type;
   typedef iterator::reference reference;
   typedef reference const_reference;

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

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

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

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

inline
AllPermutations all_permutations(int n)
{
   return AllPermutations(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;
   }
};

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;
      static const Element One(1);
      return me.get_perm()[i]==j ? One : zero();
   }
};

template <typename PermutationRef, typename Element>
class Rows< PermutationMatrix<PermutationRef,Element> >
   : public modified_container_impl< Rows< PermutationMatrix<PermutationRef,Element> >,
	       cons< Container< PermutationRef >,
	       cons< Accessor< UnitVector_factory<Element> >,
		     Hidden< PermutationMatrix<PermutationRef,Element> > > >
            > {
   typedef modified_container_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container& get_container() const
   {
      return this->hidden().get_perm();
   }
   typename _super::accessor get_accessor() const
   {
      return this->size();
   }
};

template <typename PermutationRef, typename Element>
class Cols< PermutationMatrix<PermutationRef,Element> >
   : public modified_container_impl< Cols< PermutationMatrix<PermutationRef,Element> >,
	       cons< Container< std::vector<int> >,
	       cons< Accessor< UnitVector_factory<Element> >,
		     Hidden< PermutationMatrix<PermutationRef,Element> > > >
            > {
   typedef modified_container_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::container& get_container() const
   {
      return this->hidden().get_inv_perm();
   }
   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::AllPermutations;
   using pm::all_permutations;
   using pm::PermutationMatrix;
   using pm::permutation_matrix;
}

#endif // _POLYMAKE_PERMUTATIONS_H

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