#include <map>

class MlfbAlgorithm {
public:
  MlfbAlgorithm(const mpz_class& initialNumber,
		const vector<vector<Exponent> >& terms,
		const Degree** degrees,
		/*TranspositionTable& transpositionTable,*/
		bool printProgress = false):
    _dim(terms[0].size()),
    _maximumDegreeSeen(0),
    _printProgress(printProgress),
    _degrees(degrees),
    _handler(_dim, terms),
    _solution(_dim),
    _callCounts(_dim),
    _lastProgressPos(0) {

    ExternalTerm b(_handler);
    SortedTermList termList(_handler);

    recurse(b, -initialNumber, termList);

    //*
    for (int i = 0; i < _dim; ++i)
      cout << "Level " << i + 1 << " had "
	   << _callCounts[i] << " calls. " << endl;
    //*/
  }

  Degree getCalculatedFrobeniusNumber() {
    return _maximumDegreeSeen;
  }

private:
  vector<ExternalTerm> bs;

  void recurse(const ExternalTerm& b, const Degree& degree, const SortedTermList& terms) {
    ASSERT(!terms.empty());

    int position = terms.position();

    ++_callCounts[position];

    if (position == _dim - 2) {
      baseCase(b, degree, terms);
      return;
    }

    //cout << " at " << position << "with b="<<b<<endl << terms;

    if (canSkipDueToUpperBound(terms, degree))
      return;

    ExternalTerm newB(b);

    // Set up [it, terms), which is the range we will search through.
    SortedTermList::iterator termsEnd = terms.end();
    SortedTermList::iterator it = terms.begin();
    terms.getFirstLarger(b[position], it);

    /*
    cout << "starting from ";
      _handler.write(&*it, cout);
      cout << "with b=" << b << endl;//*/

    // Set up the list of terms for the level of recursion lower than
    // this. We will incrementally update this list.
    SortedTermList nextTerms(b, _handler, position + 1);
    nextTerms.clearAndSetParent(terms);
    SortedTermList::iterator toAddEnd =
      terms.copyTermsWithLeadingZeroInto(nextTerms);    

    Degree newDegree;
    for (; it != termsEnd; ++it) {
      if (_printProgress && terms.position() <= 1)
	printProgress(terms, it);
      
      // Does not change entries prior to position, so we can keep
      //using newB in later iterations without resetting it to b.
      _handler.lcm(newB, b, &*it, position);
      newB.setExponent(position, _handler.getExponent(&*it, position) - 1);

      newDegree = degree + _degrees[position][newB[position]];

      _solution[position] = &*it;

      SortedTermList::iterator toAddBegin = toAddEnd;
      terms.getFirstLarger(newB[position], toAddEnd);
      nextTerms.insert(toAddBegin, toAddEnd);

      if (!(newB == b)) {
	SortedTermList nextTermsOpt(nextTerms.begin().getPointer(),
				    nextTerms.end().getPointer(),
				    nextTerms, newB);
	recurse(newB, newDegree, nextTermsOpt);
      } else {
	recurse(newB, newDegree, nextTerms);
      }
    }
  }

  void baseCase(const ExternalTerm& b, const Degree& degree, const SortedTermList& terms) {
#ifdef PROFILE
    if ((int)terms.size() == -1) {
      baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);
      baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);
      baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);
      baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);
      baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);baseCase(b,degree,terms);
    }
