diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index fdfa7ec41..fe3eb0b02 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -696,8 +696,8 @@ namespace dd { 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(); + lc = zero(); + rest = p; return; } if (level(p.root) < level_v) { @@ -707,12 +707,13 @@ namespace dd { } // 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); + factor_entry* e = &et->get_data(); + if (e->is_valid()) { + lc = pdd(e->m_lc, this); + rest = pdd(e->m_rest, this); return; } + unsigned const gc_generation = m_gc_generation; if (level(p.root) > level_v) { pdd lc1 = zero(), rest1 = zero(); pdd vv = mk_var(p.var()); @@ -743,10 +744,118 @@ namespace dd { rest = p; } } - e.m_lc = lc.root; - e.m_rest = rest.root; + if (gc_generation != m_gc_generation) { + // Cache was reset while factoring (due to GC), + // which means the old entry has been removed and we need to insert it again. + auto* 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; } + 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); + if (l < m) { + // no reduction + return false; + } + if (m == 0) { + // no reduction (result would still contain v^l) + 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); + auto div_pow2j = [&pow2j](rational const& r) -> rational { + rational result = r / pow2j; + SASSERT(result.is_int()); + return result; + }; + pdd aa = map_coefficients(a, div_pow2j); + pdd cc = map_coefficients(c, div_pow2j); + pdd vv = one(); + for (unsigned deg = l - m; deg-- > 0; ) { + vv *= mk_var(v); + } + r = b * cc - aa * d * vv; + return true; + } + + /** + * 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); @@ -1142,6 +1251,7 @@ namespace dd { } void pdd_manager::gc() { + m_gc_generation++; init_dmark(); m_free_nodes.reset(); SASSERT(well_formed()); diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index 1ae6ba02f..97355970a 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -200,6 +200,7 @@ namespace dd { rational m_freeze_value; rational m_mod2N; unsigned m_power_of_2 { 0 }; + unsigned m_gc_generation { 0 }; ///< will be incremented on each GC void reset_op_cache(); void init_nodes(unsigned_vector const& l2v); @@ -261,6 +262,9 @@ 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); @@ -320,6 +324,7 @@ namespace dd { pdd reduce(pdd const& a, pdd const& b); pdd subst_val(pdd const& a, vector> const& s); pdd subst_val(pdd const& a, unsigned v, rational const& val); + bool resolve(unsigned v, pdd const& p, pdd const& q, pdd& r); bool is_linear(PDD p) { return degree(p) == 1; } bool is_linear(pdd const& p); @@ -411,6 +416,7 @@ 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); } diff --git a/src/test/pdd.cpp b/src/test/pdd.cpp index 8079d5419..4fab2f276 100644 --- a/src/test/pdd.cpp +++ b/src/test/pdd.cpp @@ -408,7 +408,7 @@ public : } SASSERT( p == lc * x_pow_d + rest ); SASSERT( d == 0 || rest.degree(v) < d ); - SASSERT( d != 0 || rest.is_zero() ); + SASSERT( d != 0 || lc.is_zero() ); }; auto test_multiple = [=](pdd const& p) { @@ -427,6 +427,71 @@ public : 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); + + SASSERT(m.zero().max_pow2_divisor() == UINT_MAX); + SASSERT(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"; + SASSERT(p.max_pow2_divisor() == 10); + SASSERT((p + p).max_pow2_divisor() == 11); + SASSERT((p * p).max_pow2_divisor() == 20); + SASSERT((p + 2*b).max_pow2_divisor() == 1); + SASSERT((p + b*b*b).max_pow2_divisor() == 0); + } + + 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; + SASSERT(m.resolve(va, p, q, r)); + SASSERT(r == a*b*b*b - a*b*b); + SASSERT(!m.resolve(va, q, p, r)); + SASSERT(!m.resolve(vb, p, q, r)); + SASSERT(m.resolve(vb, q, p, r)); + SASSERT(r == a*a*a*b - a*a*b); + SASSERT(!m.resolve(vc, p, q, r)); + + p = 2*a*a*b + 13*a*a; + q = 6*a*b*b*b + 14*b*b*b; + SASSERT(m.resolve(va, p, q, r)); + SASSERT(r == (2*b+13)*2*b*b*b*a); + SASSERT(!m.resolve(va, q, p, r)); + SASSERT(!m.resolve(vb, p, q, r)); + SASSERT(m.resolve(vb, q, p, r)); + SASSERT(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; + SASSERT(!m.resolve(va, p, q, r)); + SASSERT(!m.resolve(va, q, p, r)); + SASSERT(!m.resolve(vb, p, q, r)); + SASSERT(m.resolve(vb, q, p, r)); + SASSERT(r == 2*a*a*b*b + 8*a*a*b*c + 4*a*b*b*c + 2*b*b); + SASSERT(m.resolve(vc, p, q, r)); + SASSERT(r == 2*a*a*b*b - 2*a*a*b - 3*a*b*b + a*b*b*b + 4*b); + SASSERT(m.resolve(vc, q, p, r)); + SASSERT(r == -(2*a*a*b*b - 2*a*a*b - 3*a*b*b + a*b*b*b + 4*b)); + } + }; } @@ -443,4 +508,6 @@ void tst_pdd() { dd::test::mod4_operations(); dd::test::degree_of_variables(); dd::test::factor(); + dd::test::max_pow2_divisor(); + dd::test::binary_resolve(); }