/*
 * This file is part of PARAM.
 *
 * PARAM 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 3 of the License, or
 * (at your option) any later version.
 *
 * PARAM 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 PARAM.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Copyright 2009 Ernst Moritz Hahn (emh@cs.uni-sb.de)
 */

#include <assert.h>
#include <string>
#include <sstream>
#include "CoCoA/library.H"
#include "CoCoA/PPOrdering.H"
#include "CancellatorCoCoA.h"
#include "Base.h"
#include "Polynomial.h"

using namespace CoCoA;
using namespace std;

GlobalManager manager;

namespace Rational {
  std::vector<std::string> CancellatorCoCoA::symbolStrings;
  CoCoA::SparsePolyRing *CancellatorCoCoA::polyRing(NULL);
  CoCoA::SparsePolyRing *CancellatorCoCoA::polyRingGCD(NULL);
  CoCoA::RingHom *CancellatorCoCoA::forth(NULL);
  CoCoA::RingHom *CancellatorCoCoA::back(NULL);
  int CancellatorCoCoA::cleanupInt(atexit(CancellatorCoCoA::clear));
  
  string CancellatorCoCoA::int2string(unsigned number) {
    stringstream numberStringStream;
    numberStringStream << number;
    string numberString(numberStringStream.str());
    for (unsigned index(0); index < numberString.length(); index++) {
      numberString[index] -= '0';
      numberString[index] += 'A';
    }

    return numberString;
  }


  void CancellatorCoCoA::clear() {
    if (NULL != forth) {
      delete forth;
      forth = NULL;
    }
    if (NULL != back) {
      delete back;
      back = NULL;
    }
    if (NULL != polyRing) {
      delete polyRing;
      polyRing = NULL;
    }
    if (NULL != polyRingGCD) {
      delete polyRingGCD;
      polyRingGCD = NULL;
    }
    symbolStrings.clear();
  }

  void CancellatorCoCoA::start() {
    if (0 == symbolStrings.size()) {
      return;
    }
    const unsigned numSymbols(symbolStrings.size());
    ring Q = RingQ();
    vector<CoCoA::symbol> params;
    for (unsigned symNr(0); symNr < numSymbols; symNr++) {
      params.push_back(symbol(int2string(symNr)));
    }

    PPOrdering order = NewLexOrdering(numSymbols);
    polyRing = new SparsePolyRing(NewPolyRing(Q, params, order));
    polyRingGCD = new SparsePolyRing(NewPolyRing(Q, params));

    vector<RingElem> indetImages;
    for (unsigned index(0); index < numSymbols; index++) {
      vector<long> expv(numSymbols);
      expv[index] = 1;
      RingElem c(RingQ(), 1);
      indetImages.push_back(monomial(*polyRingGCD, c, expv));
    }
    forth = new RingHom(PolyAlgebraHom(*polyRing, *polyRingGCD, indetImages));

    indetImages.clear();
    for (unsigned index(0); index < numSymbols; index++) {
      vector<long> expv(numSymbols);
      expv[index] = 1;
      RingElem c(RingQ(), 1);
      indetImages.push_back(monomial(*polyRing, c, expv));
    }
    back = new RingHom(PolyAlgebraHom(*polyRingGCD, *polyRing, indetImages));
  }

  void CancellatorCoCoA::convert(const Polynomial *poly, RingElem &result) {
    const unsigned *monomials(poly->getMonomials());
    const mpz_t *coefficients(poly->getCoefficients());
    const unsigned numSymbols(Base::getNumSymbols());
    const unsigned numTerms(poly->getNumTerms());
    ring Q = RingQ();
    vector<long> expv(numSymbols);
    RefRingElem refResult(result);
    for (unsigned termNr(0); termNr < numTerms; termNr++) {
      for (unsigned symNr(0); symNr < numSymbols; symNr++) {
        expv[symNr] = monomials[termNr * numSymbols + symNr];
      }
      ZZ z(coefficients[termNr]);
      RingElem c(Q, z);
#if 1
      PushBack(refResult, c, expv);
#else
      result += monomial(*polyRing, c, expv);
#endif
    }
  }

  void CancellatorCoCoA::convert(const RingElem &poly, Polynomial *result) {
    const unsigned numSymbols(Base::getNumSymbols());
    delete[] result->monomials;
    for (unsigned termNr(0); termNr < result->numTerms; termNr++) {
      mpz_clear(result->coefficients[termNr]);
    }
    delete[] result->coefficients;

    result->numTerms = NumTerms(poly);
    result->monomials = new unsigned[result->numTerms * numSymbols];
    result->coefficients = new mpz_t[result->numTerms];

    unsigned termNr(0);
    for (SparsePolyIter it(BeginIter(poly));
         it != EndIter(poly); it++) {
      const RingElem &cf(coeff(it));
      const RingElem &nm(num(cf));
      const PPMonoidElem &pp = PP(it);
      vector<long> expv(numSymbols);
      exponents(expv, pp);
      for (unsigned symNr(0); symNr < numSymbols; symNr++) {
        result->monomials[termNr * numSymbols + symNr] = expv[symNr];
      }
      
      mpz_init(result->coefficients[termNr]);
      ZZ z = floor(nm);
      mpz_set(result->coefficients[termNr], mpzref(z));
      termNr++;
    }
  }

  void CancellatorCoCoA::cancel(RingElem &poly1, RingElem &poly2) {
    RingHom phi = CanonicalHom(RingZ(), *polyRing);

    const RingElem poly1GCD((*forth)(poly1));
    const RingElem poly2GCD((*forth)(poly2));

    const RingElem divisorGCD(gcd(poly1GCD, poly2GCD));
    const RingElem divisor((*back)(divisorGCD));

    RingElem content(RingZ(), 0);
    for (SparsePolyIter it(BeginIter(divisor));
         it != EndIter(divisor); it++) {
      const RingElem &cf(coeff(it));
      const RingElem &nm(num(cf));
      content = gcd(content, nm);
    }
    divisor /= phi(content);
    poly1 /= divisor;
    poly2 /= divisor;
  }

  void CancellatorCoCoA::cancel(Polynomial *poly1, Polynomial *poly2) {
    if (0 == symbolStrings.size()) {
      return;
    }
    if (poly2->isConstant()) {
      return;
    }
    if (*poly1 == *poly2) {
      poly1->setToConstant(1);
      poly2->setToConstant(1);
    } else {
      RingElem cocoa1(*polyRing);
      RingElem cocoa2(*polyRing);
      convert(poly1, cocoa1);
      convert(poly2, cocoa2);
      cancel(cocoa1, cocoa2);
      convert(cocoa1, poly1);
      convert(cocoa2, poly2);
    }
  }

}