#endif

    //cout << "at basecase." << terms << endl;

    ASSERT(!terms.empty());

    int position = terms.position();

    SortedTermList::iterator termsEnd = terms.end();
    SortedTermList::iterator term = terms.begin();

    SortedTermList::iterator nextTerm = term;
    ++nextTerm;

    ExternalTerm newB(b);
    Degree newDegree;
    while (nextTerm != termsEnd) {
      ASSERT(_handler.getExponent(&*nextTerm, position) > 0);
      ASSERT(_handler.getExponent(&*term, position + 1) > 0);

      newB.setExponent(position,
		       _handler.getExponent(&*nextTerm, position) - 1);
      newB.setExponent(position + 1,
		       _handler.getExponent(&*term, position + 1) - 1);
      _solution[position] = &*nextTerm;
      _solution[position + 1] = &*term;

      newDegree = degree +
	_degrees[position][newB[position]] +
	_degrees[position + 1][newB[position + 1]];

      registerSolution(newB, newDegree);
	
      ++term;
      ++nextTerm;
    }
  }

  void improveB(ExternalTerm& b, const SortedTermList& terms) {
    int position = terms.position() + 1;

    ExternalTerm lcm(b);
    SortedTermList::iterator end = terms.end();
    for (SortedTermList::iterator it = terms.begin(); it != end; ++it)
      _handler.lcm(lcm, lcm, &*it, position);

    vector<ExternalTerm> gcds;
    gcds.reserve(_dim);
    for (int i = 0; i < _dim; ++i)
      gcds.push_back(lcm);

    for (SortedTermList::iterator it = terms.begin(); it != end; ++it) {
      for (int i = position; i < _dim; ++i) {
	Exponent exponent = _handler.getExponent(&*it, i);
	if (exponent > b[i])
	  _handler.gcd(gcds[i], gcds[i], &*it, position);
      }
    }

    for (int i = position; i < _dim; ++i) {
      gcds[i][i] -= 1;
      _handler.lcm(b, b, gcds[i].getTerm(), position);
    }
  }

  bool canSkipDueToUpperBound(const SortedTermList& terms, const Degree& degree) {

    int position = terms.position();

    if (position > _dim - 3)
      return false;

    ExternalTerm lcm(_handler);
    SortedTermList::iterator end = terms.end();
    for (SortedTermList::iterator it = terms.begin(); it != end; ++it)
      _handler.lcm(lcm, lcm, &*it, position);

    Degree upperBound = degree;
    for (int i = position; i < _dim; ++i) {
      ASSERT(lcm[i] > 0);
      upperBound += _degrees[i][lcm[i] - 1];
    }

    if (upperBound <= _maximumDegreeSeen)
      return true;


    if (false) { // too expensive to compute due to GMP being called
      Degree maxMinLoss = 0; // loss of degree
      for (SortedTermList::iterator it = terms.begin(); it != end; ++it) {
	Degree minLoss = -1;
	for (int i = position; i < _dim; ++i) {
	  Degree loss = _degrees[i][lcm[i] - 1] - _degrees[i][_handler.getExponent(&*it, i)];
	  if (minLoss == -1 || loss < minLoss)
	    minLoss = loss;
	}
	if (minLoss > maxMinLoss)
	  maxMinLoss = minLoss;
      }

      if (upperBound - maxMinLoss <= _maximumDegreeSeen)
	return true;
    }

    // This works almost as well while being much faster to compute.
    // It is still too expensive a bound.
    if (false) {
      SortedTermList::iterator bestCandidate = end;
      int maxMinValue = -1;
      for (SortedTermList::iterator it = terms.begin(); it != end; ++it) {
	int minValue = -1;
	for (int i = position; i < _dim; ++i) {
	  int value = lcm[i] - _handler.getExponent(&*it, i);
	  if (minValue == -1 || value < minValue) {
	    minValue = value;
	  }
	}
	if (maxMinValue == -1 || minValue > maxMinValue) {
	  maxMinValue = minValue;
	  bestCandidate = it;
	  
	}
      }
      ASSERT(bestCandidate != end);
	
      Degree minLoss = -1;
      for (int i = position; i < _dim; ++i) {
	Degree loss = _degrees[i][lcm[i] - 1] -
	  _degrees[i][_handler.getExponent(&*bestCandidate, i)];
	if (minLoss == -1 || loss < minLoss)
	  minLoss = loss;
      }
      

      if (upperBound - minLoss <= _maximumDegreeSeen)
	return true;
    }


    return false;
  }

  void registerSolution(const ExternalTerm& b, const Degree& degree) {
    //cout << "solution: " << b << endl;

    ++_callCounts[_dim - 1];

    //static map<ExternalTerm, vector<const Term*> > seen;
	  
    	
    //static int co = 0, dupCo = 0;

    //    if (!seen[b].empty()) {
    /* 
       cout << "!!!!!!!!!!!!" << endl;
       cout << "solution " << b;
       cout << " attained by two different label vectors." << endl;

       cout << "vector1:" << endl;
       for (unsigned int i = 0; i < seen[b].size(); ++i) {
       cout << i + 1 << ": ";
       _handler.write(seen[b][i], cout);
       cout << endl;
       }

       cout << "vector2:" << endl;
       for (unsigned int i = 0; i < _solution.size(); ++i) {
       cout << i + 1 << ": ";
       _handler.write(_solution[i], cout);
       cout << endl;
       }

       exit(0);
    */	
    //      ++dupCo;
    //}


    //    seen[b] = _solution;

    //++co;

    //cout << co << ' ' << dupCo << ' '<< ((double)dupCo)/ co << '\n';
    
    if (degree > _maximumDegreeSeen)
      _maximumDegreeSeen = degree;
  }

  void printProgress(const SortedTermList& terms,
		     SortedTermList::iterator it) {
    // Do not print progress more than once every 5 seconds.
    if (terms.position() >= _lastProgressPos && _progressTimer.getSeconds() < 5)
      return;
    _progressTimer.reset();
    _lastProgressPos = terms.position();

    // begin is the place we actually start the work from.
    SortedTermList::iterator begin = terms.begin();
    terms.getFirstLarger(0, begin);

    int doneCount = distance(begin, it);
    int all = terms.size() - distance(terms.begin(), begin);
    double doneRatio = ((double)doneCount)/all;
    
    cout << "At level=" << terms.position() << ": "
	 << setprecision(3)
	 << doneCount << '/' << all << '='
	 << doneRatio * 100.0 << "% done in " << _timer << '.' << endl;
  }

  int _dim;
  Degree _maximumDegreeSeen;
  bool _printProgress;
  const Degree** _degrees;
  Timer _timer;
  Timer _progressTimer;
  TermHandler _handler;
  vector<const Term*> _solution;
  //  TranspositionTable& _transpositionTable;
  vector<int> _callCounts;
  int _lastProgressPos;
};

