16 #ifndef polybori_groebner_linear_algebra_step_h_
17 #define polybori_groebner_linear_algebra_step_h_
35 #include <NTL/mat_GF2.h>
38 #ifdef PBORI_HAVE_M4RI
39 const int M4RI_MAXKAY = 16;
56 return strat.
index(min);
61 #if defined(PBORI_HAVE_NTL) || defined(PBORI_HAVE_M4RI)
63 typedef Exponent::idx_map_type from_term_map_type;
68 fix_point_iterate(
const GroebnerStrategy& strat,std::vector<Polynomial> extendable_system, std::vector<Polynomial>& res1,
MonomialSet& res_terms,
MonomialSet& leads_from_strat){
74 for(std::size_t i=0;i<extendable_system.size();i++){
89 MonomialSet::const_iterator it=new_terms.
begin();
90 MonomialSet::const_iterator end=new_terms.end();
97 Monomial m2=m/strat.generators[index].lead;
99 extendable_system.push_back(p2);
100 PBORI_ASSERT(current_ring.id() == strat.generators[index].lead.ring().id());
101 PBORI_ASSERT(current_ring.id() == strat.generators[index].p.ring().id());
108 res_terms=res_terms.unite(new_terms);
111 leads_from_strat=res_terms.diff(
mod_mon_set(res_terms,strat.generators.minimalLeadingTerms));
115 fill_matrix(mzd_t* mat,std::vector<Polynomial> polys, from_term_map_type from_term_map){
116 for(std::size_t i=0;i<polys.size();i++){
117 Polynomial::exp_iterator it=polys[i].expBegin();
118 Polynomial::exp_iterator end=polys[i].expEnd();
120 #ifndef PBORI_HAVE_M4RI
121 mat[i][from_term_map[*it]]=1;
123 from_term_map_type::const_iterator from_it=from_term_map.find(*it);
125 mzd_write_bit(mat,i,from_it->second,1);
133 translate_back(std::vector<Polynomial>& polys,
MonomialSet leads_from_strat,mzd_t* mat,
const std::vector<int>& ring_order2lex,
const std::vector<Exponent>& terms_as_exp,
const std::vector<Exponent>& terms_as_exp_lex,
int rank){
140 std::vector<int> p_t_i;
142 bool from_strat=
false;
144 #ifndef PBORI_HAVE_M4RI
149 if (p_t_i.size()==0){
150 if (leads_from_strat.owns(terms_as_exp[j])) {
151 from_strat=
true;
break;
154 p_t_i.push_back(ring_order2lex[j]);
158 std::vector<Exponent> p_t(p_t_i.size());
159 std::sort(p_t_i.begin(),p_t_i.end(),std::less<int>());
160 for(std::size_t j=0;j<p_t_i.size();j++){
161 p_t[j]=terms_as_exp_lex[p_t_i[j]];
172 linalg_step(std::vector<Polynomial>& polys,
MonomialSet terms,
MonomialSet leads_from_strat,
bool log,
bool optDrawMatrices=
false,
const char* matrixPrefix=
"mat"){
177 int rows=polys.size();
178 int cols=terms.size();
180 std::cout<<
"ROWS:"<<rows<<
"COLUMNS:"<<cols<<std::endl;
182 #ifndef PBORI_HAVE_M4RI
183 mat_GF2 mat(INIT_SIZE,rows,cols);
185 mzd_t* mat=mzd_init(rows,cols);
187 MatrixMonomialOrderTables tabs(terms);
189 fill_matrix(mat,polys,tabs.from_term_map);
192 #ifndef PBORI_HAVE_M4RI
197 std::ostringstream matname;
198 matname << matrixPrefix << round <<
".png"<< std::ends;
201 int rank=mzd_echelonize_m4ri(mat, TRUE, 0);
204 std::cout<<
"finished gauss"<<std::endl;
206 translate_back(polys, leads_from_strat, mat,tabs.ring_order2lex, tabs.terms_as_exp,tabs.terms_as_exp_lex,rank);
208 #ifdef PBORI_HAVE_M4RI
214 printPackedMatrixMB(mzd_t* mat){
216 for(i=0;i<mat->nrows;i++){
217 for(j=0;j<mat->ncols;j++){
218 std::cout<<(int)mzd_read_bit(mat,i,j);
220 std::cout<<std::endl;
225 transposePackedMB(mzd_t* mat){
226 return mzd_transpose(NULL,mat);
238 pbori_transpose(mzd_t* mat) {
241 return mzd_init(mat->ncols, 0);
244 return mzd_init(0, mat->nrows);
246 return mzd_transpose(NULL,mat);
254 BoolePolyRing current_ring(terms.ring());
255 PBORI_ASSERT(current_ring.id() == leads_from_strat.ring().id());
257 std::vector < Polynomial >::const_iterator start(polys.begin()),
260 while (start != finish) {
266 int unmodified_rows=polys.size();
267 int unmodified_cols=terms.size();
270 if (
PBORI_UNLIKELY( (PseudoLongProduct(unmodified_cols, unmodified_rows) >
271 Long64From32BitsPair<4u, 2820130816u>::get()) )){
272 PBoRiError error(CTypes::matrix_size_exceeded);
280 std::vector < Monomial > terms_unique_vec;
282 std::vector < std::pair < Polynomial, Monomial > >polys_lm;
284 for (std::size_t i = 0; i < polys.size(); i++) {
288 std:: sort(polys_lm.begin(), polys_lm.end(), PolyMonomialPairComparerLess());
297 polys.resize(1, current_ring);
303 std::vector < Polynomial > polys_triangular;
304 std::vector < Polynomial > polys_rest;
307 std::vector < std::pair < Polynomial, Monomial > >::iterator it = polys_lm.begin();
308 std::vector < std::pair < Polynomial, Monomial > >::iterator end = polys_lm.end();
313 polys_triangular.push_back(it->first);
316 PBORI_ASSERT(std:: find(terms_unique_vec.begin(), terms_unique_vec.end(), it->second) == terms_unique_vec.end());
318 terms_unique_vec.push_back(it->second);
319 terms_step1=terms_step1.unite(it->first.diagram());
321 polys_rest.push_back(it->first);
326 std::reverse(polys_triangular.begin(), polys_triangular.end());
327 terms_unique =
add_up_generic(terms_unique_vec, current_ring.zero());
330 from_term_map_type eliminated2row_number;
333 std::vector<int> compactified_columns2old_columns;
335 std::vector<int> row_start;
337 MatrixMonomialOrderTables step1(terms_step1);
340 int rows=polys_triangular.size();
341 int cols=terms_step1.size();
344 std::cout<<
"STEP1: ROWS:"<<rows<<
"COLUMNS:"<<cols<<std::endl;
347 mat_step1=mzd_init(rows,cols);
353 fill_matrix(mat_step1,polys_triangular,step1.from_term_map);
355 polys_triangular.clear();
358 std::ostringstream matname;
359 matname << matrixPrefix << round <<
"_step1.png" <<std::ends;
363 mzd_top_echelonize_m4ri
367 std::cout<<
"finished gauss"<<std::endl;
369 int rank=mat_step1->nrows;
373 row_start.resize(rows);
375 remaining_cols=cols-rows;
376 compactified_columns2old_columns.resize(remaining_cols);
377 for(
int i=0;i<cols;i++){
382 mzd_row_swap(mat_step1,j,pivot_row);
384 eliminated2row_number[step1.terms_as_exp[i]]=pivot_row;
385 row_start[pivot_row]=i;
393 compactified_columns2old_columns[i-pivot_row]=i;
398 std::cout<<
"finished sort"<<std::endl;
402 translate_back(polys, leads_from_strat, mat_step1,step1.ring_order2lex, step1.terms_as_exp,step1.terms_as_exp_lex,rank);
405 std::cout<<
"finished translate"<<std::endl;
409 mzd_t* transposed_step1 = pbori_transpose(mat_step1);
411 std::cout<<
"finished transpose"<<std::endl;
415 for(
int i=0;i<remaining_cols;i++){
416 int source=compactified_columns2old_columns[i];
419 if PBORI_LIKELY(i!=source) mzd_row_swap(transposed_step1,source,i);
423 std::cout<<
"finished permute"<<std::endl;
427 mzd_t* sub_step1=mzd_submatrix(NULL,transposed_step1,0,0,remaining_cols,rows);
430 std::cout<<
"finished submat"<<std::endl;
432 mzd_free(transposed_step1);
433 mat_step1 = pbori_transpose(sub_step1);
435 std::cout<<
"finished transpose"<<std::endl;
440 const int rows_step2=polys_rest.
size();
441 const int cols_step2=terms_step2.size();
442 mzd_t* mat_step2=mzd_init(rows_step2,cols_step2);
443 mzd_t* mat_step2_factor=mzd_init(rows_step2,mat_step1->nrows);
444 MatrixMonomialOrderTables step2(terms_step2);
453 for(std::size_t i=0;i<polys_rest.size();i++){
454 Polynomial p_r=polys_rest[i];
455 Polynomial p_t=p_r.diagram().intersect(terms_step2);
456 Polynomial p_u=p_r.diagram().diff(p_t.diagram());
457 Polynomial::exp_iterator it(p_u.expBegin()), end(p_u.expEnd());
461 from_term_map_type::const_iterator from_it=eliminated2row_number.find(e);
462 PBORI_ASSERT(step1.terms_as_exp[row_start[from_it->second]]==e);
464 int index=from_it->second;
465 mzd_write_bit(mat_step2_factor,i,index,1);
472 from_term_map_type::const_iterator from_it=step2.from_term_map.find(e);
474 int index=from_it->second;
475 mzd_write_bit(mat_step2,i,index,1);
481 std::cout<<
"iterate over rest polys"<<std::endl;
484 std::vector<int> remaining_col2new_col(remaining_cols);
485 for(
int i=0;i<remaining_cols;i++){
486 remaining_col2new_col[i]=step2.from_term_map[step1.terms_as_exp[compactified_columns2old_columns[i]]];
488 PBORI_ASSERT(mat_step2_factor->ncols==mat_step1->nrows);
494 std::ostringstream matname;
495 matname << matrixPrefix << round <<
"_mult_A.png"<<std::ends;
496 draw_matrix(mat_step2_factor, matname.str().c_str());
498 matname << mat_step2_factor << round <<
"_mult_B.png"<<std::ends;
502 std::cout<<
"start mult"<<std::endl;
506 if PBORI_LIKELY((mat_step2_factor->nrows!=0) && (mat_step1->nrows!=0) && (mat_step2_factor->ncols!=0) && (mat_step1->ncols!=0))
507 eliminated=mzd_mul_m4rm(NULL,mat_step2_factor,mat_step1,0);
509 eliminated=mzd_init(mat_step2_factor->nrows,mat_step1->ncols);
511 mzd_free(mat_step2_factor);
513 std::cout<<
"end mult"<<std::endl;
518 for(std::size_t i=0;i<polys_rest.size();i++){
520 for(
int j=0;j<remaining_cols;j++){
522 PBORI_ASSERT(step2.terms_as_exp[remaining_col2new_col[j]]==step1.terms_as_exp[compactified_columns2old_columns[j]]);
524 if PBORI_UNLIKELY(mzd_read_bit(mat_step2,i,remaining_col2new_col[j])==1){
525 mzd_write_bit(mat_step2,i,remaining_col2new_col[j],0);
526 }
else mzd_write_bit(mat_step2,i,remaining_col2new_col[j],1);
531 mzd_free(eliminated);
534 std::cout<<
"STEP2: ROWS:"<<rows_step2<<
"COLUMNS:"<<cols_step2<<std::endl;
538 std::ostringstream matname;
539 matname << matrixPrefix << round <<
"_step2.png"<<std::ends;
545 if ((mat_step2->ncols>0) &&( mat_step2->nrows>0)){
547 mzd_echelonize_m4ri(mat_step2,TRUE,0);
553 std::cout<<
"finished gauss"<<std::endl;
556 translate_back(polys, leads_from_strat, mat_step2,step2.ring_order2lex, step2.terms_as_exp,step2.terms_as_exp_lex,rank_step2);
563 inline std::vector<Polynomial>
564 gauss_on_polys(
const std::vector<Polynomial>& orig_system){
566 if (orig_system.empty())
569 Polynomial init(0, orig_system[0].ring());
572 std::vector<Polynomial> polys(orig_system);
573 linalg_step(polys, terms, from_strat,
false);
This class is just a wrapper for using variables for storing indices as interim data structure for Bo...
Definition: BooleExponent.h:34
Polynomial add_up_generic(const std::vector< T > &res_vec, int start, int end, Polynomial init)
Definition: add_up.h:205
polybori::BooleExponent Exponent
Definition: groebner_defs.h:33
Polynomial cheap_reductions(const ReductionStrategy &strat, Polynomial p)
Definition: nf.cc:831
size_type index(const KeyType &key) const
Retrieve index associated to key.
Definition: PolyEntryVector.h:99
#define END_NAMESPACE_PBORIGB
Definition: groebner_defs.h:16
const self & diagram() const
Access internal decision diagram.
Definition: BooleSet.h:231
const ring_type & ring() const
Access ring, where this belongs to.
Definition: BoolePolynomial.h:478
self divisorsOf(const term_type &rhs) const
Compute intersection with divisors of rhs.
Definition: BooleSet.cc:156
void draw_matrix(mzd_t *mat, const char *filename)
Definition: draw_matrix.h:36
exp_iterator expEnd() const
Finish of iteration over exponent vectors.
Definition: BooleSet.cc:109
This class reinterprets decicion diagram managers as Boolean polynomial rings, adds an ordering and v...
Definition: BoolePolyRing.h:40
BoolePolynomial Polynomial
Definition: embed.h:51
Polynomial unite_polynomials(const std::vector< Polynomial > &res_vec, int start, int end, Polynomial init)
Definition: add_up.h:153
#define BEGIN_NAMESPACE_PBORIGB
Definition: groebner_defs.h:15
exp_iterator expBegin() const
Start of iteration over exponent vectors.
Definition: BooleSet.cc:101
This class wraps the underlying decicion diagram type and defines the necessary operations.
Definition: BoolePolynomial.h:85
hash_type id() const
Get unique identifier for this ring.
Definition: BoolePolyRing.h:151
#define PBORI_LIKELY(expression)
For optimizing if-branches.
Definition: pbori_defs.h:56
polybori::BooleSet MonomialSet
Definition: groebner_defs.h:45
bool isZero() const
Test whether diagram represents the empty set.
Definition: CCuddDDFacade.h:244
MonomialSet mod_mon_set(const MonomialSet &as, const MonomialSet &vs)
Definition: nf.cc:855
size_type size() const
Returns number of terms.
Definition: BooleSet.h:242
int select_largest_degree(const ReductionStrategy &strat, const Monomial &m)
Definition: linear_algebra_step.h:49
#define PBORI_ASSERT(arg)
Definition: pbori_defs.h:118
const ring_type & ring() const
Access ring, where this belongs to.
Definition: BooleMonomial.h:235
MonomialSet add_up_lex_sorted_exponents(const BoolePolyRing &ring, std::vector< Exponent > &vec, int start, int end)
Definition: add_up.h:61
bool_type isZero() const
Check whether polynomial is constant zero.
Definition: BoolePolynomial.h:294
LeadingTerms leadingTerms
Definition: ReductionTerms.h:51
const_iterator begin() const
Start of iteration over terms.
Definition: BooleSet.cc:70
Definition: BooleSet.h:57
This class defines ReductionStrategy.
Definition: ReductionStrategy.h:34
This class is just a wrapper for using variables from cudd's decicion diagram.
Definition: BooleMonomial.h:50
This class defines LargerDegreeComparer.
Definition: LargerDegreeComparer.h:28
BooleMonomial Monomial
Definition: embed.h:53
#define PBORI_UNLIKELY(expression)
Definition: pbori_defs.h:59