PolyBoRi
|
00001 // -*- c++ -*- 00002 //***************************************************************************** 00014 //***************************************************************************** 00015 00016 #ifndef polybori_groebner_linear_algebra_step_h_ 00017 #define polybori_groebner_linear_algebra_step_h_ 00018 00019 // include basic definitions 00020 #include "groebner_defs.h" 00021 #include "add_up.h" 00022 #include "nf.h" 00023 #include "draw_matrix.h" 00024 00025 #include "GroebnerStrategy.h" 00026 #include "MatrixMonomialOrderTables.h" 00027 #include "PolyMonomialPairComparerLess.h" 00028 00029 #include "BitMask.h" 00030 #include "PseudoLongProduct.h" 00031 #include "Long64From32BitsPair.h" 00032 00033 #ifdef PBORI_HAVE_NTL 00034 #include <NTL/GF2.h> 00035 #include <NTL/mat_GF2.h> 00036 NTL_CLIENT 00037 #endif 00038 #ifdef PBORI_HAVE_M4RI 00039 const int M4RI_MAXKAY = 16; 00040 #endif 00041 00042 #include <vector> 00043 #include <utility> 00044 #include <sstream> 00045 00046 BEGIN_NAMESPACE_PBORIGB 00047 00048 inline int 00049 select_largest_degree(const ReductionStrategy& strat, const Monomial& m){ 00050 MonomialSet ms=strat.leadingTerms.divisorsOf(m); 00051 if (ms.isZero()) 00052 return -1; 00053 else { 00054 //Monomial min=*(std::min_element(ms.begin(),ms.end(), LessWeightedLengthInStrat(strat))); 00055 Exponent min=*(std::min_element(ms.expBegin(),ms.expEnd(), LargerDegreeComparer())); 00056 return strat.index(min); 00057 } 00058 } 00059 00060 00061 #if defined(PBORI_HAVE_NTL) || defined(PBORI_HAVE_M4RI) 00062 00063 typedef Exponent::idx_map_type from_term_map_type; 00064 00065 00066 00067 inline void 00068 fix_point_iterate(const GroebnerStrategy& strat,std::vector<Polynomial> extendable_system, std::vector<Polynomial>& res1,MonomialSet& res_terms,MonomialSet& leads_from_strat){ 00069 00070 BoolePolyRing current_ring(res_terms.ring()); 00071 leads_from_strat=MonomialSet(current_ring); 00072 res_terms=MonomialSet(current_ring); 00073 00074 for(std::size_t i=0;i<extendable_system.size();i++){ 00075 Polynomial p=extendable_system[i]; 00076 PBORI_ASSERT(p.ring().id() == current_ring.id()); 00077 00078 if PBORI_UNLIKELY(p.isZero()) continue; 00079 00080 p=cheap_reductions(strat.generators, p); 00081 00082 // Polynomial p=mod_mon_set(p_orig.diagram(),strat.generators.monomials); 00083 // if (strat.generators.optLL){ 00084 // Polynomial p_bak2=p; 00085 // p=ll_red_nf(p,strat.generators.llReductor); 00086 // if (p!=p_bak2) p=mod_mon_set(p.diagram(),strat.generators.monomials); 00087 // } 00088 MonomialSet new_terms=p.diagram().diff(res_terms); 00089 MonomialSet::const_iterator it=new_terms.begin(); 00090 MonomialSet::const_iterator end=new_terms.end(); 00091 while(it!=end){ 00092 Monomial m=*it; 00093 00094 int index=select_largest_degree(strat.generators, m); 00095 if PBORI_LIKELY(index>=0){ 00096 00097 Monomial m2=m/strat.generators[index].lead; 00098 Polynomial p2=m2*strat.generators[index].p; 00099 extendable_system.push_back(p2); 00100 PBORI_ASSERT(current_ring.id() == strat.generators[index].lead.ring().id()); 00101 PBORI_ASSERT(current_ring.id() == strat.generators[index].p.ring().id()); 00102 PBORI_ASSERT(current_ring.id() == m.ring().id()); 00103 PBORI_ASSERT(current_ring.id() == m2.ring().id()); 00104 PBORI_ASSERT(current_ring.id() == p2.ring().id()); 00105 } 00106 ++it; 00107 } 00108 res_terms=res_terms.unite(new_terms); 00109 res1.push_back(p); 00110 } 00111 leads_from_strat=res_terms.diff(mod_mon_set(res_terms,strat.generators.minimalLeadingTerms)); 00112 } 00113 00114 inline void 00115 fill_matrix(mzd_t* mat,std::vector<Polynomial> polys, from_term_map_type from_term_map){ 00116 for(std::size_t i=0;i<polys.size();i++){ 00117 Polynomial::exp_iterator it=polys[i].expBegin();//not order dependend 00118 Polynomial::exp_iterator end=polys[i].expEnd(); 00119 while(it!=end){ 00120 #ifndef PBORI_HAVE_M4RI 00121 mat[i][from_term_map[*it]]=1; 00122 #else 00123 from_term_map_type::const_iterator from_it=from_term_map.find(*it); 00124 PBORI_ASSERT(from_it!=from_term_map.end()); 00125 mzd_write_bit(mat,i,from_it->second,1); 00126 #endif 00127 it++; 00128 } 00129 } 00130 } 00131 00132 inline void 00133 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){ 00134 int cols=mat->ncols; 00135 // int rows=mat->nrows; /// @todo unused? 00136 00137 int i; 00138 for(i=0;i<rank;i++){ 00139 int j; 00140 std::vector<int> p_t_i; 00141 00142 bool from_strat=false; 00143 for(j=0;j<cols;j++){ 00144 #ifndef PBORI_HAVE_M4RI 00145 if (mat[i][j]==1){ 00146 #else 00147 if PBORI_UNLIKELY(mzd_read_bit(mat,i,j)==1){ 00148 #endif 00149 if (p_t_i.size()==0){ 00150 if (leads_from_strat.owns(terms_as_exp[j])) { 00151 from_strat=true;break; 00152 } 00153 } 00154 p_t_i.push_back(ring_order2lex[j]); 00155 } 00156 } 00157 if (!(from_strat)){ 00158 std::vector<Exponent> p_t(p_t_i.size()); 00159 std::sort(p_t_i.begin(),p_t_i.end(),std::less<int>()); 00160 for(std::size_t j=0;j<p_t_i.size();j++){ 00161 p_t[j]=terms_as_exp_lex[p_t_i[j]]; 00162 } 00163 polys.push_back(add_up_lex_sorted_exponents(leads_from_strat.ring(), 00164 p_t,0,p_t.size())); 00165 PBORI_ASSERT(!(polys[polys.size()-1].isZero())); 00166 } 00167 } 00168 } 00169 00170 00171 inline void 00172 linalg_step(std::vector<Polynomial>& polys, MonomialSet terms,MonomialSet leads_from_strat, bool log, bool optDrawMatrices=false, const char* matrixPrefix="mat"){ 00173 if PBORI_UNLIKELY(polys.size()==0) return; 00174 00175 static int round=0; 00176 00177 int rows=polys.size(); 00178 int cols=terms.size(); 00179 if PBORI_UNLIKELY(log){ 00180 std::cout<<"ROWS:"<<rows<<"COLUMNS:"<<cols<<std::endl; 00181 } 00182 #ifndef PBORI_HAVE_M4RI 00183 mat_GF2 mat(INIT_SIZE,rows,cols); 00184 #else 00185 mzd_t* mat=mzd_init(rows,cols); 00186 #endif 00187 MatrixMonomialOrderTables tabs(terms); 00188 00189 fill_matrix(mat,polys,tabs.from_term_map); 00190 00191 polys.clear(); 00192 #ifndef PBORI_HAVE_M4RI 00193 int rank=gauss(mat); 00194 #else 00195 if PBORI_UNLIKELY(optDrawMatrices){ 00196 ++round; 00197 std::ostringstream matname; 00198 matname << matrixPrefix << round << ".png"; 00199 draw_matrix(mat, matname.str().c_str()); 00200 } 00201 int rank=mzd_echelonize_m4ri(mat, TRUE, 0);//optimal_k_for_gauss(mat->nrows,mat->ncols,strat)); 00202 #endif 00203 if PBORI_UNLIKELY(log){ 00204 std::cout<<"finished gauss"<<std::endl; 00205 } 00206 translate_back(polys, leads_from_strat, mat,tabs.ring_order2lex, tabs.terms_as_exp,tabs.terms_as_exp_lex,rank); 00207 00208 #ifdef PBORI_HAVE_M4RI 00209 mzd_free(mat); 00210 #endif 00211 } 00212 00213 inline void 00214 printPackedMatrixMB(mzd_t* mat){ 00215 int i,j; 00216 for(i=0;i<mat->nrows;i++){ 00217 for(j=0;j<mat->ncols;j++){ 00218 std::cout<<(int)mzd_read_bit(mat,i,j); 00219 } 00220 std::cout<<std::endl; 00221 } 00222 } 00223 00224 inline mzd_t* 00225 transposePackedMB(mzd_t* mat){ 00226 return mzd_transpose(NULL,mat); 00227 /*mzd_t* res=mzd_init(mat->ncols,mat->nrows); 00228 int i,j; 00229 for(i=0;i<mat->nrows;i++){ 00230 for(j=0;j<mat->ncols;j++){ 00231 mzd_write_bit(res,j,i,mzd_read_bit(mat,i,j)); 00232 } 00233 } 00234 return res;*/ 00235 } 00236 00237 inline mzd_t* 00238 pbori_transpose(mzd_t* mat) { 00239 00240 if PBORI_UNLIKELY(mat->nrows == 0) 00241 return mzd_init(mat->ncols, 0); 00242 00243 if PBORI_UNLIKELY(mat->ncols == 0) 00244 return mzd_init(0, mat->nrows); 00245 00246 return mzd_transpose(NULL,mat); 00247 } 00248 00249 00250 00251 inline void 00252 linalg_step_modified(std::vector < Polynomial > &polys, MonomialSet terms, MonomialSet leads_from_strat, bool log, bool optDrawMatrices, const char* matrixPrefix) 00253 { 00254 BoolePolyRing current_ring(terms.ring()); 00255 PBORI_ASSERT(current_ring.id() == leads_from_strat.ring().id()); 00256 #ifndef PBORI_NDEBUG 00257 std::vector < Polynomial >::const_iterator start(polys.begin()), 00258 finish(polys.end()); 00259 00260 while (start != finish) { 00261 PBORI_ASSERT(current_ring.id() == start->ring().id()); 00262 ++start; 00263 } 00264 #endif 00265 00266 int unmodified_rows=polys.size(); 00267 int unmodified_cols=terms.size(); 00268 00270 if (PBORI_UNLIKELY( (PseudoLongProduct(unmodified_cols, unmodified_rows) > 00271 Long64From32BitsPair<4u, 2820130816u>::get()) )){ 00272 PBoRiError error(CTypes::matrix_size_exceeded); 00273 throw error; 00274 } 00275 00276 static int round=0; 00277 round++; 00278 // const int russian_k=16; /// 00279 MonomialSet terms_unique(current_ring); 00280 std::vector < Monomial > terms_unique_vec; 00281 MonomialSet terms_step1(current_ring); 00282 std::vector < std::pair < Polynomial, Monomial > >polys_lm; 00283 00284 for (std::size_t i = 0; i < polys.size(); i++) { 00285 if PBORI_LIKELY(!(polys[i].isZero())) 00286 polys_lm.push_back(std::pair < Polynomial, Monomial > (polys[i], polys[i].lead())); 00287 } 00288 std:: sort(polys_lm.begin(), polys_lm.end(), PolyMonomialPairComparerLess()); 00289 polys.clear(); 00290 00291 //special cases 00292 if PBORI_UNLIKELY(polys_lm.size() == 0) 00293 return; 00294 Monomial last(current_ring); 00295 if PBORI_UNLIKELY(polys_lm[0].second.deg() == 0) { 00296 PBORI_ASSERT(polys_lm[0].first.isOne()); 00297 polys.resize(1, current_ring); 00298 polys[0] = 1; 00299 00300 return; 00301 } 00302 00303 std::vector < Polynomial > polys_triangular; 00304 std::vector < Polynomial > polys_rest; 00305 00306 { 00307 std::vector < std::pair < Polynomial, Monomial > >::iterator it = polys_lm.begin(); 00308 std::vector < std::pair < Polynomial, Monomial > >::iterator end = polys_lm.end(); 00309 00310 while (it != end) { 00311 if PBORI_LIKELY(it->second != last) { 00312 last = it->second; 00313 polys_triangular.push_back(it->first); 00314 00315 00316 PBORI_ASSERT(std:: find(terms_unique_vec.begin(), terms_unique_vec.end(), it->second) == terms_unique_vec.end()); 00317 00318 terms_unique_vec.push_back(it->second); 00319 terms_step1=terms_step1.unite(it->first.diagram()); 00320 } else 00321 polys_rest.push_back(it->first); 00322 it++; 00323 } 00324 } 00325 polys.clear(); 00326 std::reverse(polys_triangular.begin(), polys_triangular.end()); 00327 terms_unique = add_up_generic(terms_unique_vec, current_ring.zero()); 00328 PBORI_ASSERT(terms_step1.diff(terms).isZero()); 00329 PBORI_ASSERT(polys_triangular.size()!=0); 00330 from_term_map_type eliminated2row_number; 00331 int remaining_cols; 00332 mzd_t* mat_step1; 00333 std::vector<int> compactified_columns2old_columns; 00334 int rows_step1; 00335 std::vector<int> row_start; 00336 //std::vector<Polynomial> result; 00337 MatrixMonomialOrderTables step1(terms_step1); 00338 //std::vector<Exponent> terms_as_exp_step1; 00339 { 00340 int rows=polys_triangular.size(); 00341 int cols=terms_step1.size(); 00342 rows_step1=rows; 00343 if PBORI_UNLIKELY(log){ 00344 std::cout<<"STEP1: ROWS:"<<rows<<"COLUMNS:"<<cols<<std::endl; 00345 } 00346 00347 mat_step1=mzd_init(rows,cols); 00348 //vector<Exponent> terms_as_exp_lex_step1; 00349 //vector<int> ring_order2lex_step1; 00350 //vector<int> lex_order2ring_step1; 00351 // from_term_map_type from_term_map_step1; 00352 //setup_order_tables(terms_as_exp_step1,terms_as_exp_lex_step1,ring_order2lex_step1,lex_order2ring_step1,from_term_map_step1, terms_step1); 00353 fill_matrix(mat_step1,polys_triangular,step1.from_term_map); 00354 00355 polys_triangular.clear(); 00356 00357 if PBORI_UNLIKELY(optDrawMatrices) { 00358 std::ostringstream matname; 00359 matname << matrixPrefix << round << "_step1.png"; 00360 draw_matrix(mat_step1, matname.str().c_str()); 00361 } 00362 //optimize: call back subst directly 00363 mzd_top_echelonize_m4ri 00364 (mat_step1,0); 00365 00366 if PBORI_UNLIKELY(log){ 00367 std::cout<<"finished gauss"<<std::endl; 00368 } 00369 int rank=mat_step1->nrows; 00370 00371 //sort rows 00372 int pivot_row=0; 00373 row_start.resize(rows); 00374 PBORI_ASSERT(cols>=rows); 00375 remaining_cols=cols-rows; 00376 compactified_columns2old_columns.resize(remaining_cols); 00377 for(int i=0;i<cols;i++){ 00378 int j=pivot_row; 00379 for(;j<rows;j++){ 00380 if PBORI_UNLIKELY(mzd_read_bit(mat_step1,j,i)==1){ 00381 if (j!=pivot_row) 00382 mzd_row_swap(mat_step1,j,pivot_row); 00383 00384 eliminated2row_number[step1.terms_as_exp[i]]=pivot_row; 00385 row_start[pivot_row]=i; 00386 pivot_row++; 00387 00388 break; 00389 } 00390 } 00391 if PBORI_UNLIKELY(j==rows){ 00392 PBORI_ASSERT(i>=pivot_row); 00393 compactified_columns2old_columns[i-pivot_row]=i; 00394 } 00395 00396 } 00397 if PBORI_UNLIKELY(log){ 00398 std::cout<<"finished sort"<<std::endl; 00399 } 00400 PBORI_ASSERT(pivot_row==rows); 00401 00402 translate_back(polys, leads_from_strat, mat_step1,step1.ring_order2lex, step1.terms_as_exp,step1.terms_as_exp_lex,rank); 00403 00404 if PBORI_UNLIKELY(log){ 00405 std::cout<<"finished translate"<<std::endl; 00406 } 00407 00408 //delete columns 00409 mzd_t* transposed_step1 = pbori_transpose(mat_step1); 00410 if PBORI_UNLIKELY(log){ 00411 std::cout<<"finished transpose"<<std::endl; 00412 } 00413 mzd_free(mat_step1); 00414 00415 for(int i=0;i<remaining_cols;i++){ 00416 int source=compactified_columns2old_columns[i]; 00417 PBORI_ASSERT(i<=source); 00418 PBORI_ASSERT(source<=transposed_step1->nrows); 00419 if PBORI_LIKELY(i!=source) mzd_row_swap(transposed_step1,source,i); 00420 00421 } 00422 if PBORI_UNLIKELY(log){ 00423 std::cout<<"finished permute"<<std::endl; 00424 } 00425 00426 //cols, rows arguments are swapped, as matrix is transposed 00427 mzd_t* sub_step1=mzd_submatrix(NULL,transposed_step1,0,0,remaining_cols,rows); 00428 00429 if PBORI_UNLIKELY(log){ 00430 std::cout<<"finished submat"<<std::endl; 00431 } 00432 mzd_free(transposed_step1); 00433 mat_step1 = pbori_transpose(sub_step1); 00434 if PBORI_UNLIKELY(log){ 00435 std::cout<<"finished transpose"<<std::endl; 00436 } 00437 mzd_free(sub_step1); 00438 } 00439 MonomialSet terms_step2=terms.diff(terms_unique); 00440 const int rows_step2=polys_rest.size(); 00441 const int cols_step2=terms_step2.size(); 00442 mzd_t* mat_step2=mzd_init(rows_step2,cols_step2); 00443 mzd_t* mat_step2_factor=mzd_init(rows_step2,mat_step1->nrows); 00444 MatrixMonomialOrderTables step2(terms_step2); 00445 // vector<Exponent> step2.terms_as_exp; 00446 // vector<Exponent> step2.terms_as_exp_lex; 00447 // vector<int> step2.ring_order2lex; 00448 // vector<int> step2.lex_order2ring; 00449 // from_term_map_type step2.from_term_map; 00450 // setup_order_tables(step2.terms_as_exp,step2.terms_as_exp_lex,step2.ring_order2lex,step2.lex_order2ring,step2.from_term_map, terms_step2); 00451 00452 00453 for(std::size_t i=0;i<polys_rest.size();i++){ 00454 Polynomial p_r=polys_rest[i]; 00455 Polynomial p_t=p_r.diagram().intersect(terms_step2); 00456 Polynomial p_u=p_r.diagram().diff(p_t.diagram()); 00457 Polynomial::exp_iterator it(p_u.expBegin()), end(p_u.expEnd()); 00458 00459 while(it!=end){ 00460 Exponent e=*it; 00461 from_term_map_type::const_iterator from_it=eliminated2row_number.find(e); 00462 PBORI_ASSERT(step1.terms_as_exp[row_start[from_it->second]]==e); 00463 PBORI_ASSERT(from_it!=eliminated2row_number.end()); 00464 int index=from_it->second;//...translate e->line number; 00465 mzd_write_bit(mat_step2_factor,i,index,1); 00466 it++; 00467 } 00468 it=p_t.expBegin(); 00469 end=p_t.expEnd(); 00470 while(it!=end){ 00471 Exponent e=*it; 00472 from_term_map_type::const_iterator from_it=step2.from_term_map.find(e); 00473 PBORI_ASSERT(from_it!=step2.from_term_map.end()); 00474 int index=from_it->second; 00475 mzd_write_bit(mat_step2,i,index,1); 00476 it++; 00477 } 00478 } 00479 00480 if PBORI_UNLIKELY(log){ 00481 std::cout<<"iterate over rest polys"<<std::endl; 00482 } 00483 00484 std::vector<int> remaining_col2new_col(remaining_cols); 00485 for(int i=0;i<remaining_cols;i++){ 00486 remaining_col2new_col[i]=step2.from_term_map[step1.terms_as_exp[compactified_columns2old_columns[i]]]; 00487 } 00488 PBORI_ASSERT(mat_step2_factor->ncols==mat_step1->nrows); 00489 PBORI_ASSERT(mat_step1->nrows==terms_unique.size()); 00490 PBORI_ASSERT(mat_step1->ncols==remaining_cols); 00491 00492 if PBORI_UNLIKELY(optDrawMatrices) 00493 { 00494 std::ostringstream matname; 00495 matname << matrixPrefix << round << "_mult_A.png"; 00496 draw_matrix(mat_step2_factor, matname.str().c_str()); 00497 matname.clear(); 00498 matname << mat_step2_factor << round << "_mult_B.png"; 00499 draw_matrix(mat_step1,matname.str().c_str()); 00500 } 00501 if PBORI_UNLIKELY(log){ 00502 std::cout<<"start mult"<<std::endl; 00503 } 00504 00505 mzd_t* eliminated; 00506 if PBORI_LIKELY((mat_step2_factor->nrows!=0) && (mat_step1->nrows!=0) && (mat_step2_factor->ncols!=0) && (mat_step1->ncols!=0)) 00507 eliminated=mzd_mul_m4rm(NULL,mat_step2_factor,mat_step1,0); 00508 else 00509 eliminated=mzd_init(mat_step2_factor->nrows,mat_step1->ncols); 00510 00511 mzd_free(mat_step2_factor); 00512 if PBORI_UNLIKELY(log){ 00513 std::cout<<"end mult"<<std::endl; 00514 } 00515 mzd_free(mat_step1); 00516 PBORI_ASSERT(polys_rest.size()==eliminated->nrows); 00517 PBORI_ASSERT(mat_step2->nrows==eliminated->nrows); 00518 for(std::size_t i=0;i<polys_rest.size();i++){ 00519 PBORI_ASSERT(remaining_cols==eliminated->ncols); 00520 for(int j=0;j<remaining_cols;j++){ 00521 if PBORI_UNLIKELY(mzd_read_bit(eliminated,i,j)==1){ 00522 PBORI_ASSERT(step2.terms_as_exp[remaining_col2new_col[j]]==step1.terms_as_exp[compactified_columns2old_columns[j]]); 00523 00524 if PBORI_UNLIKELY(mzd_read_bit(mat_step2,i,remaining_col2new_col[j])==1){ 00525 mzd_write_bit(mat_step2,i,remaining_col2new_col[j],0); 00526 } else mzd_write_bit(mat_step2,i,remaining_col2new_col[j],1); 00527 } 00528 } 00529 } 00530 00531 mzd_free(eliminated); 00532 00533 if PBORI_UNLIKELY(log){ 00534 std::cout<<"STEP2: ROWS:"<<rows_step2<<"COLUMNS:"<<cols_step2<<std::endl; 00535 } 00536 if PBORI_UNLIKELY(optDrawMatrices) 00537 { 00538 std::ostringstream matname; 00539 matname << matrixPrefix << round << "_step2.png"; 00540 draw_matrix(mat_step2, matname.str().c_str()); 00541 } 00542 00543 00544 int rank_step2; 00545 if ((mat_step2->ncols>0) &&( mat_step2->nrows>0)){ 00546 rank_step2= 00547 mzd_echelonize_m4ri(mat_step2,TRUE,0); 00548 //mzd_echelonize_pluq(mat_step2,TRUE); 00549 } else 00550 rank_step2=0; 00551 00552 if PBORI_UNLIKELY(log){ 00553 std::cout<<"finished gauss"<<std::endl; 00554 } 00555 00556 translate_back(polys, leads_from_strat, mat_step2,step2.ring_order2lex, step2.terms_as_exp,step2.terms_as_exp_lex,rank_step2); 00557 mzd_free(mat_step2); 00558 00559 00560 } 00561 00562 00563 inline std::vector<Polynomial> 00564 gauss_on_polys(const std::vector<Polynomial>& orig_system){ 00565 00566 if (orig_system.empty()) 00567 return orig_system; 00568 00569 Polynomial init(0, orig_system[0].ring()); 00570 MonomialSet terms=unite_polynomials(orig_system, init); 00571 MonomialSet from_strat(init.ring());//no strat 00572 std::vector<Polynomial> polys(orig_system); 00573 linalg_step(polys, terms, from_strat, false); 00574 return polys; 00575 } 00576 #endif 00577 00578 END_NAMESPACE_PBORIGB 00579 00580 #endif /* polybori_groebner_linear_algebra_step_h_ */