// GNU General Public License Agreement
// Copyright (C) 2004-2010 CodeCogs, Zyba Ltd, Broadwood, Holford, TA5 1DU, England.
//
// 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 CodeCogs. 
// You must retain a copy of this licence in all copies. 
//
// 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.
// ---------------------------------------------------------------------------------
//! Calculates the inverse of the given permutation.

#ifndef MATHS_COMBINATORICS_PERMUTATIONS_INVERSE_H
#define MATHS_COMBINATORICS_PERMUTATIONS_INVERSE_H

#include <assert.h>
#include <vector>

namespace Maths
{

namespace Combinatorics
{

namespace Permutations
{

//! Calculates the inverse of the given permutation.

std::vector<int>inverse(int n, int *p)
{
  std::vector<int> inv(p, p + n);

  for (int i = 1; i <= n; inv[i - 1] *= -1, i++)
    for (int j = inv[i - 1]; j > i; inv[j - 1] *= -1, j = -inv[j - 1]);
  for (int i = 1; i <= n; i++)
  {
    int j = -inv[i - 1];
    if (j >= 0)
    {
      int k = i;
      while (true)
      {
        int aux = inv[j - 1];
        inv[j - 1] = k;
        if (aux < 0) break;
        k = j, j = aux;
      }
    }
  }

  return inv;
}

}

}

}

#endif

#include <algorithm>
#include <iostream>

  template
  < typename T
  >
  std::ostream&
operator<<
  ( std::ostream& sout
  , std::vector<T> const& a_container
  )
{
   sout<<"{ ";
   unsigned const size=a_container.size();
   for(unsigned i=0; i<size; ++i)
   {
      if(0<i) sout<<", ";
      sout<<a_container[i];
   }
   sout<<"} ";
   return sout;
}

using namespace Maths::Combinatorics::Permutations;

int main()
{
    int const rank=5;
    typedef std::vector<int> permut_t;
    permut_t permut_o(rank);
    std::iota( permut_o.begin(), permut_o.end(), 1);
    {
        permut_t permut_p(permut_o);
        unsigned count=0;
        do
        {
            std::cout<<":count="<<count<<"\n";
            std::cout<<":permut_p="<<permut_p<<"\n";
            permut_t permut_i(rank);
            permut_i=inverse(rank,permut_p.data());
          #if 1
            permut_t permut_pi(rank);
            for( int i=0; i<rank; ++i)
            {
                int pi=permut_p[permut_i[i]-1];
                if(permut_o[i]!=pi)
                {
                    std::cout<<":inversion failed:permut_p[permut_i["<<i<<"]]="<<pi<<"\n";
                }
            }
          #endif
        }while 
          (  ++count< 130
          && std::next_permutation
             ( permut_p.begin()
             , permut_p.end()
             )
          );
    }
    return 0;
}