PolyBoRi
minimal_elements.h
Go to the documentation of this file.
00001 // -*- c++ -*-
00002 //*****************************************************************************
00014 //*****************************************************************************
00015 
00016 #ifndef polybori_groebner_minimal_elements_h_
00017 #define polybori_groebner_minimal_elements_h_
00018 
00019 // include definitions
00020 #include "groebner_defs.h"
00021 #include "contained_variables.h"
00022 //#include "nf.h"
00023 
00024 #include <vector>
00025 #include <algorithm>
00026 
00027 BEGIN_NAMESPACE_PBORIGB
00028 MonomialSet mod_mon_set(const MonomialSet& as, const MonomialSet &vs); // nf.cc
00029 MonomialSet mod_deg2_set(const MonomialSet& as, const MonomialSet &vs);  // groebner_alg.cc
00030 MonomialSet mod_var_set(const MonomialSet& as, const MonomialSet& vs);   // groebner_alg.cc
00031 
00032 inline MonomialSet
00033 minimal_elements_internal(const MonomialSet& s){
00034     if (s.isZero()) return s;
00035     if (Polynomial(s).isOne()) return s;
00036     MonomialSet::navigator nav=s.navigation();
00037     int i=*nav;
00038     
00039     
00040     if (Polynomial(s).hasConstantPart()) return MonomialSet(Polynomial(true, s.ring()));
00041     int l=s.length();
00042     if (l<=1) {
00043         return s;
00044     }
00045     
00046     if(l==2){
00047         MonomialSet::const_iterator it=s.begin();
00048         Monomial a=*it;
00049         Monomial b=*(++it);
00050         if (a.reducibleBy(b)) return MonomialSet(b.diagram());
00051         else return s;
00052     }
00053     
00054     MonomialSet s0_raw=s.subset0(i);
00055     MonomialSet s0=minimal_elements_internal(s0_raw);
00056     MonomialSet s1=minimal_elements_internal(s.subset1(i).diff(s0_raw));
00057     if (!(s0.isZero())){
00058       s1=s1.diff(s0.unateProduct(Polynomial(s1).usedVariablesExp().divisors(s.ring())));
00059         
00060     }
00061     return s0.unite(s1.change(i));
00062 
00063 }
00064 
00065 inline MonomialSet
00066 minimal_elements_internal2(MonomialSet s){
00067     if (s.isZero()) return s;
00068     if (Polynomial(s).isOne()) return s;
00069     
00070     
00071     
00072     
00073     if (Polynomial(s).hasConstantPart()) 
00074       return MonomialSet(Polynomial(true, s.ring()));
00075     MonomialSet result(s.ring());
00076     std::vector<idx_type> cv=contained_variables(s);
00077     if ((cv.size()>0) && (s.length()==cv.size())){
00078         return s;
00079     } else {
00080     
00081         MonomialSet::size_type z;
00082         MonomialSet cv_set(s.ring());
00083         for(z=cv.size()-1;z>=0;z--){
00084             Monomial mv=Variable(cv[z], s.ring());
00085             cv_set=cv_set.unite(mv.diagram());
00086         }
00087         for(z=0;z<cv.size();z++){
00088             s=s.subset0(cv[z]);
00089         }
00090         result=cv_set;
00091     }
00092     
00093     if (s.isZero()) return result;
00094     PBORI_ASSERT(!(Polynomial(s).hasConstantPart()));
00095     
00096     
00097     
00098     MonomialSet::navigator nav=s.navigation();
00099     idx_type i;
00100     #if 1
00101     
00102     //first index of ZDD
00103     i=*nav;
00104     #else
00105     
00106     //first index of least lex. term
00107     while(!(nav.isConstant())){
00108         i=*nav;
00109         nav.incrementElse();
00110     }
00111     #endif
00112     
00113     
00114     
00115     /*MonomialSet s0=minimal_elements_internal2(s.subset0(i));
00116     MonomialSet s1=s.subset1(i);
00117     if ((s0!=s1)&&(!(s1.diff(s0).isZero()))){
00118         s1=minimal_elements_internal2(s1.unite(s0)).diff(s0);
00119     } else return s0;
00120     return s0.unite(s1.change(i));*/
00121     
00122     MonomialSet s0_raw=s.subset0(i);
00123     MonomialSet s0=minimal_elements_internal2(s0_raw);
00124     MonomialSet s1=minimal_elements_internal2(s.subset1(i).diff(s0_raw));
00125     if (!(s0.isZero())){
00126         s1=s1.diff(s0.unateProduct(Polynomial(s1).usedVariables().divisors()));
00127         
00128     }
00129     return s0.unite(s1.change(i)).unite(result);
00130 
00131 }
00132 
00133 inline std::vector<Exponent>
00134 minimal_elements_internal3(MonomialSet s){
00135     std::vector<Exponent> result;
00136     if (s.isZero()) return result;
00137     if ((Polynomial(s).isOne()) || (Polynomial(s).hasConstantPart())){
00138         result.push_back(Exponent());
00139         return result;
00140     }
00141     std::vector<idx_type> cv=contained_variables(s);
00142     for(MonomialSet::size_type i=0;i<cv.size();i++){
00143             s=s.subset0(cv[i]);
00144             Exponent t;
00145             t.insert(cv[i]);
00146             result.push_back(t);
00147     }
00148     
00149     if (s.isZero()){
00150         return result;
00151     } else {
00152         std::vector<Exponent> exponents;
00153         //Pol sp=s;
00154         exponents.insert(exponents.end(), s.expBegin(),s.expEnd());
00155         int nvars=s.ring().nVariables();
00156         std::vector<std::vector<int> > occ_vecs(nvars);
00157         for(MonomialSet::size_type i=0;i<exponents.size()-1;i++){
00158             Exponent::const_iterator it=((const Exponent&) exponents[i]).begin();
00159             Exponent::const_iterator end=((const Exponent&) exponents[i]).end();
00160             while(it!=end){
00161                 occ_vecs[*it].push_back(i);
00162                 it++;
00163             }
00164         }
00165         //now exponents are ordered
00166         /*std::vector<std::set<int> > occ_sets(nvars);
00167         for(i=occ_sets.size()-1;i>=0;i--){
00168             occ_sets[i].insert(occ_vecs[i].begin(),occ_vecs[i].end());
00169         }*/
00170         std::vector<bool> still_minimal(exponents.size());
00171         for(MonomialSet::size_type i=exponents.size()-1;i>=0;i--){
00172             still_minimal[i]=true;
00173         }
00174 
00175         //lex smalles is last so backwards
00176         for(MonomialSet::size_type i=exponents.size()-1;i>=0;i--){
00177             if (still_minimal[i]){
00178                 //we assume, that each exponents has deg>0
00179                 Exponent::const_iterator it=((const Exponent&) exponents[i]).begin();
00180                 Exponent::const_iterator end=((const Exponent&) exponents[i]).end();
00181                 int first_index=*it;
00182                 std::vector<int> occ_set=occ_vecs[first_index];
00183                 it++;
00184                 while(it!=end){
00185                     
00186                     std::vector<int> occ_set_next;
00187                     set_intersection(
00188                         occ_set.begin(),
00189                         occ_set.end(),
00190                         occ_vecs[*it].begin(), 
00191                         occ_vecs[*it].end(),
00192                         std::back_insert_iterator<std::vector<int> >(occ_set_next));
00193                     occ_set=occ_set_next;
00194                     it++;
00195                 }
00196                 
00197                 std::vector<int>::const_iterator oc_it= occ_set.begin();
00198                 std::vector<int>::const_iterator  oc_end= occ_set.end();
00199                 while(oc_it!=oc_end){
00200                     still_minimal[*oc_it]=false;
00201                     oc_it++;
00202                 }
00203                 
00204                 it=((const Exponent&) exponents[i]).begin();
00205                 while(it!=end){
00206                     std::vector<int> occ_set_difference;
00207                     set_difference(
00208                         occ_vecs[*it].begin(),
00209                         occ_vecs[*it].end(),
00210                         occ_set.begin(),
00211                         occ_set.end(),
00212                         std::back_insert_iterator<std::vector<int> >(occ_set_difference));
00213                     occ_vecs[*it]=occ_set_difference;
00214                     it++;
00215                 }
00216                 result.push_back(exponents[i]);
00217             }
00218         }
00219         
00220     }
00221     return result;
00222 }
00223 
00224 inline MonomialSet
00225 minimal_elements(const MonomialSet& s){
00226 #if 0
00227     //return minimal_elements_internal2(s);
00228     return minimal_elements_cudd_style_unary(s);
00229 #else
00230 #if 1
00231   return s.minimalElements();
00232 #else
00233   MonomialSet minElts = minimal_elements_internal(s);
00234 
00235   if (minElts!=s.minimalElements()){
00236 
00237     std::cout <<"ERROR minimalElements"<<std::endl;
00238     std::cout <<"orig "<<s<< std::endl;
00239     std::cout <<"correct "<<minElts<< std::endl;
00240     std::cout <<"wrong "<<s.minimalElements()<< std::endl;
00241   }
00242 
00243 
00244   return minElts;
00245 #endif
00246 #endif
00247 }
00248 
00249 
00250 
00251 inline MonomialSet
00252 minimal_elements_cudd_style_unary(MonomialSet m){
00253 
00254   if (m.isZero()) return m;
00255   
00256   if (m.ownsOne()) return Polynomial(1, m.ring()).diagram();
00257 
00258   MonomialSet::navigator m_nav=m.navigation();
00259   MonomialSet::navigator ms0=m_nav.elseBranch();
00260   MonomialSet::navigator ms1=m_nav.thenBranch();
00261 
00262   
00263   typedef PBORI::CacheManager<CCacheTypes::minimal_elements>
00264     cache_mgr_type;
00265 
00266 
00267 
00268   cache_mgr_type cache_mgr(m.ring());
00269   PBORI::BoolePolynomial::navigator cached =
00270     cache_mgr.find(m_nav);
00271 
00272 
00273       
00274   if (cached.isValid() ){
00275     return cache_mgr.generate(cached);
00276   }
00277   
00278   MonomialSet minimal_else=minimal_elements_cudd_style_unary(cache_mgr.generate(ms0));
00279   MonomialSet minimal_then=minimal_elements_cudd_style_unary(mod_mon_set(cache_mgr.generate(ms1),minimal_else));
00280   MonomialSet result(m.ring());
00281   if ((minimal_else.navigation()==ms0) &&(minimal_then.navigation()==ms1)) result=m;
00282   else
00283     result= MonomialSet(*m_nav,minimal_then,minimal_else);//result0.unite(result1.change(index));
00284 
00285   cache_mgr.insert(m_nav, result.navigation());
00286   return result;
00287 }
00288 
00289 inline MonomialSet
00290 do_minimal_elements_cudd_style(MonomialSet m, MonomialSet mod){
00291   Polynomial p_mod=mod;
00292   if (m.isZero()) return m;
00293   if (mod.ownsOne())
00294     return MonomialSet(mod.ring());
00295   if (m.ownsOne()) return Polynomial(1,m.ring()).diagram();
00296   MonomialSet mod_cv=contained_variables_cudd_style(mod);
00297   m=mod_var_set(m,mod_cv);
00298   mod=mod_var_set(mod,mod_cv);
00299   MonomialSet mod_d2=contained_deg2_cudd_style(mod);
00300   mod=mod_deg2_set(mod,mod_d2);
00301   m=mod_deg2_set(m,mod_d2);
00302   MonomialSet cv=contained_variables_cudd_style(m);
00303   MonomialSet cv_orig=cv;
00304   cv=cv.diff(mod);
00305   mod=mod_var_set(mod,cv_orig);
00306   m=mod_var_set(m,cv_orig);
00307   m=m.diff(mod);
00308   if (m.isZero()) return cv;
00309   bool cv_empty=cv.isZero();
00310   
00311   MonomialSet result(m.ring());
00312   int index=*m.navigation();
00313   
00314   
00315   
00316   
00317   if (!mod.isZero())
00318   {
00319     MonomialSet::navigator nav_mod=mod.navigation();
00320     while((!(nav_mod.isConstant())) && (index>*nav_mod)){
00321       nav_mod.incrementElse();
00322      
00323     }
00324     mod=MonomialSet(nav_mod, m.ring());
00325   }
00326   
00327   MonomialSet::navigator m_nav=m.navigation();
00328   MonomialSet::navigator ms0=m_nav.elseBranch();
00329   MonomialSet::navigator ms1=m_nav.thenBranch();
00330   MonomialSet::navigator mod_nav=mod.navigation();
00331   
00332   typedef PBORI::CacheManager<CCacheTypes::minimal_mod>
00333     cache_mgr_type;
00334 
00335 
00336 
00337   cache_mgr_type cache_mgr(m.ring());
00338   PBORI::BoolePolynomial::navigator cached =
00339     cache_mgr.find(m_nav, mod_nav);
00340 
00341 
00342       
00343   if (cached.isValid() ){
00344     return cv.unite((MonomialSet)cache_mgr.generate(cached));
00345   }
00346   
00347   if (mod.isZero()){
00348     
00349     MonomialSet result0=do_minimal_elements_cudd_style(cache_mgr.generate(ms0),
00350                                                        mod); 
00351     MonomialSet result1= do_minimal_elements_cudd_style(cache_mgr.generate(ms1),
00352                                                         result0);
00353     result= MonomialSet(index,result1,result0);//result0.unite(result1.change(index));
00354      
00355   } else {
00356       if (*mod_nav>index){
00357         MonomialSet
00358           result0=do_minimal_elements_cudd_style(cache_mgr.generate(ms0), mod);
00359         MonomialSet result1= do_minimal_elements_cudd_style(
00360           cache_mgr.generate(ms1),result0.unite(mod));
00361         if (result1.isZero()) {result=result0;}
00362         else
00363           {result= MonomialSet(index,result1,result0);}
00364       } else {
00365       PBORI_ASSERT(index==*mod_nav);
00366       MonomialSet::navigator mod0=mod_nav.elseBranch();
00367       MonomialSet::navigator mod1=mod_nav.thenBranch();
00368       MonomialSet
00369         result0=do_minimal_elements_cudd_style(cache_mgr.generate(ms0), cache_mgr.generate(mod0));
00370       //MonomialSet mod1=mod.subset1(index);
00371       MonomialSet result1=
00372         do_minimal_elements_cudd_style(cache_mgr.generate(ms1), 
00373                                        result0.unite(cache_mgr.generate(ms0).unite(cache_mgr.generate(mod1))));
00374       result= MonomialSet(index,result1,result0);//result0.unite(result1.change(index));
00375     }
00376   }
00377   cache_mgr.insert(m.navigation(), mod.navigation(), result.navigation());
00378   if (cv_empty) return result;
00379   else
00380     return cv.unite(result);
00381 }
00382 
00383 inline MonomialSet
00384 minimal_elements_cudd_style(MonomialSet m){
00385   return do_minimal_elements_cudd_style(m, MonomialSet(m.ring()));
00386 }
00387 
00388 
00389 
00390 
00391 
00392 inline MonomialSet
00393 minimal_elements_multiplied_monoms(MonomialSet m, Monomial lm){
00394   
00395     if (m.divisorsOf(lm).isZero()){
00396        m = m.existAbstract(lm);
00397        m = minimal_elements_cudd_style_unary(m);
00398 
00399        return m.unateProduct(lm.set());
00400     }
00401     return lm.set() ;
00402 }
00403 
00404 inline std::vector<Monomial>
00405 minimal_elements_multiplied(MonomialSet m, Monomial lm){ 
00406 
00407   MonomialSet result = minimal_elements_multiplied_monoms(m, lm);
00408   return std::vector<Monomial>(result.begin(), result.end());
00409 }
00410   
00411 
00412 
00413 //#define MIN_ELEMENTS_BINARY 1
00414 #ifdef MIN_ELEMENTS_BINARY
00415 inline MonomialSet
00416 minimal_elements_divided(MonomialSet m, Monomial lm, MonomialSet mod){
00417 
00418     if (m.divisorsOf(lm).isZero()){
00419         m=divide_monomial_divisors_out(m,lm);
00420         //mod=divide_monomial_divisors_out(mod,lm);
00421         return do_minimal_elements_cudd_style(m,mod);
00422     }
00423     return m.ring().one();
00424 }
00425 
00426 
00427 #else
00428 
00429 inline MonomialSet
00430 minimal_elements_divided(MonomialSet m, Monomial lm, MonomialSet mod){
00431 
00432     if (m.divisorsOf(lm).isZero()){
00433 
00434         m=m.existAbstract(lm);
00435         mod=mod.existAbstract(lm);
00436         //mod=divide_monomial_divisors_out(mod,lm);
00437         m=mod_mon_set(m,mod);
00438 
00439         return minimal_elements_cudd_style_unary(m);
00440     }
00441     return m.ring().one();
00442 }
00443 #endif
00444 
00445 inline void
00446 minimal_elements_divided(MonomialSet m, Monomial lm, MonomialSet mod,
00447                          std::vector<Exponent>& result){
00448   MonomialSet elements = minimal_elements_divided(m, lm, mod);
00449   result.resize(elements.length());
00450   std::copy(elements.expBegin(), elements.expEnd(), result.begin());
00451 }
00452 
00453 
00454 END_NAMESPACE_PBORIGB
00455 
00456 #endif /* polybori_groebner_minimal_elements_h_ */