f4_reduction.hpp

#ifndef __F4_REDUCTION_HPP__
#define __F4_REDUCTION_HPP__

/*****************************************************************************\
* This file is part of DynGB.                                                 *
*                                                                             *
* DynGB 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 of the License, or           *
* (at your option) any later version.                                         *
*                                                                             *
* DynGB 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 DynGB. If not, see <http://www.gnu.org/licenses/>.               *
\*****************************************************************************/

#include "system_constants.hpp"

#include "fields.hpp"
#include "monomial.hpp"
#include "monomial_ordering.hpp"
#include "polynomial.hpp"
#include "critical_pair.hpp"
#include "f4_hash.hpp"

#include <set>
using std::set;
#include <list>
using std::list;
#include <map>
using std::map;
#include <vector>
using std::vector;
#include <cstdlib>
#include <utility>
using std::pair;
#include <thread>
using std::thread;
#include <mutex>
using std::mutex;

/**
  @ingroup GBComputation
  @brief equivalent to @c buchberger(), but for Faug&egrave;re&rsquo;s F4 algorithm
  @param F generators of a polynomial ideal
  @return a Gr&ouml;bner basis of the ideal generated by @f$ F @f$
    with respect to the ordering already assigned to its polynomials
*/
list<Constant_Polynomial *> f4_control(const list<Abstract_Polynomial *> &F);

/**
  @brief used to compare monomials for STL containers such as @c map
*/
struct MonCmp {
  /**
    @brief returns @c true iff the monomial pointed to by @p t is smaller than
      the monomial pointed to by @p u
    @param t monomial we're comparing
    @param u monomial we're comparing
    @return @c true iff the monomial pointed to by @p t is smaller than
      the monomial pointed to by @p u
  */
  bool operator() (const Monomial * t, const Monomial * u) const {
    return *t < *u;
  }
};