// returns the maximal id in this position.
inline Exponent compressRange(const vector<vector<mpz_class> >& input,
			      vector<vector<Exponent> >& output,
			      const mpz_class& degree,
			      Degree*& precomputedDegrees,
			      int position) {
  // Collect the exponents.
  vector<mpz_class> exponents;
  for (unsigned int i = 0; i < input.size(); ++i)
    exponents.push_back(input[i][position]);

  // Sort the exponents and remove duplicates.
  sort(exponents.begin(), exponents.end());
  vector<mpz_class>::iterator uniqueEnd =
    unique(exponents.begin(), exponents.end());
  exponents.erase(uniqueEnd, exponents.end());

  // Construct a map from the exponent values that we will need to
  // look at to the id numbers we will be working on from now on.
  //
  // The set of id numbers must be such that every exponent that
  // appears has an id. Also, each exponent minus one must have an id,
  // as the partial solution will use such values. The ids must have
  // the same sorted order as the numbers they are ids for.
  //
  // This of course makes it impossible to see how much an entry
  // should contribute towards the degree just from looking at the id,
  // so we must precompute the degree that corresponds to each id.
  map<mpz_class, int> rangeReduction;
  rangeReduction[0] = 0;

  precomputedDegrees = new Degree[exponents.size() * 2 + 1];
  precomputedDegrees[0] = 0;

  // range is the largest id we have assigned so far.
  unsigned int range = 0;
  mpz_class lastExponent = 0;

  ASSERT(exponents[0] == 0);
  for (unsigned int i = 1; i < exponents.size(); ++i) {
    const mpz_class& exponent = exponents[i];

    // We must only assign exponent - 1 a new id if 
    // we have not already done so.
    if (exponent > lastExponent + 1) {
      ++range;
      rangeReduction[exponent - 1] = range;
      precomputedDegrees[range] = (exponent - 1) * degree;
    }

    ++range;
    rangeReduction[exponent] = range;
    precomputedDegrees[range] = exponent * degree;

    lastExponent = exponent;
  }

  // Apply the range compression map.
  for (unsigned int i = 0; i < input.size(); ++i)
    output[i][position] = rangeReduction[input[i][position]];

  return range;
}

inline mpz_class computeFrobeniusNumber(const vector<vector<mpz_class> >& terms,
					const vector<mpz_class>& degrees,
					bool printProgress = false) {
  ASSERT(!terms.empty());

  if (degrees.size() == 2)
    return degrees[0] * degrees[1] - degrees[0] - degrees[1];

  // Before we can use MlfbAlgorithm, we need to compress the range of
  // the inputs so that we are sure the exponents will fit in 32 bits.

  unsigned int dimension = terms[0].size();

  vector<vector<Exponent> > rangeReducedTerms(terms.size());
  for (unsigned int i = 0; i < terms.size(); ++i)
    rangeReducedTerms[i].resize(dimension);

  Degree** precomputedDegrees = new Degree*[dimension];

  //  TranspositionTable transpositionTable(dimension);

  for (unsigned int position = 0; position < dimension; ++position) {
    /*Exponent range = */compressRange(terms,
				       rangeReducedTerms,
				       degrees[position + 1],
				       precomputedDegrees[position],
				       position);
    //    transpositionTable.setRange(position, range);
  }

  // Now we can use MlfbAlgorithm with the compressed range.
  MlfbAlgorithm algo(degrees[0],
		     rangeReducedTerms,
		     (const Degree**)precomputedDegrees,
		     //		     transpositionTable,
		     printProgress);

  // Clean up the allocated memory.
  for (unsigned int i = 0; i < dimension; ++i)
    delete[] precomputedDegrees[i];
  delete[] precomputedDegrees;

  return algo.getCalculatedFrobeniusNumber();
}