PolyBoRi
|
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_ */