/**
  @brief Implementation of Faug&egrave;re&rsquo;s F4 algorithm.
  @ingroup GBComputation
  @details Currently computes a Gr&ouml;bner basis by selecting several
    s-polynomials of lowest lcm degree. Data is stored in a semi-sparse matrix
    format, with each row a contiguous array of entries (@c A).
    Each row of @c A begins at the position indicated
    by the corresponding @c offset,
    and its first non-zero entry appears at the position indicated by
    the corresponding @c head.
    So the leading coefficient of the polynomial in row @c k appears in
    <c>A[head[k]]</c> and the leading monomial appears in
    <c>M[offset[k]+head[k]]</c>.
    Starting from <c>head[k]</c>, row @c k is actually dense.
    While this is not sparse, it does succeed in saving more space
    than one might expect.
*/
class F4_Reduction_Data {
public:
  /** @name Construction */
  ///@{
  /**
    @brief encapsulation of one step of the F4 algorithm for the polynomials
        indicated by @p P and @p B
    @param P list of critical pairs that will create a new basis; matrix will
        have this many rows
    @param B list of polynomials currently in the basis
  */
  F4_Reduction_Data(
      const list<Critical_Pair_Basic *> & P,
      const list<Abstract_Polynomial *> & B
  );
  /**
    @brief adds monomials of @f$ ug @f$ to @c M_builder
    @param g polynomial, such as an s-polynomial generator or a reducer
    @param u monomial to multiply to @p g to add monomials
    @param new_row if @c true, adds the leading monomial; otherwise, adds only
        trailing monomials
  */
  void add_monomials(
      const Abstract_Polynomial *g,
      const Monomial & u,
      bool new_row = false
  );
  /**
    @brief creates the matrix
    @details called internally by constructor, which has already performed some
        setup; inadvisable to use elsewhere
    @param P list of critical pairs that will generate the matrix
  */
  void initialize_many(const list<Critical_Pair_Basic *> & P);
  ///@}
  /** @name Destruction */
  ///@{
  /**
    @brief releases space for the matrix and deletes any strategies not already
    set to @c nullptr
  */
  ~F4_Reduction_Data();
  ///@}
  /** @name Basic properties */
  ///@{
  /**
    @brief returns @c true iff all the entries are 0
  */
  bool is_zero();
  /**
    @brief returns the strategies currently in use
  */
  vector<Poly_Sugar_Data *> get_strategies() { return strategies; }
  /** @brief basic properties */
  unsigned number_of_rows() { return A.size(); }
  ///@}
  /** @name Conversion */
  ///@{
  /**
    @brief converts @c this to a vector of Constant_Polynomial
      and returns the result
  */
  vector<Constant_Polynomial *> finalize();
  ///@}
  /** @name Modification */
  ///@{
  /** @brief clears the strategy; do this if you have saved it elsewhere */
  void clear_strategy(unsigned i) { strategies[i] = nullptr; }
  ///@}
  /** @name Computation */
  ///@{
  /** @brief reduces polynomials */
  void reduce_by_old();
  /** @brief reduces polynomials */
  void reduce_by_new();
  ///@}
  /** @name I/O */
  ///@{
  /** @brief lists the reducers selected for each column, in order */
  void list_reducers();
  /**
    @brief prints the matrix
    @param show_data whether to show the monomials that correspond to each column
  */
  void print_matrix(bool show_data=false);
  /**
    @brief prints the indicated row of the matrix as an unformatted polynomial
  **/
  void print_row(unsigned);
  /**
    @brief prints the pairs of monomials and reducers that are being built
  **/
  void print_builder();
  ///@}
protected:
  /*void check_consistency() {
    for (unsigned i = 0; i < num_rows; ++i) {
      const auto & Ai = A[i];
      unsigned n = 0;
      for (auto a : Ai)
        if (a != 0) ++n;
      if (n != nonzero_entries[i]) cout << "row " << i << " inconsistent\n";
    }
  }*/
  /**
    @brief creates rows of the matrix indexed by the specified pairs,
      starting at the specified row
  */
  void initialize_some_rows(const list<Critical_Pair_Basic *> &, unsigned);
  /**
    @brief reduces the specified set of rows, suitable for multithreading
    @details Under the current implementation, we need two buffers for the rows
        of a matrix. We intend to change that eventually to 1.
    @param A_buffer buffer to store the dense representation of a row
    @param R_buffer buffer to store the dense representation of a row
    @param rows which rows this thread will reduce
  */
  void reduce_my_rows(
      vector<COEF_TYPE> & A_buffer,
      vector<COEF_TYPE> & R_buffer,
      const set<unsigned> & rows
  );
  /**
    @brief reduces the specified set of rows by the specified row,
      suitable for multithreading
    @param buffer storage for the dense representation of a row
    @param i which row to use as reductor
    @param F the ground field, needed not just for the modulus but for inverses
    @param rows which rows this thread will reduce
  */
  void reduce_my_new_rows(
      vector<COEF_TYPE> & buffer,
      unsigned i,
      const Prime_Field & F,
      const set<unsigned> & rows
  );
  /**
    @brief convert a row from sparse to dense reprsentation
    @param target dense representation of row
    @param source sparse representation of row
  */
  void densify_row(
      vector<COEF_TYPE> & target,
      const vector<pair<unsigned, COEF_TYPE> > & source
  );
  /**
    @brief reduce dense row using the sparse source, modulo the indicated value
    @param target dense representation of a polynomial being reduced
    @param source sparse representation of a reductor
    @param mod modulus for all operations
  */
  void reduce_buffer(
      vector<COEF_TYPE> & target,
      const vector<pair<unsigned, COEF_TYPE> > & source,
      COEF_TYPE mod
  );
  /**
    @brief transform row from dense to sparse format
    @param buffer dense representation of row
    @param source sparse representation of row
  */
  void sparsify_row(
      const vector<COEF_TYPE> & buffer, vector<pair<unsigned, COEF_TYPE> > & source
  );
  /**
    @brief reduce dense row using the sparse source, modulo the indicated value,
    starting at the indicated position
    @details The difference between this and reduce_buffer is that in this case
    you may not want to start from the first position, because you haven't
    yet added the value \$ a \$. Thus, you can start with the second term,
    rather than waste time (however little) starting with the first.
    @param target dense representation of a polynomial being reduced
    @param source sparse representation of a reductor
    @param mod modulus for all operations
    @param a scalar multiple of source
  */
  void add_reductor(
      vector<COEF_TYPE> & target, COEF_TYPE a,
      const vector<pair<unsigned, COEF_TYPE> > & source, COEF_TYPE mod
  );
  /**
    @brief build the sparse representation of the @b reduced polynomial chosen
        to reduce monomial @p i
    @param buffer where to store the dense representation while building the
        reductor
    @param i index in @c M of the monomial to be reduced
  */
  void build_reducer(vector<COEF_TYPE> & buffer, unsigned i);
  /** @brief number of columns in the polynomial */
  unsigned num_cols;
  /** @brief number of rows in the matrix */
  unsigned num_rows;
  /** @brief monomials for each matrix */
  vector<Monomial *> M;
  /** @brief hash table of the monomials in @c M */
  F4_Hash M_table;
  /** @brief locks on the reducers */
  vector<mutex> red_lock;
  /**
    @brief coefficient data in sparse representation; each vector entry is a
      subrow of the dense matrix
  */
  vector<vector<pair<unsigned, COEF_TYPE> > > A;
  /** @brief buffers for rows of A during reduction; one per thread */
  vector<vector<COEF_TYPE> > A_buf;
  /** @brief storage of monomials and reducers while preprocessing */
  map<Monomial *, Abstract_Polynomial *, MonCmp> M_builder;
  /** @brief finalized list of indices of reducers for the corresponding monomials of @c f */
  vector<Abstract_Polynomial *> R;
  /** @brief reducers actually generated */
  vector<vector<pair<unsigned, COEF_TYPE> > > R_built;
  /** @brief buffer to build reducers for the matrix */
  vector<vector<COEF_TYPE> > R_buf;
  /** @brief count of threads that have actually read a generated actually built */
  vector<unsigned> num_readers;
  /** @brief current basis of ideal */
  const list<Abstract_Polynomial *> & G;
  /** @brief how the monomials are ordered */
  const Monomial_Ordering * mord;
  /** @brief polynomial ring */
  Polynomial_Ring & Rx;
  /** @brief strategy data for each polynomial */
  vector<Poly_Sugar_Data *> strategies;
};

#endif