From de6a0ab1a7991e8bc713621320b319131f57e1a6 Mon Sep 17 00:00:00 2001 From: Jakob Rath Date: Mon, 1 Aug 2022 14:49:06 +0200 Subject: [PATCH] PDD operations --- src/math/dd/dd_pdd.cpp | 511 ++++++++++++++++++++++++++++++++++++++--- src/math/dd/dd_pdd.h | 123 ++++++++-- src/test/pdd.cpp | 301 +++++++++++++++++++++++- 3 files changed, 877 insertions(+), 58 deletions(-) diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index b8946ebf6..cca2d5e04 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -34,6 +34,7 @@ namespace dd { s = mod2_e; m_semantics = s; m_mod2N = rational::power_of_two(power_of_2); + m_max_value = m_mod2N - 1; m_power_of_2 = power_of_2; unsigned_vector l2v; for (unsigned i = 0; i < num_vars; ++i) l2v.push_back(i); @@ -50,6 +51,7 @@ namespace dd { void pdd_manager::reset(unsigned_vector const& level2var) { reset_op_cache(); + m_factor_cache.reset(); m_node_table.reset(); m_nodes.reset(); m_free_nodes.reset(); @@ -126,7 +128,7 @@ namespace dd { * Example: 2^4*x + 2 is non-zero for every x. */ - bool pdd_manager::is_non_zero(PDD p) { + bool pdd_manager::is_never_zero(PDD p) { if (is_val(p)) return !is_zero(p); if (m_semantics != mod2N_e) @@ -163,7 +165,11 @@ namespace dd { return true; } - pdd pdd_manager::subst_val(pdd const& p, vector> const& _s) { + pdd pdd_manager::subst_val(pdd const& p, pdd const& s) { + return pdd(apply(p.root, s.root, pdd_subst_val_op), this); + } + + pdd pdd_manager::subst_val0(pdd const& p, vector> const& _s) { typedef std::pair pr; vector s(_s); std::function compare_level = @@ -171,10 +177,16 @@ namespace dd { std::sort(s.begin(), s.end(), compare_level); pdd r(one()); for (auto const& q : s) - r = (r*mk_var(q.first)) + q.second; - return pdd(apply(p.root, r.root, pdd_subst_val_op), this); + r = (r * mk_var(q.first)) + q.second; + return subst_val(p, r); } + pdd pdd_manager::subst_add(pdd const& s, unsigned v, rational const& val) { + pdd v_val = mk_var(v) + val; + return pdd(apply(s.root, v_val.root, pdd_subst_add_op), this); + } + + pdd_manager::PDD pdd_manager::apply(PDD arg1, PDD arg2, pdd_op op) { bool first = true; SASSERT(well_formed()); @@ -216,6 +228,7 @@ namespace dd { if (is_val(p) && is_val(q)) return imk_val(val(p) - val(q)); if (m_semantics != mod2_e) break; op = pdd_add_op; + Z3_fallthrough; case pdd_add_op: if (is_zero(p)) return q; if (is_zero(q)) return p; @@ -239,12 +252,16 @@ namespace dd { break; case pdd_subst_val_op: while (!is_val(q) && !is_val(p)) { - if (level(p) == level(q)) break; - if (level(p) < level(q)) q = lo(q); - else p = lo(p); + if (level(p) < level(q)) q = hi(q); + else break; } if (is_val(p) || is_val(q)) return p; break; + case pdd_subst_add_op: + if (is_one(p)) return q; + SASSERT(!is_val(p)); + SASSERT(!is_val(q)); + break; default: UNREACHABLE(); break; @@ -345,7 +362,8 @@ namespace dd { push(make_node(level_p, lo(n), read(1))); r = make_node(level_p, bd, read(1)); npop = 7; - } else { + } + else { push(make_node(level_p, n, ac)); r = make_node(level_p, bd, read(1)); npop = 6; @@ -387,12 +405,33 @@ namespace dd { case pdd_subst_val_op: SASSERT(!is_val(p)); SASSERT(!is_val(q)); - SASSERT(level_p == level_q); + SASSERT(level_p >= level_q); push(apply_rec(lo(p), q, pdd_subst_val_op)); // lo := subst(lo(p), s) push(apply_rec(hi(p), q, pdd_subst_val_op)); // hi := subst(hi(p), s) - push(apply_rec(lo(q), read(1), pdd_mul_op)); // hi := hi*s[var(p)] - r = apply_rec(read(1), read(3), pdd_add_op); // r := hi + lo := subst(lo(p),s) + s[var(p)]*subst(hi(p),s) - npop = 3; + + if (level_p > level_q) { + r = make_node(level_p, read(2), read(1)); + npop = 2; + } + else { + push(apply_rec(lo(q), read(1), pdd_mul_op)); // hi := hi*s[var(p)] + r = apply_rec(read(1), read(3), pdd_add_op); // r := hi + lo := subst(lo(p),s) + s[var(p)]*subst(hi(p),s) + npop = 3; + } + break; + case pdd_subst_add_op: + SASSERT(!is_val(p)); + SASSERT(!is_val(q)); + SASSERT(level_p != level_q); + if (level_p < level_q) { + r = make_node(level_q, lo(q), p); // p*hi(q) + lo(q) + npop = 0; + } + else { + push(apply_rec(hi(p), q, pdd_subst_add_op)); // hi := add_subst(hi(p), q) + r = make_node(level_p, lo(p), read(1)); // r := hi*var(p) + lo(p) + npop = 1; + } break; default: r = null_pdd; @@ -405,6 +444,7 @@ namespace dd { return r; } + pdd pdd_manager::minus(pdd const& a) { if (m_semantics == mod2_e) { return a; @@ -442,6 +482,113 @@ namespace dd { return r; } + /** + * Divide PDD by a constant value. + * + * IMPORTANT: Performs regular numerical division. + * For semantics 'mod2N_e', this means that 'c' must be an integer + * and all coefficients of 'a' must be divisible by 'c'. + * + * NOTE: Why do we not just use 'mul(a, inv(c))' instead? + * In case of semantics 'mod2N_e', an invariant is that all PDDs have integer coefficients. + * But such a multiplication would create nodes with non-integral coefficients. + */ + pdd pdd_manager::div(pdd const& a, rational const& c) { + pdd res(zero_pdd, this); + VERIFY(try_div(a, c, res)); + return res; + } + + bool pdd_manager::try_div(pdd const& a, rational const& c, pdd& out_result) { + if (m_semantics == free_e) { + // Don't cache separately for the free semantics; + // use 'mul' so we can share results for a/c and a*(1/c). + out_result = mul(inv(c), a); + return true; + } + SASSERT(c.is_int()); + bool first = true; + SASSERT(well_formed()); + scoped_push _sp(*this); + while (true) { + try { + PDD res = div_rec(a.root, c, null_pdd); + if (res != null_pdd) + out_result = pdd(res, this); + SASSERT(well_formed()); + return res != null_pdd; + } + catch (const mem_out &) { + try_gc(); + if (!first) throw; + first = false; + } + } + } + + /// Returns null_pdd if one of the coefficients is not divisible by c. + pdd_manager::PDD pdd_manager::div_rec(PDD a, rational const& c, PDD c_pdd) { + SASSERT(m_semantics != free_e); + SASSERT(c.is_int()); + if (is_zero(a)) + return zero_pdd; + if (is_val(a)) { + rational r = val(a) / c; + if (r.is_int()) + return imk_val(r); + else + return null_pdd; + } + if (c_pdd == null_pdd) + c_pdd = imk_val(c); + op_entry* e1 = pop_entry(a, c_pdd, pdd_div_const_op); + op_entry const* e2 = m_op_cache.insert_if_not_there(e1); + if (check_result(e1, e2, a, c_pdd, pdd_div_const_op)) + return e2->m_result; + push(div_rec(lo(a), c, c_pdd)); + push(div_rec(hi(a), c, c_pdd)); + PDD l = read(2); + PDD h = read(1); + PDD res = null_pdd; + if (l != null_pdd && h != null_pdd) + res = make_node(level(a), l, h); + pop(2); + e1->m_result = res; + return res; + } + + pdd pdd_manager::pow(pdd const &p, unsigned j) { + return pdd(pow(p.root, j), this); + } + + pdd_manager::PDD pdd_manager::pow(PDD p, unsigned j) { + if (j == 0) + return one_pdd; + else if (j == 1) + return p; + else if (is_zero(p)) + return zero_pdd; + else if (is_one(p)) + return one_pdd; + else if (is_val(p)) + return imk_val(power(val(p), j)); + else + return pow_rec(p, j); + } + + pdd_manager::PDD pdd_manager::pow_rec(PDD p, unsigned j) { + SASSERT(j > 0); + if (j == 1) + return p; + // j even: pow(p,2*j') = pow(p*p,j') + // j odd: pow(p,2*j'+1) = p*pow(p*p,j') + PDD q = pow_rec(apply(p, p, pdd_mul_op), j / 2); + if (j & 1) { + q = apply(q, p, pdd_mul_op); + } + return q; + } + // // produce polynomial where a is reduced by b. // all monomials in a that are divisible by lm(b) @@ -691,27 +838,37 @@ namespace dd { * factor p into lc*v^degree + rest * such that degree(rest, v) < degree * Initial implementation is very naive - * - memoize intermediary results */ void pdd_manager::factor(pdd const& p, unsigned v, unsigned degree, pdd& lc, pdd& rest) { unsigned level_v = m_var2level[v]; if (degree == 0) { - lc = p; - rest = zero(); - } - else if (level(p.root) < level_v) { lc = zero(); rest = p; + return; } - else if (level(p.root) > level_v) { + if (level(p.root) < level_v) { + lc = zero(); + rest = p; + return; + } + // Memoize nontrivial cases + auto* et = m_factor_cache.insert_if_not_there2({p.root, v, degree}); + factor_entry* e = &et->get_data(); + if (e->is_valid()) { + lc = pdd(e->m_lc, this); + rest = pdd(e->m_rest, this); + return; + } + + if (level(p.root) > level_v) { pdd lc1 = zero(), rest1 = zero(); pdd vv = mk_var(p.var()); factor(p.hi(), v, degree, lc, rest); factor(p.lo(), v, degree, lc1, rest1); - lc += lc1; - rest += rest1; lc *= vv; rest *= vv; + lc += lc1; + rest += rest1; } else { unsigned d = 0; @@ -733,8 +890,228 @@ namespace dd { rest = p; } } + et = m_factor_cache.insert_if_not_there2({p.root, v, degree}); + e = &et->get_data(); + e->m_lc = lc.root; + e->m_rest = rest.root; } + bool pdd_manager::factor(pdd const& p, unsigned v, unsigned degree, pdd& lc) { + pdd rest = lc; + factor(p, v, degree, lc, rest); + return rest.is_zero(); + } + + /** + * Apply function f to all coefficients of the polynomial. + * The function should be of type + * rational const& -> rational + * rational const& -> unsigned + * and should always return integers. + * + * NOTE: the operation is not cached. + */ + template + pdd pdd_manager::map_coefficients(pdd const& p, Fn f) { + if (p.is_val()) { + return mk_val(rational(f(p.val()))); + } else { + pdd x = mk_var(p.var()); + pdd lo = map_coefficients(p.lo(), f); + pdd hi = map_coefficients(p.hi(), f); + return x*hi + lo; + } + } + + /** + * Perform S-polynomial reduction on p by q, + * treating monomial with v as leading. + * + * p = a v^l + b = a' 2^j v^l + b + * q = c v^m + d = c' 2^j v^m + d + * such that + * deg(v, p) = l, i.e., v does not divide a and there is no v^l in b + * deg(v, q) = m, i.e., v does not divide c and there is no v^m in d + * l >= m + * j maximal, i.e., not both of a', c' are divisible by 2 + * + * Then we reduce p by q: + * + * r = c' p - a' v^(l-m) q + * = b c' - a' d v^(l-m) + */ + bool pdd_manager::resolve(unsigned v, pdd const& p, pdd const& q, pdd& r) { + unsigned const l = p.degree(v); + unsigned const m = q.degree(v); + // no reduction + if (l < m || m == 0) + return false; + + pdd a = zero(); + pdd b = zero(); + pdd c = zero(); + pdd d = zero(); + p.factor(v, l, a, b); + q.factor(v, m, c, d); + unsigned const j = std::min(max_pow2_divisor(a), max_pow2_divisor(c)); + SASSERT(j != UINT_MAX); // should only be possible when both l and m are 0 + rational const pow2j = rational::power_of_two(j); + pdd const aa = div(a, pow2j); + pdd const cc = div(c, pow2j); + pdd vv = pow(mk_var(v), l - m); + r = b * cc - aa * d * vv; + return true; + } + + /** + * Reduce polynomial a with respect to b by eliminating terms using v + * + * a := a1*v^l + a2 + * b := b1*v^m + b2 + * l >= m + * q, r := quot_rem(a1, b1) + * that is: + * q * b1 + r = a1 + * r = 0 + * result := reduce(v, q*b2*v^{l-m}, b) + reduce(v, a2, b) + */ + pdd pdd_manager::reduce(unsigned v, pdd const& a, pdd const& b) { + unsigned const m = b.degree(v); + // no reduction + if (m == 0) + return a; + + pdd b1 = zero(); + pdd b2 = zero(); + b.factor(v, m, b1, b2); + + // TODO - generalize this case to when leading coefficient is not a value + if (m_semantics == mod2N_e && b1.is_val() && b1.val().is_odd() && !b1.is_one()) { + rational b_inv; + VERIFY(b1.val().mult_inverse(m_power_of_2, b_inv)); + b1 = 1; + b2 *= b_inv; + } + + return reduce(v, a, m, b1, b2); + } + + pdd pdd_manager::reduce(unsigned v, pdd const& a, unsigned m, pdd const& b1, pdd const& b2) { + SASSERT(m > 0); + unsigned const l = a.degree(v); + if (l < m) + return a; + + pdd a1 = zero(); + pdd a2 = zero(); + pdd q = zero(); + pdd r = zero(); + a.factor(v, l, a1, a2); + + quot_rem(a1, b1, q, r); + if (r.is_zero()) { + SASSERT(q * b1 == a1); + a1 = -q * b2; + if (l > m) + a1 = reduce(v, a1 * pow(mk_var(v), l - m), m, b1, b2); + } + else + a1 = a1 * pow(mk_var(v), l); + a2 = reduce(v, a2, m, b1, b2); + + return a1 + a2; + } + + /** + * quotient/remainder of 'a' divided by 'b' + * a := x*hi + lo + * x > level(b): + * hi = q1*b + r1 + * lo = q2*b + r2 + * x*hi + lo = (x*q1 + q2)*b + (x*r1 + r2) + * q := x*q1 + q2 + * r := x*r1 + r2 + * Some special cases. + * General multi-variable polynomial division is TBD. + */ + void pdd_manager::quot_rem(pdd const& a, pdd const& b, pdd& q, pdd& r) { + if (level(a.root) > level(b.root)) { + pdd q1(*this), q2(*this), r1(*this), r2(*this); + quot_rem(a.hi(), b, q1, r1); + quot_rem(a.lo(), b, q2, r2); + q = mk_var(a.var()) * q1 + q2; + r = mk_var(a.var()) * r1 + r2; + } + else if (level(a.root) < level(b.root)) { + q = zero(); + r = a; + } + else if (a == b) { + q = one(); + r = zero(); + } + else if (a.is_val() && b.is_val() && divides(b.val(), a.val())) { + q = mk_val(a.val() / b.val()); + r = zero(); + } + else if (a.is_val() || b.is_val()) { + q = zero(); + r = a; + } + else { + SASSERT(level(a.root) == level(b.root)); + pdd q1(*this), q2(*this), r1(*this), r2(*this); + quot_rem(a.hi(), b.hi(), q1, r1); + quot_rem(a.lo(), b.lo(), q2, r2); + if (q1 == q2 && r1.is_zero() && r2.is_zero()) { + q = q1; + r = zero(); + } + else { + q = zero(); + r = a; + } + } + } + + /** + * Returns the largest j such that 2^j divides p. + */ + unsigned pdd_manager::max_pow2_divisor(PDD p) { + init_mark(); + unsigned min_j = UINT_MAX; + m_todo.push_back(p); + while (!m_todo.empty()) { + PDD r = m_todo.back(); + m_todo.pop_back(); + if (is_marked(r)) { + continue; + } + set_mark(r); + if (is_zero(r)) { + // skip + } + else if (is_val(r)) { + rational const& c = val(r); + if (c.is_odd()) { + m_todo.reset(); + return 0; + } else { + unsigned j = c.trailing_zeros(); + min_j = std::min(j, min_j); + } + } + else { + m_todo.push_back(lo(r)); + m_todo.push_back(hi(r)); + } + } + return min_j; + } + + unsigned pdd_manager::max_pow2_divisor(pdd const& p) { + return max_pow2_divisor(p.root); + } bool pdd_manager::is_linear(pdd const& p) { return is_linear(p.root); @@ -764,6 +1141,32 @@ namespace dd { } } + /** Determine whether p contains at most one variable. */ + bool pdd_manager::is_univariate(PDD p) { + unsigned const lvl = level(p); + while (!is_val(p)) { + if (!is_val(lo(p))) + return false; + if (level(p) != lvl) + return false; + p = hi(p); + } + return true; + } + + /** + * Push coefficients of univariate polynomial in order of ascending degree. + * Example: a*x^2 + b*x + c ==> [ c, b, a ] + */ + void pdd_manager::get_univariate_coefficients(PDD p, vector& coeff) { + SASSERT(is_univariate(p)); + while (!is_val(p)) { + coeff.push_back(val(lo(p))); + p = hi(p); + } + coeff.push_back(val(p)); + } + /* \brief determine if v occurs as a leaf variable. */ @@ -832,9 +1235,8 @@ namespace dd { if (m_semantics == mod2N_e && (r < 0 || r >= m_mod2N)) return imk_val(mod(r, m_mod2N)); const_info info; - if (!m_mpq_table.find(r, info)) { + if (!m_mpq_table.find(r, info)) init_value(info, r); - } return info.m_node_index; } @@ -1005,7 +1407,7 @@ namespace dd { unsigned pdd_manager::degree(PDD p, unsigned v) { init_mark(); - unsigned level_v = level(v); + unsigned level_v = m_var2level[v]; unsigned max_d = 0, d = 0; m_todo.push_back(p); while (!m_todo.empty()) { @@ -1017,15 +1419,17 @@ namespace dd { else if (level(r) < level_v) m_todo.pop_back(); else if (level(r) == level_v) { - d = 1; - while (!is_val(hi(r)) && level(hi(r)) == level_v) { + d = 0; + do { ++d; + set_mark(r); r = hi(r); - } + } while (!is_val(r) && level(r) == level_v); max_d = std::max(d, max_d); m_todo.pop_back(); } else { + set_mark(r); m_todo.push_back(lo(r)); m_todo.push_back(hi(r)); } @@ -1128,6 +1532,7 @@ namespace dd { } void pdd_manager::gc() { + m_gc_generation++; init_dmark(); m_free_nodes.reset(); SASSERT(well_formed()); @@ -1167,6 +1572,8 @@ namespace dd { m_op_cache.insert(e); } + m_factor_cache.reset(); + m_node_table.reset(); // re-populate node cache for (unsigned i = m_nodes.size(); i-- > 2; ) { @@ -1221,14 +1628,32 @@ namespace dd { rational c = abs(m.first); m.second.reverse(); if (!c.is_one() || m.second.empty()) { - out << c; + if (m_semantics == mod2N_e && mod(-c, m_mod2N) < c) + out << -mod(-c, m_mod2N); + else + out << c; if (!m.second.empty()) out << "*"; } - bool f = true; + unsigned v_prev = UINT_MAX; + unsigned pow = 0; for (unsigned v : m.second) { - if (!f) out << "*"; - f = false; - out << "v" << v; + if (v == v_prev) { + pow++; + continue; + } + if (v_prev != UINT_MAX) { + out << "v" << v_prev; + if (pow > 1) + out << "^" << pow; + out << "*"; + } + pow = 1; + v_prev = v; + } + if (v_prev != UINT_MAX) { + out << "v" << v_prev; + if (pow > 1) + out << "^" << pow; } } if (first) out << "0"; @@ -1291,6 +1716,30 @@ namespace dd { return *this; } + pdd& pdd::operator=(unsigned k) { + m.dec_ref(root); + root = m.mk_val(k).root; + m.inc_ref(root); + return *this; + } + + pdd& pdd::operator=(rational const& k) { + m.dec_ref(root); + root = m.mk_val(k).root; + m.inc_ref(root); + return *this; + } + + rational const& pdd::leading_coefficient() const { + pdd p = *this; + while (!p.is_val()) + p = p.hi(); + return p.val(); + } + + pdd pdd::shl(unsigned n) const { + return (*this) * rational::power_of_two(n); + } /** * \brief substitute variable v by r. diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index 9cd454698..aef0eb6f7 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -23,6 +23,11 @@ Abstract: - try_spoly(a, b, c) returns true if lt(a) and lt(b) have a non-trivial overlap. c is the resolvent (S polynomial). + - reduce(v, a, b) reduces 'a' using b = 0 with respect to eliminating v. + Given b = v^l*b1 + b2, where l is the leading coefficient of v in b + Given a = v^m*a1 + b2, where m >= l is the leading coefficient of v in b. + reduce(a1, b1)*v^{m - l} + reduce(v, a2, b); + Author: Nikolaj Bjorner (nbjorner) 2019-12-17 @@ -63,7 +68,9 @@ namespace dd { pdd_mul_op = 5, pdd_reduce_op = 6, pdd_subst_val_op = 7, - pdd_no_op = 8 + pdd_subst_add_op = 8, + pdd_div_const_op = 9, + pdd_no_op = 10 }; struct node { @@ -140,9 +147,44 @@ namespace dd { typedef ptr_hashtable op_table; + struct factor_entry { + factor_entry(PDD p, unsigned v, unsigned degree): + m_p(p), + m_v(v), + m_degree(degree), + m_lc(UINT_MAX), + m_rest(UINT_MAX) + {} + + factor_entry(): m_p(0), m_v(0), m_degree(0), m_lc(UINT_MAX), m_rest(UINT_MAX) {} + + PDD m_p; // input + unsigned m_v; // input + unsigned m_degree; // input + PDD m_lc; // output + PDD m_rest; // output + + bool is_valid() { return m_lc != UINT_MAX; } + + unsigned hash() const { return mk_mix(m_p, m_v, m_degree); } + }; + + struct hash_factor_entry { + unsigned operator()(factor_entry const& e) const { return e.hash(); } + }; + + struct eq_factor_entry { + bool operator()(factor_entry const& a, factor_entry const& b) const { + return a.m_p == b.m_p && a.m_v == b.m_v && a.m_degree == b.m_degree; + } + }; + + typedef hashtable factor_table; + svector m_nodes; vector m_values; op_table m_op_cache; + factor_table m_factor_cache; node_table m_node_table; mpq_table m_mpq_table; svector m_pdd_stack; @@ -164,7 +206,9 @@ namespace dd { unsigned_vector m_free_values; rational m_freeze_value; rational m_mod2N; - unsigned m_power_of_2 { 0 }; + unsigned m_power_of_2 = 0; + rational m_max_value; + unsigned m_gc_generation = 0; ///< will be incremented on each GC void reset_op_cache(); void init_nodes(unsigned_vector const& l2v); @@ -177,6 +221,9 @@ namespace dd { PDD apply(PDD arg1, PDD arg2, pdd_op op); PDD apply_rec(PDD arg1, PDD arg2, pdd_op op); PDD minus_rec(PDD p); + PDD div_rec(PDD p, rational const& c, PDD c_pdd); + PDD pow(PDD p, unsigned j); + PDD pow_rec(PDD p, unsigned j); PDD reduce_on_match(PDD a, PDD b); bool lm_occurs(PDD p, PDD q) const; @@ -206,7 +253,8 @@ namespace dd { inline bool is_one(PDD p) const { return p == one_pdd; } inline bool is_val(PDD p) const { return m_nodes[p].is_val(); } inline bool is_internal(PDD p) const { return m_nodes[p].is_internal(); } - bool is_non_zero(PDD p); + inline bool is_var(PDD p) const { return !is_val(p) && is_zero(lo(p)) && is_one(hi(p)); } + bool is_never_zero(PDD p); inline unsigned level(PDD p) const { return m_nodes[p].m_level; } inline unsigned var(PDD p) const { return m_level2var[level(p)]; } inline PDD lo(PDD p) const { return m_nodes[p].m_lo; } @@ -226,8 +274,14 @@ namespace dd { unsigned degree(pdd const& p) const; unsigned degree(PDD p) const; unsigned degree(PDD p, unsigned v); + unsigned max_pow2_divisor(PDD p); + unsigned max_pow2_divisor(pdd const& p); + template pdd map_coefficients(pdd const& p, Fn f); void factor(pdd const& p, unsigned v, unsigned degree, pdd& lc, pdd& rest); + bool factor(pdd const& p, unsigned v, unsigned degree, pdd& lc); + + pdd reduce(unsigned v, pdd const& a, unsigned m, pdd const& b1, pdd const& b2); bool var_is_leaf(PDD p, unsigned v); @@ -258,7 +312,7 @@ namespace dd { struct mem_out {}; - pdd_manager(unsigned nodes, semantics s = free_e, unsigned power_of_2 = 0); + pdd_manager(unsigned num_vars, semantics s = free_e, unsigned power_of_2 = 0); ~pdd_manager(); semantics get_semantics() const { return m_semantics; } @@ -278,13 +332,21 @@ namespace dd { pdd sub(pdd const& a, pdd const& b); pdd mul(pdd const& a, pdd const& b); pdd mul(rational const& c, pdd const& b); + pdd div(pdd const& a, rational const& c); + bool try_div(pdd const& a, rational const& c, pdd& out_result); pdd mk_or(pdd const& p, pdd const& q); pdd mk_xor(pdd const& p, pdd const& q); pdd mk_xor(pdd const& p, unsigned q); pdd mk_not(pdd const& p); pdd reduce(pdd const& a, pdd const& b); - pdd subst_val(pdd const& a, vector> const& s); + pdd subst_val0(pdd const& a, vector> const& s); pdd subst_val(pdd const& a, unsigned v, rational const& val); + pdd subst_val(pdd const& a, pdd const& s); + pdd subst_add(pdd const& s, unsigned v, rational const& val); + bool resolve(unsigned v, pdd const& p, pdd const& q, pdd& r); + pdd reduce(unsigned v, pdd const& a, pdd const& b); + void quot_rem(pdd const& a, pdd const& b, pdd& q, pdd& r); + pdd pow(pdd const& p, unsigned j); bool is_linear(PDD p) { return degree(p) == 1; } bool is_linear(pdd const& p); @@ -294,6 +356,9 @@ namespace dd { bool is_monomial(PDD p); + bool is_univariate(PDD p); + void get_univariate_coefficients(PDD p, vector& coeff); + // create an spoly r if leading monomials of a and b overlap bool try_spoly(pdd const& a, pdd const& b, pdd& r); @@ -308,6 +373,8 @@ namespace dd { unsigned num_vars() const { return m_var2pdd.size(); } unsigned power_of_2() const { return m_power_of_2; } + rational const& max_value() const { return m_max_value; } + rational const& two_to_N() const { return m_mod2N; } unsigned_vector const& free_vars(pdd const& p); @@ -330,21 +397,30 @@ namespace dd { pdd(pdd const& other): root(other.root), m(other.m) { m.inc_ref(root); } pdd(pdd && other) noexcept : root(0), m(other.m) { std::swap(root, other.root); } pdd& operator=(pdd const& other); + pdd& operator=(unsigned k); + pdd& operator=(rational const& k); ~pdd() { m.dec_ref(root); } pdd lo() const { return pdd(m.lo(root), m); } pdd hi() const { return pdd(m.hi(root), m); } unsigned index() const { return root; } unsigned var() const { return m.var(root); } rational const& val() const { SASSERT(is_val()); return m.val(root); } + rational const& leading_coefficient() const; bool is_val() const { return m.is_val(root); } bool is_one() const { return m.is_one(root); } bool is_zero() const { return m.is_zero(root); } bool is_linear() const { return m.is_linear(root); } + bool is_var() const { return m.is_var(root); } + /** Polynomial is of the form a * x + b for numerals a, b. */ + bool is_unilinear() const { return !is_val() && lo().is_val() && hi().is_val(); } bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); } bool is_offset() const { return !is_val() && lo().is_val() && hi().is_one(); } bool is_binary() const { return m.is_binary(root); } bool is_monomial() const { return m.is_monomial(root); } - bool is_non_zero() const { return m.is_non_zero(root); } + bool is_univariate() const { return m.is_univariate(root); } + void get_univariate_coefficients(vector& coeff) const { m.get_univariate_coefficients(root, coeff); } + vector get_univariate_coefficients() const { vector coeff; m.get_univariate_coefficients(root, coeff); return coeff; } + bool is_never_zero() const { return m.is_never_zero(root); } bool var_is_leaf(unsigned v) const { return m.var_is_leaf(root, v); } pdd operator-() const { return m.minus(*this); } @@ -359,20 +435,28 @@ namespace dd { pdd operator*(rational const& other) const { return m.mul(other, *this); } pdd operator+(rational const& other) const { return m.add(other, *this); } pdd operator~() const { return m.mk_not(*this); } + pdd shl(unsigned n) const; pdd rev_sub(rational const& r) const { return m.sub(m.mk_val(r), *this); } + pdd div(rational const& other) const { return m.div(*this, other); } + bool try_div(rational const& other, pdd& out_result) const { return m.try_div(*this, other, out_result); } + pdd pow(unsigned j) const { return m.pow(*this, j); } pdd reduce(pdd const& other) const { return m.reduce(*this, other); } bool different_leading_term(pdd const& other) const { return m.different_leading_term(*this, other); } - void factor(unsigned v, unsigned degree, pdd& lc, pdd& rest) { m.factor(*this, v, degree, lc, rest); } + void factor(unsigned v, unsigned degree, pdd& lc, pdd& rest) const { m.factor(*this, v, degree, lc, rest); } + bool factor(unsigned v, unsigned degree, pdd& lc) const { return m.factor(*this, v, degree, lc); } + bool resolve(unsigned v, pdd const& other, pdd& result) { return m.resolve(v, *this, other, result); } + pdd reduce(unsigned v, pdd const& other) const { return m.reduce(v, *this, other); } /** * \brief factor out variables */ std::pair var_factors() const; - pdd subst_val(vector> const& s) const { return m.subst_val(*this, s); } + pdd subst_val0(vector> const& s) const { return m.subst_val0(*this, s); } + pdd subst_val(pdd const& s) const { return m.subst_val(*this, s); } pdd subst_val(unsigned v, rational const& val) const { return m.subst_val(*this, v, val); } + pdd subst_add(unsigned var, rational const& val) { return m.subst_add(*this, var, val); } - /** * \brief substitute variable v by r. */ @@ -381,6 +465,7 @@ namespace dd { std::ostream& display(std::ostream& out) const { return m.display(out, *this); } bool operator==(pdd const& other) const { return root == other.root; } bool operator!=(pdd const& other) const { return root != other.root; } + unsigned hash() const { return root; } unsigned power_of_2() const { return m.power_of_2(); } @@ -388,11 +473,15 @@ namespace dd { double tree_size() const { return m.tree_size(*this); } unsigned degree() const { return m.degree(*this); } unsigned degree(unsigned v) const { return m.degree(root, v); } + unsigned max_pow2_divisor() const { return m.max_pow2_divisor(root); } unsigned_vector const& free_vars() const { return m.free_vars(*this); } + void swap(pdd& other) { std::swap(root, other.root); } pdd_iterator begin() const; pdd_iterator end() const; + + pdd_manager& manager() const { return m; } }; inline pdd operator*(rational const& r, pdd const& b) { return b * r; } @@ -403,23 +492,25 @@ namespace dd { inline pdd operator+(int x, pdd const& b) { return b + rational(x); } inline pdd operator+(pdd const& b, int x) { return b + rational(x); } - inline pdd operator^(unsigned x, pdd const& b) { return b + x; } - inline pdd operator^(bool x, pdd const& b) { return b + x; } + inline pdd operator^(unsigned x, pdd const& b) { return b ^ x; } + inline pdd operator^(bool x, pdd const& b) { return b ^ x; } inline pdd operator-(rational const& r, pdd const& b) { return b.rev_sub(r); } inline pdd operator-(int x, pdd const& b) { return rational(x) - b; } inline pdd operator-(pdd const& b, int x) { return b + (-rational(x)); } inline pdd operator-(pdd const& b, rational const& r) { return b + (-r); } - + inline pdd& operator&=(pdd & p, pdd const& q) { p = p & q; return p; } inline pdd& operator^=(pdd & p, pdd const& q) { p = p ^ q; return p; } - inline pdd& operator*=(pdd & p, pdd const& q) { p = p * q; return p; } inline pdd& operator|=(pdd & p, pdd const& q) { p = p | q; return p; } + inline pdd& operator*=(pdd & p, pdd const& q) { p = p * q; return p; } inline pdd& operator-=(pdd & p, pdd const& q) { p = p - q; return p; } inline pdd& operator+=(pdd & p, pdd const& q) { p = p + q; return p; } - inline pdd& operator+=(pdd & p, rational const& v) { p = p + v; return p; } - inline pdd& operator-=(pdd & p, rational const& v) { p = p - v; return p; } - inline pdd& operator*=(pdd & p, rational const& v) { p = p * v; return p; } + inline pdd& operator*=(pdd & p, rational const& q) { p = p * q; return p; } + inline pdd& operator-=(pdd & p, rational const& q) { p = p - q; return p; } + inline pdd& operator+=(pdd & p, rational const& q) { p = p + q; return p; } + + inline void swap(pdd& p, pdd& q) { p.swap(q); } std::ostream& operator<<(std::ostream& out, pdd const& b); diff --git a/src/test/pdd.cpp b/src/test/pdd.cpp index 9139e499f..a0946d81d 100644 --- a/src/test/pdd.cpp +++ b/src/test/pdd.cpp @@ -4,7 +4,7 @@ namespace dd { class test { -public : +public: static void hello_world() { pdd_manager m(3); @@ -310,17 +310,288 @@ public : sub1.push_back(std::make_pair(va, rational(1))); sub2.push_back(std::make_pair(va, rational(2))); sub3.push_back(std::make_pair(va, rational(3))); - std::cout << "sub 0 " << p.subst_val(sub0) << "\n"; - std::cout << "sub 1 " << p.subst_val(sub1) << "\n"; - std::cout << "sub 2 " << p.subst_val(sub2) << "\n"; - std::cout << "sub 3 " << p.subst_val(sub3) << "\n"; + std::cout << "sub 0 " << p.subst_val0(sub0) << "\n"; + std::cout << "sub 1 " << p.subst_val0(sub1) << "\n"; + std::cout << "sub 2 " << p.subst_val0(sub2) << "\n"; + std::cout << "sub 3 " << p.subst_val0(sub3) << "\n"; - std::cout << "expect 1 " << (2*a + 1).is_non_zero() << "\n"; - std::cout << "expect 1 " << (2*a*b + 2*b + 1).is_non_zero() << "\n"; - std::cout << "expect 0 " << (2*a + 2).is_non_zero() << "\n"; - SASSERT((2*a + 1).is_non_zero()); - SASSERT((2*a + 2*b + 1).is_non_zero()); - SASSERT(!(2*a*b + 3*b + 2).is_non_zero()); + std::cout << "expect 1 " << (2*a + 1).is_never_zero() << "\n"; + std::cout << "expect 1 " << (2*a*b + 2*b + 1).is_never_zero() << "\n"; + std::cout << "expect 0 " << (2*a + 2).is_never_zero() << "\n"; + VERIFY((2*a + 1).is_never_zero()); + VERIFY((2*a + 2*b + 1).is_never_zero()); + VERIFY(!(2*a*b + 3*b + 2).is_never_zero()); + } + + static void degree_of_variables() { + std::cout << "degree of variables\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 3); + unsigned va = 0; + unsigned vb = 1; + unsigned vc = 2; + pdd a = m.mk_var(va); + pdd b = m.mk_var(vb); + pdd c = m.mk_var(vc); + + VERIFY(a.var() == va); + VERIFY(b.var() == vb); + + VERIFY(a.degree(va) == 1); + VERIFY(a.degree(vb) == 0); + VERIFY(a.degree(vc) == 0); + VERIFY(c.degree(vc) == 1); + VERIFY(c.degree(vb) == 0); + VERIFY(c.degree(va) == 0); + + { + pdd p = a * a * a; + VERIFY(p.degree(va) == 3); + VERIFY(p.degree(vb) == 0); + } + + { + pdd p = b * a; + VERIFY(p.degree(va) == 1); + VERIFY(p.degree(vb) == 1); + VERIFY(p.degree(vc) == 0); + } + + { + pdd p = (a*a*b + b*a*b + b + a*c)*a + b*b*c; + VERIFY(p.degree(va) == 3); + VERIFY(p.degree(vb) == 2); + VERIFY(p.degree(vc) == 1); + } + + { + // check that skipping marked nodes works (1) + pdd p = b*a + c*a*a*a; + VERIFY(p.degree(va) == 3); + } + + { + // check that skipping marked nodes works (2) + pdd p = (b+c)*(a*a*a); + VERIFY(p.degree(va) == 3); + } + + { + // check that skipping marked nodes works (3) + pdd p = a*a*a*b*b*b*c + a*a*a*b*b*b; + VERIFY(p.degree(va) == 3); + } + } + + static void factor() { + std::cout << "factor\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 3); + + unsigned const va = 0; + unsigned const vb = 1; + unsigned const vc = 2; + unsigned const vd = 3; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + pdd const c = m.mk_var(vc); + pdd const d = m.mk_var(vd); + + auto test_one = [&m](pdd const& p, unsigned v, unsigned d) { + pdd lc = m.zero(); + pdd rest = m.zero(); + std::cout << "Factoring p = " << p << " by v" << v << "^" << d << "\n"; + p.factor(v, d, lc, rest); + std::cout << " lc = " << lc << "\n"; + std::cout << " rest = " << rest << "\n"; + pdd x = m.mk_var(v); + pdd x_pow_d = m.one(); + for (unsigned i = 0; i < d; ++i) { + x_pow_d *= x; + } + VERIFY( p == lc * x_pow_d + rest ); + VERIFY( d == 0 || rest.degree(v) < d ); + VERIFY( d != 0 || lc.is_zero() ); + }; + + auto test_multiple = [=](pdd const& p) { + for (auto v : {va, vb, vc, vd}) { + for (unsigned d = 0; d <= 5; ++d) { + test_one(p, v, d); + } + } + }; + + test_multiple( b ); + test_multiple( b*b*b ); + test_multiple( b + c ); + test_multiple( a*a*a*a*a + a*a*a*b + a*a*b*b + c ); + test_multiple( c*c*c*c*c + b*b*b*c + 3*b*c*c + a ); + test_multiple( (a + b) * (b + c) * (c + d) * (d + a) ); + } + + static void max_pow2_divisor() { + std::cout << "max_pow2_divisor\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 256); + + unsigned const va = 0; + unsigned const vb = 1; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + + VERIFY(m.zero().max_pow2_divisor() == UINT_MAX); + VERIFY(m.one().max_pow2_divisor() == 0); + pdd p = (1 << 20)*a*b + 1024*b*b*b; + std::cout << p << " divided by 2^" << p.max_pow2_divisor() << "\n"; + VERIFY(p.max_pow2_divisor() == 10); + VERIFY(p.div(rational::power_of_two(10)) == 1024*a*b + b*b*b); + VERIFY((p + p).max_pow2_divisor() == 11); + VERIFY((p * p).max_pow2_divisor() == 20); + VERIFY((p + 2*b).max_pow2_divisor() == 1); + VERIFY((p + b*b*b).max_pow2_divisor() == 0); + } + + static void try_div() { + std::cout << "try_div\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 256); + pdd const a = m.mk_var(0); + pdd const b = m.mk_var(1); + + pdd const p = 5*a + 15*a*b; + VERIFY_EQ(p.div(rational(5)), a + 3*a*b); + pdd res = a; + VERIFY(!p.try_div(rational(3), res)); + } + + static void binary_resolve() { + std::cout << "binary resolve\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 4); + + unsigned const va = 0; + unsigned const vb = 1; + unsigned const vc = 2; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + pdd const c = m.mk_var(vc); + + pdd r = m.zero(); + + pdd p = a*a*b - a*a; + pdd q = a*b*b - b*b; + VERIFY(m.resolve(va, p, q, r)); + VERIFY(r == a*b*b*b - a*b*b); + VERIFY(!m.resolve(va, q, p, r)); + VERIFY(!m.resolve(vb, p, q, r)); + VERIFY(m.resolve(vb, q, p, r)); + VERIFY(r == a*a*a*b - a*a*b); + VERIFY(!m.resolve(vc, p, q, r)); + + p = 2*a*a*b + 13*a*a; + q = 6*a*b*b*b + 14*b*b*b; + VERIFY(m.resolve(va, p, q, r)); + VERIFY(r == (2*b+13)*2*b*b*b*a); + VERIFY(!m.resolve(va, q, p, r)); + VERIFY(!m.resolve(vb, p, q, r)); + VERIFY(m.resolve(vb, q, p, r)); + VERIFY(r == 9*a*a*a*b*b + 5*a*a*b*b); + + p = a*a*b - a*a + 4*a*c + 2; + q = 3*b*b - b*b*b + 8*b*c; + VERIFY(!m.resolve(va, p, q, r)); + VERIFY(!m.resolve(va, q, p, r)); + VERIFY(!m.resolve(vb, p, q, r)); + VERIFY(m.resolve(vb, q, p, r)); + VERIFY(r == 2*a*a*b*b + 8*a*a*b*c + 4*a*b*b*c + 2*b*b); + VERIFY(m.resolve(vc, p, q, r)); + VERIFY(r == 2*a*a*b*b - 2*a*a*b - 3*a*b*b + a*b*b*b + 4*b); + VERIFY(m.resolve(vc, q, p, r)); + VERIFY(r == -(2*a*a*b*b - 2*a*a*b - 3*a*b*b + a*b*b*b + 4*b)); + } + + static void pow() { + std::cout << "pow\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 5); + + unsigned const va = 0; + unsigned const vb = 1; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + + VERIFY(a.pow(0) == m.one()); + VERIFY(a.pow(1) == a); + VERIFY(a.pow(2) == a*a); + VERIFY(a.pow(7) == a*a*a*a*a*a*a); + VERIFY((3*a*b).pow(3) == 27*a*a*a*b*b*b); + } + + static void subst_val() { + std::cout << "subst_val\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 2); + + unsigned const va = 0; + unsigned const vb = 1; + unsigned const vc = 2; + unsigned const vd = 3; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + pdd const c = m.mk_var(vc); + pdd const d = m.mk_var(vd); + + { + pdd const p = 2*a + b + 1; + VERIFY(p.subst_val(va, rational(0)) == b + 1); + } + + { + pdd const p = a + 2*b; + VERIFY(p.subst_val(va, rational(0)) == 2*b); + VERIFY(p.subst_val(va, rational(2)) == 2*b + 2); + VERIFY(p.subst_val(vb, rational(0)) == a); + VERIFY(p.subst_val(vb, rational(1)) == a + 2); + VERIFY(p.subst_val(vb, rational(2)) == a); + VERIFY(p.subst_val(vb, rational(3)) == a + 2); + VERIFY(p.subst_val(va, rational(0)).subst_val(vb, rational(3)) == 2*m.one()); + } + + { + pdd const p = a + b + c + d; + vector> sub; + sub.push_back({vb, rational(2)}); + sub.push_back({vc, rational(3)}); + VERIFY(p.subst_val0(sub) == a + d + 1); + } + + { + pdd const p = (a + b) * (b + c) * (c + d); + vector> sub; + sub.push_back({vb, rational(2)}); + sub.push_back({vc, rational(3)}); + VERIFY(p.subst_val0(sub) == (a + 2) * (d + 3)); + sub.push_back({va, rational(3)}); + sub.push_back({vd, rational(2)}); + VERIFY(p.subst_val0(sub) == m.one()); + } + } + + static void univariate() { + std::cout << "univariate\n"; + pdd_manager m(4, pdd_manager::mod2N_e, 4); + + unsigned const va = 0; + unsigned const vb = 1; + pdd const a = m.mk_var(va); + pdd const b = m.mk_var(vb); + + pdd p = a*a*b - a*a; + VERIFY(!p.is_univariate()); + + pdd q = 3*a*a*a + 1*a + 2; + VERIFY(q.is_univariate()); + vector coeff; + q.get_univariate_coefficients(coeff); + VERIFY_EQ(coeff.size(), 4); + VERIFY_EQ(coeff[0], 2); + VERIFY_EQ(coeff[1], 1); + VERIFY_EQ(coeff[2], 0); + VERIFY_EQ(coeff[3], 3); } static void factors() { @@ -393,5 +664,13 @@ void tst_pdd() { dd::test::order(); dd::test::order_lm(); dd::test::mod4_operations(); + dd::test::degree_of_variables(); + dd::test::factor(); + dd::test::max_pow2_divisor(); + dd::test::try_div(); + dd::test::binary_resolve(); + dd::test::pow(); + dd::test::subst_val(); + dd::test::univariate(); dd::test::factors(); }