1 /*++ 2 Copyright (c) 2019 Microsoft Corporation 3 4 Abstract: 5 6 Solver core based on pdd representation of polynomials 7 8 Author: 9 Nikolaj Bjorner (nbjorner) 10 Lev Nachmanson (levnach) 11 12 --*/ 13 14 #include "math/grobner/pdd_solver.h" 15 #include "math/grobner/pdd_simplifier.h" 16 #include "util/uint_set.h" 17 #include <math.h> 18 19 20 namespace dd { 21 22 /*** 23 A simple algorithm maintains two sets (S, A), 24 where S is m_processed, and A is m_to_simplify. 25 Initially S is empty and A contains the initial equations. 26 27 Each step proceeds as follows: 28 - pick a in A, and remove a from A 29 - simplify a using S 30 - simplify S using a 31 - for s in S: 32 b = superpose(a, s) 33 add b to A 34 - add a to S 35 - simplify A using a 36 37 Alternate: 38 39 - Fix a variable ordering x1 > x2 > x3 > .... 40 In each step: 41 - pick a in A with *highest* variable x_i in leading term of *lowest* degree. 42 - simplify a using S 43 - simplify S using a 44 - if a does not contains x_i, put it back into A and pick again (determine whether possible) 45 - for s in S: 46 b = superpose(a, s) 47 add b to A 48 - add a to S 49 - simplify A using a 50 51 - the variable ordering should be chosen from a candidate model M, 52 in a way that is compatible with weights that draw on the number of occurrences 53 in polynomials with violated evaluations and with the maximal slack (degree of freedom). 54 55 weight(x) := < count { p in P | x in p, M(p) != 0 }, min_{p in P | x in p} slack(p,x) > 56 57 slack is computed from interval assignments to variables, and is an interval in which x can possibly move 58 (an over-approximation). 59 60 */ 61 solver(reslimit & lim,pdd_manager & m)62 solver::solver(reslimit& lim, pdd_manager& m) : 63 m(m), 64 m_limit(lim), 65 m_conflict(nullptr) 66 {} 67 ~solver()68 solver::~solver() { 69 reset(); 70 } 71 72 adjust_cfg()73 void solver::adjust_cfg() { 74 auto & cfg = m_config; 75 IF_VERBOSE(3, verbose_stream() << "start saturate\n"; display_statistics(verbose_stream())); 76 cfg.m_eqs_threshold = static_cast<unsigned>(cfg.m_eqs_growth * ceil(log(1 + m_to_simplify.size()))* m_to_simplify.size()); 77 cfg.m_expr_size_limit = 0; 78 cfg.m_expr_degree_limit = 0; 79 for (equation* e: m_to_simplify) { 80 cfg.m_expr_size_limit = std::max(cfg.m_expr_size_limit, (unsigned)e->poly().tree_size()); 81 cfg.m_expr_degree_limit = std::max(cfg.m_expr_degree_limit, e->poly().degree()); 82 } 83 cfg.m_expr_size_limit *= cfg.m_expr_size_growth; 84 cfg.m_expr_degree_limit *= cfg.m_expr_degree_growth;; 85 86 IF_VERBOSE(3, verbose_stream() << "set m_config.m_eqs_threshold " << m_config.m_eqs_threshold << "\n"; 87 verbose_stream() << "set m_config.m_expr_size_limit to " << m_config.m_expr_size_limit << "\n"; 88 verbose_stream() << "set m_config.m_expr_degree_limit to " << m_config.m_expr_degree_limit << "\n"; 89 ); 90 91 } saturate()92 void solver::saturate() { 93 simplify(); 94 if (done()) { 95 return; 96 } 97 init_saturate(); 98 TRACE("dd.solver", display(tout);); 99 try { 100 while (!done() && step()) { 101 TRACE("dd.solver", display(tout);); 102 DEBUG_CODE(invariant();); 103 IF_VERBOSE(3, display_statistics(verbose_stream())); 104 } 105 DEBUG_CODE(invariant();); 106 } 107 catch (pdd_manager::mem_out) { 108 IF_VERBOSE(2, verbose_stream() << "mem-out\n"); 109 // don't reduce further 110 } 111 } 112 done()113 void solver::scoped_process::done() { 114 pdd p = e->poly(); 115 SASSERT(!p.is_val()); 116 if (p.degree() == 1) { 117 g.push_equation(solved, e); 118 } 119 else { 120 g.push_equation(processed, e); 121 } 122 e = nullptr; 123 } 124 ~scoped_process()125 solver::scoped_process::~scoped_process() { 126 if (e) { 127 pdd p = e->poly(); 128 SASSERT(!p.is_val()); 129 g.push_equation(processed, e); 130 } 131 } 132 simplify()133 void solver::simplify() { 134 simplifier s(*this); 135 s(); 136 } 137 138 superpose(equation const & eq)139 void solver::superpose(equation const & eq) { 140 for (equation* target : m_processed) { 141 superpose(eq, *target); 142 } 143 } 144 145 /* 146 Use a set of equations to simplify eq 147 */ simplify_using(equation & eq,equation_vector const & eqs)148 void solver::simplify_using(equation& eq, equation_vector const& eqs) { 149 bool simplified, changed_leading_term; 150 do { 151 simplified = false; 152 for (equation* p : eqs) { 153 if (try_simplify_using(eq, *p, changed_leading_term)) { 154 simplified = true; 155 } 156 if (canceled() || eq.poly().is_val()) { 157 break; 158 } 159 } 160 } 161 while (simplified && !eq.poly().is_val()); 162 163 TRACE("dd.solver", display(tout << "simplification result: ", eq);); 164 } 165 166 /* 167 Use the given equation to simplify equations in set 168 */ simplify_using(equation_vector & set,equation const & eq)169 void solver::simplify_using(equation_vector& set, equation const& eq) { 170 struct scoped_update { 171 equation_vector& set; 172 unsigned i, j, sz; 173 scoped_update(equation_vector& set): set(set), i(0), j(0), sz(set.size()) {} 174 void nextj() { 175 set[j] = set[i]; 176 set[i]->set_index(j++); 177 } 178 ~scoped_update() { 179 for (; i < sz; ++i) { 180 nextj(); 181 } 182 set.shrink(j); 183 } 184 }; 185 186 scoped_update sr(set); 187 for (; sr.i < sr.sz; ++sr.i) { 188 equation& target = *set[sr.i]; 189 bool changed_leading_term = false; 190 bool simplified = true; 191 simplified = !done() && try_simplify_using(target, eq, changed_leading_term); 192 193 if (simplified && is_trivial(target)) { 194 retire(&target); 195 } 196 else if (simplified && check_conflict(target)) { 197 // pushed to solved 198 } 199 else if (simplified && changed_leading_term) { 200 SASSERT(target.state() == processed); 201 push_equation(to_simplify, target); 202 if (!m_var2level.empty()) { 203 m_levelp1 = std::max(m_var2level[target.poly().var()]+1, m_levelp1); 204 } 205 } 206 else { 207 sr.nextj(); 208 } 209 } 210 } 211 212 /* 213 simplify target using source. 214 return true if the target was simplified. 215 set changed_leading_term if the target is in the m_processed set and the leading term changed. 216 */ try_simplify_using(equation & dst,equation const & src,bool & changed_leading_term)217 bool solver::try_simplify_using(equation& dst, equation const& src, bool& changed_leading_term) { 218 if (&src == &dst) { 219 return false; 220 } 221 m_stats.incr_simplified(); 222 pdd t = src.poly(); 223 pdd r = dst.poly().reduce(t); 224 if (r == dst.poly()){ 225 return false; 226 } 227 if (is_too_complex(r)) { 228 m_too_complex = true; 229 return false; 230 } 231 TRACE("dd.solver", 232 tout << "reduce: " << dst.poly() << "\n"; 233 tout << "using: " << t << "\n"; 234 tout << "to: " << r << "\n";); 235 changed_leading_term = dst.state() == processed && m.different_leading_term(r, dst.poly()); 236 dst = r; 237 dst = m_dep_manager.mk_join(dst.dep(), src.dep()); 238 update_stats_max_degree_and_size(dst); 239 return true; 240 } 241 simplify_using(equation & dst,equation const & src,bool & changed_leading_term)242 void solver::simplify_using(equation & dst, equation const& src, bool& changed_leading_term) { 243 if (&src == &dst) return; 244 m_stats.incr_simplified(); 245 pdd t = src.poly(); 246 pdd r = dst.poly().reduce(t); 247 changed_leading_term = dst.state() == processed && m.different_leading_term(r, dst.poly()); 248 TRACE("dd.solver", 249 tout << "reduce: " << dst.poly() << "\n"; 250 tout << "using: " << t << "\n"; 251 tout << "to: " << r << "\n";); 252 253 if (r == dst.poly()) return; 254 dst = r; 255 dst = m_dep_manager.mk_join(dst.dep(), src.dep()); 256 update_stats_max_degree_and_size(dst); 257 } 258 259 /* 260 let eq1: ab+q=0, and eq2: ac+e=0, then qc - eb = 0 261 */ superpose(equation const & eq1,equation const & eq2)262 void solver::superpose(equation const& eq1, equation const& eq2) { 263 TRACE("dd.solver_d", display(tout << "eq1=", eq1); display(tout << "eq2=", eq2);); 264 pdd r(m); 265 if (m.try_spoly(eq1.poly(), eq2.poly(), r) && !r.is_zero()) { 266 if (is_too_complex(r)) { 267 m_too_complex = true; 268 } 269 else { 270 m_stats.m_superposed++; 271 add(r, m_dep_manager.mk_join(eq1.dep(), eq2.dep())); 272 } 273 } 274 } 275 step()276 bool solver::step() { 277 m_stats.m_compute_steps++; 278 IF_VERBOSE(3, if (m_stats.m_compute_steps % 100 == 0) verbose_stream() << "compute steps = " << m_stats.m_compute_steps << "\n";); 279 equation* e = pick_next(); 280 if (!e) return false; 281 scoped_process sd(*this, e); 282 equation& eq = *e; 283 SASSERT(eq.state() == to_simplify); 284 simplify_using(eq, m_processed); 285 if (is_trivial(eq)) { sd.e = nullptr; retire(e); return true; } 286 if (check_conflict(eq)) { sd.e = nullptr; return false; } 287 m_too_complex = false; 288 simplify_using(m_processed, eq); 289 if (done()) return false; 290 TRACE("dd.solver", display(tout << "eq = ", eq);); 291 superpose(eq); 292 simplify_using(m_to_simplify, eq); 293 if (done()) return false; 294 if (!m_too_complex) sd.done(); 295 return true; 296 } 297 init_saturate()298 void solver::init_saturate() { 299 unsigned_vector const& l2v = m.get_level2var(); 300 m_level2var.resize(l2v.size()); 301 m_var2level.resize(l2v.size()); 302 for (unsigned i = 0; i < l2v.size(); ++i) { 303 m_level2var[i] = l2v[i]; 304 m_var2level[l2v[i]] = i; 305 } 306 m_levelp1 = m_level2var.size(); 307 } 308 pick_next()309 solver::equation* solver::pick_next() { 310 while (m_levelp1 > 0) { 311 unsigned v = m_level2var[m_levelp1-1]; 312 equation* eq = nullptr; 313 for (equation* curr : m_to_simplify) { 314 SASSERT(curr->idx() != UINT_MAX); 315 pdd const& p = curr->poly(); 316 if (curr->state() == to_simplify && p.var() == v) { 317 if (!eq || is_simpler(*curr, *eq)) 318 eq = curr; 319 } 320 } 321 if (eq) { 322 pop_equation(eq); 323 return eq; 324 } 325 --m_levelp1; 326 } 327 return nullptr; 328 } 329 equations()330 solver::equation_vector const& solver::equations() { 331 m_all_eqs.reset(); 332 for (equation* eq : m_solved) m_all_eqs.push_back(eq); 333 for (equation* eq : m_to_simplify) m_all_eqs.push_back(eq); 334 for (equation* eq : m_processed) m_all_eqs.push_back(eq); 335 return m_all_eqs; 336 } 337 reset()338 void solver::reset() { 339 for (equation* e : m_solved) dealloc(e); 340 for (equation* e : m_to_simplify) dealloc(e); 341 for (equation* e : m_processed) dealloc(e); 342 m_solved.reset(); 343 m_processed.reset(); 344 m_to_simplify.reset(); 345 m_stats.reset(); 346 m_level2var.reset(); 347 m_var2level.reset(); 348 m_conflict = nullptr; 349 } 350 add(pdd const & p,u_dependency * dep)351 void solver::add(pdd const& p, u_dependency * dep) { 352 if (p.is_zero()) return; 353 equation * eq = alloc(equation, p, dep); 354 if (check_conflict(*eq)) { 355 return; 356 } 357 push_equation(to_simplify, eq); 358 359 if (!m_var2level.empty()) { 360 m_levelp1 = std::max(m_var2level[p.var()]+1, m_levelp1); 361 } 362 update_stats_max_degree_and_size(*eq); 363 } 364 canceled()365 bool solver::canceled() { 366 return m_limit.is_canceled(); 367 } 368 done()369 bool solver::done() { 370 return 371 m_to_simplify.size() + m_processed.size() >= m_config.m_eqs_threshold || 372 m_stats.simplified() >= m_config.m_max_simplified || 373 canceled() || 374 m_stats.m_compute_steps > m_config.m_max_steps || 375 m_conflict != nullptr; 376 } 377 get_queue(equation const & eq)378 solver::equation_vector& solver::get_queue(equation const& eq) { 379 switch (eq.state()) { 380 case processed: return m_processed; 381 case to_simplify: return m_to_simplify; 382 case solved: return m_solved; 383 } 384 UNREACHABLE(); 385 return m_to_simplify; 386 } 387 del_equation(equation * eq)388 void solver::del_equation(equation* eq) { 389 pop_equation(eq); 390 retire(eq); 391 } 392 retire(equation * eq)393 void solver::retire(equation* eq) { 394 #if 0 395 // way to check if retired equations are ever accessed. 396 eq->set_index(UINT_MAX); 397 #else 398 dealloc(eq); 399 #endif 400 } 401 402 pop_equation(equation & eq)403 void solver::pop_equation(equation& eq) { 404 equation_vector& v = get_queue(eq); 405 unsigned idx = eq.idx(); 406 if (idx != v.size() - 1) { 407 equation* eq2 = v.back(); 408 eq2->set_index(idx); 409 v[idx] = eq2; 410 } 411 v.pop_back(); 412 } 413 push_equation(eq_state st,equation & eq)414 void solver::push_equation(eq_state st, equation& eq) { 415 SASSERT(st != to_simplify || !eq.poly().is_val()); 416 SASSERT(st != processed || !eq.poly().is_val()); 417 eq.set_state(st); 418 equation_vector& v = get_queue(eq); 419 eq.set_index(v.size()); 420 v.push_back(&eq); 421 } 422 update_stats_max_degree_and_size(const equation & e)423 void solver::update_stats_max_degree_and_size(const equation& e) { 424 m_stats.m_max_expr_size = std::max(m_stats.m_max_expr_size, e.poly().tree_size()); 425 m_stats.m_max_expr_degree = std::max(m_stats.m_max_expr_degree, e.poly().degree()); 426 } 427 collect_statistics(statistics & st) const428 void solver::collect_statistics(statistics& st) const { 429 st.update("dd.solver.steps", m_stats.m_compute_steps); 430 st.update("dd.solver.simplified", m_stats.simplified()); 431 st.update("dd.solver.superposed", m_stats.m_superposed); 432 st.update("dd.solver.processed", m_processed.size()); 433 st.update("dd.solver.solved", m_solved.size()); 434 st.update("dd.solver.to_simplify", m_to_simplify.size()); 435 st.update("dd.solver.degree", m_stats.m_max_expr_degree); 436 st.update("dd.solver.size", m_stats.m_max_expr_size); 437 } 438 display(std::ostream & out,const equation & eq) const439 std::ostream& solver::display(std::ostream & out, const equation & eq) const { 440 out << eq.poly() << "\n"; 441 if (m_print_dep) m_print_dep(eq.dep(), out); 442 return out; 443 } 444 display(std::ostream & out) const445 std::ostream& solver::display(std::ostream& out) const { 446 out << "solved\n"; for (auto e : m_solved) display(out, *e); 447 out << "processed\n"; for (auto e : m_processed) display(out, *e); 448 out << "to_simplify\n"; for (auto e : m_to_simplify) display(out, *e); 449 return display_statistics(out); 450 } 451 display_statistics(std::ostream & out) const452 std::ostream& solver::display_statistics(std::ostream& out) const { 453 statistics st; 454 collect_statistics(st); 455 st.display(out); 456 out << "\n----\n"; 457 return out; 458 } 459 invariant() const460 void solver::invariant() const { 461 return; // disabling the check that seem incorrect after changing in the order of pdd expressions 462 // equations in processed have correct indices 463 // they are labled as processed 464 unsigned i = 0; 465 for (auto* e : m_processed) { 466 VERIFY(e->state() == processed); 467 VERIFY(e->idx() == i); 468 VERIFY(!e->poly().is_val()); 469 ++i; 470 } 471 472 i = 0; 473 uint_set head_vars; 474 for (auto* e : m_solved) { 475 VERIFY(e->state() == solved); 476 VERIFY(e->idx() == i); 477 ++i; 478 pdd p = e->poly(); 479 if (p.degree() == 1) { 480 unsigned v = p.var(); 481 SASSERT(!head_vars.contains(v)); 482 head_vars.insert(v); 483 } 484 } 485 486 if (!head_vars.empty()) { 487 for (auto * e : m_to_simplify) { 488 for (auto v : m.free_vars(e->poly())) VERIFY(!head_vars.contains(v)); 489 } 490 for (auto * e : m_processed) { 491 for (auto v : m.free_vars(e->poly())) VERIFY(!head_vars.contains(v)); 492 } 493 } 494 495 // equations in to_simplify have correct indices 496 // they are labeled as non-processed 497 i = 0; 498 for (auto* e : m_to_simplify) { 499 VERIFY(e->idx() == i); 500 VERIFY(e->state() == to_simplify); 501 pdd const& p = e->poly(); 502 VERIFY(!p.is_val()); 503 ++i; 504 } 505 } 506 } 507 508