/* 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_SMITH_NORMAL_FORM_H
#define _POLYMAKE_SMITH_NORMAL_FORM_H "$Project: polymake $$Id: Smith_normal_form.h 7315 2006-04-02 21:37:53Z gawrilow $"

#include <SparseMatrix.h>
#include <Bitset.h>
#include <list>

namespace polymake { namespace topaz {

template <typename E>
class dummy_companion_logger {
public:
   void from_right(const SparseMatrix2x2<E>&) { }
   void from_left(const SparseMatrix2x2<E>&) { }
};

template <typename Logger> class TransposedLogger;

template <typename Logger> inline
TransposedLogger<Logger>& transpose_logger(Logger& l)
{
   return reinterpret_cast<TransposedLogger<Logger>&>(l);
}

template <typename E>
class TransposedLogger< dummy_companion_logger<E> > : public dummy_companion_logger<E> { };

template <typename Iterator>
typename pm::object_traits< typename pm::iterator_traits<Iterator>::value_type >::persistent_type
gcd_of_sequence(Iterator src)
{
   typename pm::object_traits< typename pm::iterator_traits<Iterator>::value_type >::persistent_type g=std::abs(*src);
   while (g>1 && (++src, !src.at_end()))
      g=gcd(g,*src);
   return g;
}

#ifdef POLY_DEBUG
template <typename Output, typename Matrix> inline
void dump_matrix(Output& os, const GenericMatrix<Matrix>& M)
{
   os << std::setw(6) << M;
}
#endif

template <typename _Matrix, typename CompanionLogger>
int Smith_normal_form_steps(_Matrix& M, CompanionLogger& Logger
#ifdef POLY_DEBUG
			    , pm::SparseMatrixStatistics<typename _Matrix::element_type>& stat,
			    bool do_dump, const char *line_name
#endif
			    )
{
   const bool Logger_dummy=pm::derived_from_instance<CompanionLogger, dummy_companion_logger>::value;
   typedef typename _Matrix::element_type E;
   E pivot_elem(1), g(1), g2(1);	// these are working variables in the following loops,
   SparseMatrix2x2<E> U;                // but they shouldn't be created and destroyed all over the time

   int r=0, c=0, skipped=0;
   const int R=M.rows();
   while (r<R) {
      bool can_reduce=true, can_eliminate=false;
      typename Entire<typename _Matrix::row_type>::iterator e=entire(M.row(r));
      if (e.at_end()  ||  M.row(r).size()==1 && M.col(e.index()).size()==1) {
	 ++r; ++skipped; continue;
      }
      skipped=0;
      g=gcd_of_sequence(e);
      do {
	 if (abs_equal(g,*e)) {
	    can_reduce=false;
	    if (g==1 || g==gcd_of_sequence(entire(M.col(e.index())))) {
	       can_eliminate=true;
	       c=e.index();
	       pivot_elem=*e;
	       break;
	    }
	 }
	 ++e;
      } while (!e.at_end());

      if (can_reduce) {
	 e=entire(M.row(r));
	 pivot_elem=*e;
	 c=e.index();
	 for (++e; !e.at_end(); ++e) {
	    gcd_ext(pivot_elem,*e,g2,U.a_ij,U.a_jj);
	    if (U.a_jj) {
	       if (U.a_ij) {
		  U.i=c; U.j=e.index();
		  U.a_ii=-div_exact(*e,g2);  U.a_ji=div_exact(pivot_elem,g2);
		  M.multiply_from_right(U);
		  Logger.from_right(U);
	       }
	       pivot_elem=*e;
	       c=e.index();
	       if (g==g2) break;
	    }
	 }
	 can_eliminate= g==1 || g==gcd_of_sequence(entire(M.col(c)));

#ifdef POLY_DEBUG
	 stat.gather(M);
	 if (do_dump) {
	    cout << "Smith_normal_form(reducing " << line_name << " " << r << " to " << g << " ):\n";
	    dump_matrix(cout,M);
	 }
#endif
      }

      if (can_eliminate) {
	 if (!Logger_dummy) {
	    U.i=r; U.a_ii=1; U.a_jj=1; U.a_ij=0;
	 }
	 int next_r=-1;
	 for (typename Entire<typename _Matrix::col_type>::iterator e=entire(M.col(c)); !e.at_end(); ) {
	    if ((U.j=e.index())==r) { ++e; continue; }
	    if (next_r<0) next_r=U.j;
	    if (abs_equal(pivot_elem,*e)) {
	       if (pm::sign(pivot_elem)==pm::sign(*e)) {
		  if (!Logger_dummy) U.a_ji=-1;
		  ++e;
		  M.row(U.j) -= M.row(r);
	       } else {
		  if (!Logger_dummy) U.a_ji=1;
		  ++e;
		  M.row(U.j) += M.row(r);
	       }
	    } else {
	       U.a_ji=-div_exact(*e,pivot_elem);
	       ++e;
	       M.row(U.j) += M.row(r) * U.a_ji;
	    }
	    Logger.from_left(U);
	 }
	 if (!Logger_dummy) {
	    U.i=c; U.a_ji=0;
	 }
	 for (typename Entire<typename _Matrix::row_type>::iterator e=entire(M.row(r)); !e.at_end(); ) {
	    if ((U.j=e.index())==c) { ++e; continue; }
	    if (!Logger_dummy) {
	       U.a_ij=-div_exact(*e,pivot_elem);
	       Logger.from_right(U);
	    }
	    M.row(r).erase(e++);
	 }

#ifdef POLY_DEBUG
	 stat.gather(M);
	 if (do_dump) {
	    cout << "Smith_normal_form(eliminating " << line_name << " " << r << " [" << c << "] ):\n";
	    dump_matrix(cout,M);
	 }
#endif
	 if (next_r>=0) r=next_r;
	 else ++r;
      } else {
	 ++r;
      }
   }
   return skipped;
}

template <typename E, typename CompanionLogger>
int Smith_normal_form(SparseMatrix<E>& M, std::list< std::pair<E,int> >& torsion, CompanionLogger& Logger)
{
   const bool Logger_dummy=pm::derived_from_instance<CompanionLogger, dummy_companion_logger>::value;

#ifdef POLY_DEBUG
   const bool do_dump=getenv("SMITH_DUMP");
   if (do_dump) {
      cout << "Smith_normal_form(initial):\n";
      dump_matrix(cout,M);
   }
   pm::SparseMatrixStatistics<E> stat;
   stat.gather(M);
   cout << "Smith_normal_form(initial statistics):\n" << stat;

   while (Smith_normal_form_steps(M, Logger , stat, do_dump, "row") < M.rows() &&
	  Smith_normal_form_steps(T(M), transpose_logger(Logger) , stat, do_dump, "col") < M.cols()) ;

   cout << "Smith_normal_form(final statistics):\n" << stat << endl;
#else
   while (Smith_normal_form_steps(M, Logger) < M.rows() &&
	  Smith_normal_form_steps(T(M), transpose_logger(Logger)) < M.cols()) ;
#endif
   int rank=0;
   torsion.clear();
   for (typename Entire< Rows< SparseMatrix<E> > >::iterator r=entire(rows(M)); !r.at_end(); ++r) {
      if (!r->empty()) {
	 typename SparseMatrix<E>::row_type::iterator e=r->begin();
	 ++rank;
	 if (!abs_equal(*e,1))
	    torsion.push_back(std::make_pair(abs(*e), e.index()));
      }
   }

   E g, p, q, k, l;
   for (typename std::list< std::pair<E,int> >::iterator t=torsion.begin(), t_end=torsion.end();
	t != t_end; ++t) {
      typename std::list< std::pair<E,int> >::iterator t2=t;
      for (++t2; t2 != t_end; ) {
	 if (Logger_dummy) {
	    g=gcd(t->first, t2->first);
	 } else {
	    gcd_ext(t->first, t2->first, g, p, q);
	 }
	 if (t->first == g) {
	    std::swap(*t, *t2);  ++t2;
	 } else if (t2->first == g) {
	    ++t2;
	 } else {
	    l=div_exact(t2->first,g);
	    if (!Logger_dummy) {
	       k=-div_exact(t->first,g);
	       Logger.from_left(SparseMatrix2x2<E>(M.col(t->second).begin().index(), M.col(t2->second).begin().index(),
						   l, p*k,
						   1,   q) );
	       Logger.from_right(SparseMatrix2x2<E>(t->second, t2->second,
						    q*l, p,
						    k,   1) );
	    }
	    t->first *= l;
	    if (g==1) {
	       t2=torsion.erase(t2);
	    } else {
	       t2->first=g;  ++t2;
	    }
	 }
      }
   }
   return rank;
}

template <typename E, typename CompanionLogger>
int eliminate_ones(SparseMatrix<E>& M, Bitset& elim_rows, Bitset& elim_cols, CompanionLogger& Logger)
{
   const bool Logger_dummy=pm::derived_from_instance<CompanionLogger, dummy_companion_logger>::value;

#ifdef POLY_DEBUG
   const bool do_dump=getenv("SMITH_DUMP");
   if (do_dump) {
      cout << "eliminate_ones(initial):\n";
      dump_matrix(cout,M);
   }
   pm::SparseMatrixStatistics<E> stat;
   stat.gather(M);
   cout << "eliminate_ones(initial statistics):\n" << stat;
#endif

   SparseMatrix2x2<E> U;
   const int R=M.rows(), C=M.cols();
   elim_rows.clear(); elim_rows.reserve(R);
   elim_cols.clear(); elim_cols.reserve(C);
   int count=0;

   for (int c=0, c_last=C;  c != c_last;  c==C && c_last<C && (c=0)) {    // repeat until matrix is empty or contains no +1/-1 more
      if (M.col(c).empty() || elim_cols.contains(c)) { ++c; continue; }

      typename Entire<typename SparseMatrix<E>::col_type>::iterator e=entire(M.col(c));
      while (!abs_equal(*e,1) && !(++e).at_end()) ;
      if (e.at_end()) { ++c; continue; }

      int r=e.index();

      if (!Logger_dummy) {
	 U.i=r; U.a_ii=1; U.a_jj=1; U.a_ij=0;
      }
      for (typename Entire<typename SparseMatrix<E>::col_type>::iterator e2=entire(M.col(c)); !e2.at_end(); ) {
	 if ((U.j=e2.index())==r) { ++e2; continue; }
	 if (abs_equal(*e,*e2)) {
	    if (pm::sign(*e)==pm::sign(*e2)) {
	       if (!Logger_dummy) U.a_ji=-1;
	       ++e2;
	       M.row(U.j) -= M.row(r);
	    } else {
	       if (!Logger_dummy) U.a_ji=1;
	       ++e2;
	       M.row(U.j) += M.row(r);
	    }
	 } else {
	    U.a_ji= (*e)>0 ? -(*e2) : (*e2);
	    ++e2;
	    M.row(U.j) += U.a_ji * M.row(r);
	 }
	 Logger.from_left(U);
      }
      if (Logger_dummy) {
	 M.row(r).clear();
	 ++count;
      } else {
	 U.i=c; U.a_ji=0;
	 for (typename Entire<typename SparseMatrix<E>::row_type>::iterator e2=entire(M.row(r)); !e2.at_end(); ) {
	    if ((U.j=e2.index())==c) { ++e2; continue; }
	    U.a_ij= (*e)>0 ? -(*e2) : (*e2);
	    Logger.from_right(U);
	    M.row(r).erase(e2++);
	 }
      }
      elim_rows+=r;  elim_cols+=c;
#ifdef POLY_DEBUG
      stat.gather(M);
      if (do_dump) {
	 cout << "eliminate_ones(eliminating col " << c << " [" << r << "] ):\n";
	 dump_matrix(cout,M);
      }
#endif
      c_last=c++;
   }

#ifdef POLY_DEBUG
   cout << "eliminate_ones(final statistics):\n" << stat << endl;
#endif
   return count;
}

} }

#endif // _POLYMAKE_SMITH_NORMAL_FORM_H

// Local Variables:
// mode:C++
// c-basic-offset:3
// End:
