/* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -I../include -g -c -Wall TmpFGLM.C" -*- */
  //   Copyright (c)  2006  Stefan Kaspar
  
  //   This file is part of the source of CoCoALib, the CoCoA Library.
  
  //   CoCoALib is free software; you can redistribute it and/or modify
  //   it under the terms of the GNU General Public License (version 3)
  //   as published by the Free Software Foundation.  A copy of the full
  //   licence may be found in the file COPYING in this directory.
  
  //   CoCoALib 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.
  
  //   You should have received a copy of the GNU General Public License
  //   along with CoCoA; if not, see <http://www.gnu.org/licenses/>.
  
  #ifdef HAVE_CONFIG_H
  #include "config.h"
  #endif
  #include "first.h"
  #ifdef USE_GMP_REPLACEMENTS
  #undef HAVE_LIBCOCOA
  #endif
  #ifdef HAVE_LIBCOCOA
  #include "TmpFGLM.H"
  #include "TmpLESystemSolver.H"
  #include "CoCoA/DenseMatrix.H"
  #include "CoCoA/symbol.H"
  #include "CoCoA/QBGenerator.H"
  #include "CoCoA/RingDistrMPolyInlPP.H"
  #include "CoCoA/RingHom.H"
  #include "CoCoA/SparsePolyRing.H"
  
  // #include <cstddef>
  using std::size_t;
  #include <set>
  using std::set;
  // #include <list> // Included by QBGenerator.H
  using std::list;
  #include <map>
  using std::map;
  // #include <memory> // Included by SparsePolyRing.H
  using std::auto_ptr;
  #include <utility>
  using std::make_pair;
  // #include <vector>
  using std::vector;
  
  namespace CoCoADortmund
  {
    using namespace CoCoA;
  
    // Used for constructing matrices to solve linear equation systems
    // (Kind of sparse matrix representation)
    struct MatrixMapEntry
    {
      MatrixMapEntry(ConstRefRingElem r): rhs(r) {}
      std::vector<std::size_t> VarIndices; // Column numbers in linear equation matrix
      std::vector<RingElem> coeffs; // Matrix components in linear equation matrix
      RingElem rhs; // Right hand side component in linear equation
    };
  
  
    // Embed element p of PPM into another PPM with different term ordering
    PPMonoidElem PPIntoOtherPPM(ConstRefPPMonoidElem p, const PPMonoid& OtherPPM)
    {
      vector<long> expv;
      exponents(expv, p);
      return PPMonoidElem(OtherPPM, expv);
    }
  
    // Update matrix map during main loop
    // !! Weakly exception safe, writable parameters might get rendered useless !!
    void UpdateMatrixMap(map<PPMonoidElem, struct MatrixMapEntry>& MatrixMap, size_t& NumCols, QBGenerator& NewQB, ConstRefRingElem p, ConstRefPPMonoidElem t)
    {
      for (SparsePolyIter m = BeginIter(p); !IsEnded(m); ++m)
      {
        map<PPMonoidElem, struct MatrixMapEntry>::iterator entry = MatrixMap.find(PP(m));
        // ToDo: Avoid this tedious if-else block
        if (entry == MatrixMap.end())
        {
          ring K = CoeffRing(AsSparsePolyRing(owner(p)));
          struct MatrixMapEntry NewEntry(zero(K));
          NewEntry.VarIndices.push_back(NumCols);
          NewEntry.coeffs.push_back(coeff(m));
          MatrixMap.insert(make_pair(PP(m), NewEntry));
        }
        else
        {
          entry->second.VarIndices.push_back(NumCols);
          entry->second.coeffs.push_back(coeff(m));
        }
      }
      ++NumCols;
  
      // Update quotient basis NewQB
      NewQB.myCornerPPIntoQB(t);
    }
  
    // FGLM implementation
    // exception safe
    void FGLMBasisConversion(vector<RingElem>& NewGB, const vector<RingElem>& OldGB, const PPOrdering& NewOrdering)
    {
      if (OldGB.empty())
        CoCoA_ERROR(ERR::nonstandard, "FGLMBasisConversion: empty Groebner Basis vector");
  
      // Check if generated ideal is zero-dimensional
      const ideal I(AsSparsePolyRing(owner(OldGB.front())), OldGB);
      if (!IsZeroDim(I))
        CoCoA_ERROR(ERR::nonstandard, "FGLMBasisConversion: ideal must be 0-dimensional");
  
      // Initialization of objects needed for computation
      const SparsePolyRing Kx = AsSparsePolyRing(owner(OldGB.front()));
      const ring K = CoeffRing(Kx);
      const PPMonoid PPMon = PPM(Kx);
      const RingElem FieldOne(one(K));
  
      // Adjust K[x_1, ..., x_n] and PPM(K[x_1, ..., x_n]) to correct term ordering
      const PPMonoid PPMAdjusted = NewPPMonoid(symbols(Kx), NewOrdering);
      const SparsePolyRing KxAdjusted = NewPolyRing(CoeffRing(Kx), PPMAdjusted);
  
  //    // Identity mapping Kx -> KxAdjusted
  //    // (Not yet used)
  //    RingHom KxToKxAdjusted = PolyRingHom(Kx, KxAdjusted, CoeffEmbeddingHom(KxAdjusted), indets(KxAdjusted));
  
      // These will hold the (temporary) basis elements
      vector<RingElem> NewGBTmp; // Holds elements in K[x_1, ..., x_n] w.r.t. new term ordering
      QBGenerator NewQB(PPMAdjusted); // Holds the new QB in PPM with new term ordering
      map<PPMonoidElem, struct MatrixMapEntry> MatrixMap; // Used to create the matrices for the linear dependency check
      size_t NumCols = 1;
  
      // FGLM algorithm start
      PPMonoidElem t(PPMAdjusted);
      t = one(PPMAdjusted);
      NewQB.myCornerPPIntoQB(t);
      struct MatrixMapEntry entry(zero(K));
      entry.VarIndices.push_back(0);
      entry.coeffs.push_back(FieldOne);
      MatrixMap.insert(make_pair(PPIntoOtherPPM(t, PPMon), entry)); // Want to keep the keys (PPs) in PPM with old term ordering
  
      while (!NewQB.myCorners().empty())
      {	
        t = NewQB.myCorners().front(); // t is element of PPMAdjusted
        RingElem h(Kx);
        h = NR(monomial(Kx, FieldOne, PPIntoOtherPPM(t, PPMon)), OldGB);
        vector< map<PPMonoidElem, struct MatrixMapEntry>::iterator > ResetPos;
        // Prepare creation of linear equation system and check if a solution can
        // exist (simple check if equations like "0 = c" would occur with c not
        // equal to 0)
        bool RemainderMightBeIndependent = true;
        for (SparsePolyIter m = BeginIter(h); !IsEnded(m); ++m)
        {
          map<PPMonoidElem, struct MatrixMapEntry>::iterator entry = MatrixMap.find(PP(m));
          // If h contains a term that is not already in the matrix map
          // then h is linearly independent of the previously computed remainders
          if (entry == MatrixMap.end())
          {
            RemainderMightBeIndependent = false; // Only reset-operations after matrix creation code below need to be carried out
            UpdateMatrixMap(MatrixMap, NumCols, NewQB, h, t);
            break;
          }
          else
          {
            // Right hand side component equals coeff(m)
            entry->second.rhs = -coeff(m);
            ResetPos.push_back(entry);
          }
        }
  
        // If it is not yet clear if the current remainder is linearly dependent on
        // the previously computed remainders we have to create a linear equation
        // system and try to solve it
        if (RemainderMightBeIndependent)
        {
          // Create a system of linear equations from matrix map
          matrix M = NewDenseMat(K, MatrixMap.size(), NumCols),
                 b = NewDenseMat(K, MatrixMap.size(), 1);
          size_t j = 0; // Current row
  
          for (map<PPMonoidElem, struct MatrixMapEntry>::iterator entry = MatrixMap.begin(); entry != MatrixMap.end(); ++entry, ++j)
          {
            vector<size_t>& VarIndices = entry->second.VarIndices;
            vector<RingElem>& coeffs = entry->second.coeffs;
  
            // Set matrix components
            for (size_t k = 0; k < VarIndices.size(); ++k)
              SetEntry(M, j, VarIndices[k], coeffs[k]);
  
            // Set right hand side component
            SetEntry(b, j, 0, entry->second.rhs);
          }
  
          // Check if M*x = b has a solution
          matrix x = NewDenseMat(K, NumCols, 1);
          if (LESystemSolver(x, M, b))
          {
            // Compute new Groebner Basis polynomial
            RingElem g(KxAdjusted); // g in K[x_1, ..., x_n] with correct term ordering
            g = monomial(KxAdjusted, FieldOne, t);
            const vector<PPMonoidElem>& QB = NewQB.myQB();
            for (size_t i = 0; i < NumCols; ++i)
              g += monomial(KxAdjusted, x(i, 0), QB[i]);
            NewGBTmp.push_back(g);
  
            // Update quotient basis NewQB
            NewQB.myCornerPPIntoAvoidSet(t);
          }
          else
            UpdateMatrixMap(MatrixMap, NumCols, NewQB, h, t);
        }
  
        // Reset right hand side components in matrix map
        for (size_t i = 0; i < ResetPos.size(); ++i)
          ResetPos[i]->second.rhs = zero(K);
      }
  
      // Swap computed Groebner Basis
      swap(NewGB, NewGBTmp);
    }
  
  } // End of namespace CoCoA
  
  #endif