PolyBoRi
linear_algebra_step.h
Go to the documentation of this file.
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_ */