From 7ad95aa5d25a75afbf5bd179c48706688a21c6ee Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Fri, 10 Jan 2020 12:02:34 -0800 Subject: [PATCH] Nikolaj fixes pdd_manager::reduce() to work with the changed order Signed-off-by: Lev Nachmanson --- src/math/dd/dd_pdd.cpp | 105 ++++++---- src/math/dd/dd_pdd.h | 6 +- src/math/grobner/pdd_solver.cpp | 4 +- src/math/grobner/pdd_solver.h | 1 + src/test/pdd.cpp | 337 +++++++++++++++++--------------- 5 files changed, 246 insertions(+), 207 deletions(-) diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index 4d492d8fc..aeb90909a 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -184,7 +184,8 @@ namespace dd { case pdd_reduce_op: if (is_zero(q)) return p; if (is_val(p)) return p; - if (!is_val(q) && level(p) < level(q)) return p; + if (degree(p) < degree(q)) return p; + if (level(first_leading(q)) > level(p)) return p; break; case pdd_subst_val_op: while (!is_val(q) && !is_val(p)) { @@ -310,18 +311,25 @@ namespace dd { } break; case pdd_reduce_op: - if (level_p > level_q) { + if (level(first_leading(q)) < level_p) { push(apply_rec(lo(p), q, op)); push(apply_rec(hi(p), q, op)); - if (read(2) == lo(p) && read(1) == hi(p)) { + PDD plo = read(2), phi = read(1); + if (plo == lo(p) && phi == hi(p)) { r = p; } - else { - r = make_node(level_p, read(2), read(1)); + else if (level(plo) < level_p && level(phi) <= level(p)) { + r = make_node(level_p, plo, phi); } + else { + push(apply_rec(phi, m_var2pdd[var(p)], pdd_mul_op)); + push(apply_rec(read(1), plo, pdd_add_op)); + r = read(1); + npop = 4; + } } else { - SASSERT(level_p == level_q); + SASSERT(level(first_leading(q)) == level_p); r = reduce_on_match(p, q); npop = 0; } @@ -384,11 +392,16 @@ namespace dd { return r; } - // q = lt(a)/lt(b), return a - b*q + // + // produce polynomial where a is reduced by b. + // all monomials in a that are divisible by lm(b) + // are replaced by (b - lt(b))/lm(b) + // pdd_manager::PDD pdd_manager::reduce_on_match(PDD a, PDD b) { - SASSERT(level(a) == level(b) && !is_val(a) && !is_val(b)); + SASSERT(level(first_leading(b)) == level(a)); + SASSERT(!is_val(a) && !is_val(b)); push(a); - while (lm_divides(b, a)) { + while (lm_occurs(b, a)) { push(lt_quotient(b, a)); push(apply_rec(read(1), b, pdd_mul_op)); push(apply_rec(a, read(1), pdd_add_op)); @@ -400,47 +413,65 @@ namespace dd { return a; } - // true if leading monomial of p divides leading monomial of q - bool pdd_manager::lm_divides(PDD p, PDD q) const { + // true if leading monomial of p divides a monomial of q + bool pdd_manager::lm_occurs(PDD p, PDD q) const { p = first_leading(p); - q = first_leading(q); while (true) { if (is_val(p)) return true; if (is_val(q)) return false; - if (level(p) > level(q)) return false; - if (level(p) == level(q)) { + if (level(p) > level(q)) { + return false; + } + else if (level(p) == level(q)) { p = next_leading(p); - q = next_leading(q); + q = hi(q); + } + else if (lm_occurs(p, hi(q))) { + return true; } else { - q = next_leading(q); + q = lo(q); } } } - // return minus quotient -r, such that lt(q) = lt(p)*r - // assume lm_divides(p, q) + // return minus quotient -r, such that lt(p)*r + // is a monomial in q. + // assume lm_occurs(p, q) pdd_manager::PDD pdd_manager::lt_quotient(PDD p, PDD q) { - SASSERT(lm_divides(p, q)); + SASSERT(lm_occurs(p, q)); p = first_leading(p); - q = first_leading(q); SASSERT(is_val(p) || !is_val(q)); - if (is_val(p)) { - if (is_val(q)) { - SASSERT(!val(p).is_zero()); - return imk_val(-val(q) / val(p)); + + while (!is_val(p)) { + SASSERT(level(p) <= level(q)); + SASSERT(lm_occurs(p, q)); + if (level(p) == level(q)) { + p = next_leading(p); + q = lm_occurs(p, hi(q)) ? hi(q) : lo(q); + } + else if (lm_occurs(p, hi(q))) { + return lt_quotient_hi(p, q); + } + else { + q = lo(q); } } - else if (level(p) == level(q)) { - return lt_quotient(next_leading(p), next_leading(q)); + SASSERT(!is_zero(p)); + if (is_val(q)) { + return imk_val(-val(q) / val(p)); } + else { + return lt_quotient_hi(p, q); + } + } - SASSERT(!is_val(q)); - push(lt_quotient(p, next_leading(q))); + pdd_manager::PDD pdd_manager::lt_quotient_hi(PDD p, PDD q) { + SASSERT(lm_occurs(p, hi(q))); + push(lt_quotient(p, hi(q))); PDD r = apply_rec(m_var2pdd[var(q)], read(1), pdd_mul_op); pop(1); return r; - } // @@ -606,19 +637,6 @@ namespace dd { return p; } - /* - Determine whether p is a linear polynomials. - A linear polynomial is of the form x*v1 + y*v2 + .. + vn, - where v1, v2, .., vn are values. - */ - bool pdd_manager::is_linear(PDD p) { - while (true) { - if (is_val(p)) return true; - if (!is_val(hi(p))) return false; - p = lo(p); - } - } - bool pdd_manager::is_linear(pdd const& p) { return is_linear(p.root); } @@ -853,6 +871,9 @@ namespace dd { } unsigned pdd_manager::degree(PDD p) const { + if (p == zero_pdd || p == one_pdd) { + return 0; + } if (is_dmarked(p)) { return m_degree[p]; } diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index fe9b0ecf9..3934f549e 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -46,7 +46,6 @@ namespace dd { public: enum semantics { free_e, mod2_e, zero_one_vars_e }; private: - friend test; friend pdd; friend pdd_iterator; @@ -178,8 +177,9 @@ namespace dd { PDD minus_rec(PDD p); PDD reduce_on_match(PDD a, PDD b); - bool lm_divides(PDD p, PDD q) const; + bool lm_occurs(PDD p, PDD q) const; PDD lt_quotient(PDD p, PDD q); + PDD lt_quotient_hi(PDD p, PDD q); PDD imk_val(rational const& r); void init_value(const_info& info, rational const& r); @@ -279,7 +279,7 @@ namespace dd { pdd reduce(pdd const& a, pdd const& b); pdd subst_val(pdd const& a, vector> const& s); - bool is_linear(PDD p); + bool is_linear(PDD p) { return degree(p) == 1; } bool is_linear(pdd const& p); bool is_binary(PDD p); diff --git a/src/math/grobner/pdd_solver.cpp b/src/math/grobner/pdd_solver.cpp index 2693c003b..2c0f4de53 100644 --- a/src/math/grobner/pdd_solver.cpp +++ b/src/math/grobner/pdd_solver.cpp @@ -110,7 +110,7 @@ namespace dd { void solver::scoped_process::done() { pdd p = e->poly(); SASSERT(!p.is_val()); - if (p.hi().is_val()) { + if (p.degree() == 1) { g.push_equation(solved, e); } else { @@ -462,7 +462,7 @@ namespace dd { VERIFY(e->idx() == i); ++i; pdd p = e->poly(); - if (!p.is_val() && p.hi().is_val()) { + if (p.degree() == 1) { unsigned v = p.var(); SASSERT(!head_vars.contains(v)); head_vars.insert(v); diff --git a/src/math/grobner/pdd_solver.h b/src/math/grobner/pdd_solver.h index 35f918779..b8fc04b0a 100644 --- a/src/math/grobner/pdd_solver.h +++ b/src/math/grobner/pdd_solver.h @@ -194,6 +194,7 @@ private: scoped_process(solver& g, equation* e): g(g), e(e) {} ~scoped_process(); }; + void update_stats_max_degree_and_size(const equation& e); }; diff --git a/src/test/pdd.cpp b/src/test/pdd.cpp index ae07fb476..216b3d48c 100644 --- a/src/test/pdd.cpp +++ b/src/test/pdd.cpp @@ -1,147 +1,163 @@ #include "math/dd/dd_pdd.h" namespace dd { -static void test1() { - pdd_manager m(3); - pdd v0 = m.mk_var(0); - pdd v1 = m.mk_var(1); - pdd v2 = m.mk_var(2); - std::cout << v0 << "\n"; - std::cout << v1 << "\n"; - std::cout << v2 << "\n"; - pdd c1 = v0 * v1 * v2; - pdd c2 = v2 * v0 * v1; - std::cout << c1 << "\n"; - SASSERT(c1 == c2); - c1 = v0 + v1 + v2; - c2 = v2 + v1 + v0; - std::cout << c1 << "\n"; - SASSERT(c1 == c2); - - c1 = (v0+v1) * v2; - c2 = (v0*v2) + (v1*v2); - std::cout << c1 << "\n"; - SASSERT(c1 == c2); - c1 = (c1 + 3) + 1; - c2 = (c2 + 1) + 3; - std::cout << c1 << "\n"; - SASSERT(c1 == c2); - c1 = v0 - v1; - c2 = v1 - v0; - std::cout << c1 << " " << c2 << "\n"; - - c1 = v1*v2; - c2 = (v0*v2) + (v2*v2); - pdd c3 = m.zero(); - VERIFY(m.try_spoly(c1, c2, c3)); - std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n"; - - c1 = v1*v2; - c2 = (v0*v2) + (v1*v1); - VERIFY(m.try_spoly(c1, c2, c3)); - std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n"; - - c1 = (v0*v1) - (v0*v0); - c2 = (v0*v1*(v2 + v0)) + v2; - c3 = c2.reduce(c1); - std::cout << c1 << " " << c2 << " reduce: " << c3 << "\n"; -} - -static void test2() { - std::cout << "\ntest2\n"; - // a(b^2)cd + abc + bcd + bc + cd + 3 reduce by bc - pdd_manager m(4); - pdd a = m.mk_var(0); - pdd b = m.mk_var(1); - pdd c = m.mk_var(2); - pdd d = m.mk_var(3); - pdd e = (a * b * b * c * d) + (2*a*b*c) + (b*c*d) + (b*c) + (c*d) + 3; - std::cout << e << "\n"; - pdd f = b * c; - pdd r_ef = m.reduce(e, f); - m.display(std::cout); - std::cout << "result of reduce " << e << " by " << f << " is " << r_ef << "\n"; - pdd r_fe = m.reduce(f, e); - std::cout << "result of reduce " << f << " by " << e << " is " << r_fe << "\n" ; - VERIFY(r_fe == f); -} - -static void test3() { - std::cout << "\ntest3\n"; - pdd_manager m(4); - pdd a = m.mk_var(0); - pdd b = m.mk_var(1); - pdd c = m.mk_var(2); - pdd d = m.mk_var(3); - - pdd e = a + c; - for (unsigned i = 0; i < 5; i++) { - e = e * e; - } - e = e * b; - std::cout << e << "\n"; -} - -static void test_reset() { - std::cout << "\ntest reset\n"; - pdd_manager m(4); - pdd a = m.mk_var(0); - pdd b = m.mk_var(1); - pdd c = m.mk_var(2); - pdd d = m.mk_var(3); - std::cout << (a + b)*(c + d) << "\n"; - - unsigned_vector l2v; - for (unsigned i = 0; i < 4; ++i) - l2v.push_back(3 - i); - m.reset(l2v); - a = m.mk_var(0); - b = m.mk_var(1); - c = m.mk_var(2); - d = m.mk_var(3); - std::cout << (a + b)*(c + d) << "\n"; -} - -static void test5() { - std::cout << "\ntest5\n"; - pdd_manager m(2); - pdd a = m.mk_var(0); - pdd b = m.mk_var(1); - - pdd e = (a - b) * ( a + b); - pdd f = a * a - b * b; - SASSERT(e == f); - - e = (a - b)* (a - b); - f = a * a - 2 * a * b + b * b; - SASSERT(e == f); - e = a - 3; - e = e * e; - f = a * a - 6 * a + 9; - SASSERT(e == f); - e = 2 * a - 3; - e = e * e; - f = 4 * a * a - 12 * a + 9; - SASSERT(e == f); -} - - -void test_iterator() { - std::cout << "test iterator\n"; - pdd_manager m(4); - pdd a = m.mk_var(0); - pdd b = m.mk_var(1); - pdd c = m.mk_var(2); - pdd d = m.mk_var(3); - pdd p = (a + b)*(c + 3*d) + 2; - std::cout << p << "\n"; - for (auto const& m : p) { - std::cout << m << "\n"; - } -} class test { public : + + static void hello_world() { + pdd_manager m(3); + pdd v0 = m.mk_var(0); + pdd v1 = m.mk_var(1); + pdd v2 = m.mk_var(2); + std::cout << v0 << "\n"; + std::cout << v1 << "\n"; + std::cout << v2 << "\n"; + pdd c1 = v0 * v1 * v2; + pdd c2 = v2 * v0 * v1; + std::cout << c1 << "\n"; + VERIFY(c1 == c2); + + c1 = v0 + v1 + v2; + c2 = v2 + v1 + v0; + std::cout << c1 << "\n"; + VERIFY(c1 == c2); + + c1 = (v0+v1) * v2; + c2 = (v0*v2) + (v1*v2); + std::cout << c1 << "\n"; + VERIFY(c1 == c2); + c1 = (c1 + 3) + 1; + c2 = (c2 + 1) + 3; + std::cout << c1 << "\n"; + VERIFY(c1 == c2); + c1 = v0 - v1; + c2 = v1 - v0; + std::cout << c1 << " " << c2 << "\n"; + + c1 = v1*v2; + c2 = (v0*v2) + (v2*v2); + pdd c3 = m.zero(); + VERIFY(m.try_spoly(c1, c2, c3)); + std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n"; + + c1 = v1*v2; + c2 = (v0*v2) + (v1*v1); + VERIFY(m.try_spoly(c1, c2, c3)); + std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n"; + + c1 = (v0*v1) - (v0*v0); + c2 = (v0*v1*(v2 + v0)) + v2; + c3 = c2.reduce(c1); + std::cout << c1 << " " << c2 << " reduce: " << c3 << "\n"; + VERIFY(c3 == c3.reduce(c1)); + } + + static void reduce() { + std::cout << "\nreduce\n"; + // a(b^2)cd + abc + bcd + bc + cd + 3 reduce by bc + pdd_manager m(4); + pdd a = m.mk_var(0); + pdd b = m.mk_var(1); + pdd c = m.mk_var(2); + pdd d = m.mk_var(3); + pdd e = (a * b * b * c * d) + (2*a*b*c) + (b*c*d) + (b*c) + (c*d) + 3; + std::cout << e << "\n"; + pdd f = b * c; + pdd r_ef = m.reduce(e, f); + m.display(std::cout); + std::cout << "result of reduce " << e << " by " << f << " is " << r_ef << "\n"; + VERIFY(r_ef == (c*d) + 3); + pdd r_fe = m.reduce(f, e); + std::cout << "result of reduce " << f << " by " << e << " is " << r_fe << "\n" ; + VERIFY(r_fe == f); + + // test that b*c*d is treated as the leading monomial + f = b*c*d - d*d; + r_ef = m.reduce(e, f); + std::cout << "result of reduce " << e << " by " << f << " is " << r_ef << "\n"; + VERIFY(r_ef == (a*b*d*d) + (2*a*b*c) + (d*d) + (b*c) + (c*d) + 3); + + pdd k = d*d + 3*b; + VERIFY(m.reduce(f, k) == b*c*d + 3*b); + + } + + + + static void large_product() { + std::cout << "\nlarge_product\n"; + pdd_manager m(4); + pdd a = m.mk_var(0); + pdd b = m.mk_var(1); + pdd c = m.mk_var(2); + pdd d = m.mk_var(3); + + pdd e = a + c; + for (unsigned i = 0; i < 5; i++) { + e = e * e; + } + e = e * b; + std::cout << e << "\n"; + } + + static void reset() { + std::cout << "\ntest reset\n"; + pdd_manager m(4); + pdd a = m.mk_var(0); + pdd b = m.mk_var(1); + pdd c = m.mk_var(2); + pdd d = m.mk_var(3); + std::cout << (a + b)*(c + d) << "\n"; + + unsigned_vector l2v; + for (unsigned i = 0; i < 4; ++i) + l2v.push_back(3 - i); + m.reset(l2v); + a = m.mk_var(0); + b = m.mk_var(1); + c = m.mk_var(2); + d = m.mk_var(3); + std::cout << (a + b)*(c + d) << "\n"; + } + + static void canonize() { + std::cout << "\ncanonize\n"; + pdd_manager m(2); + pdd a = m.mk_var(0); + pdd b = m.mk_var(1); + + pdd e = (a - b) * ( a + b); + pdd f = a * a - b * b; + VERIFY(e == f); + + e = (a - b)* (a - b); + f = a * a - 2 * a * b + b * b; + VERIFY(e == f); + e = a - 3; + e = e * e; + f = a * a - 6 * a + 9; + VERIFY(e == f); + e = 2 * a - 3; + e = e * e; + f = 4 * a * a - 12 * a + 9; + VERIFY(e == f); + } + + + static void iterator() { + std::cout << "test iterator\n"; + pdd_manager m(4); + pdd a = m.mk_var(0); + pdd b = m.mk_var(1); + pdd c = m.mk_var(2); + pdd d = m.mk_var(3); + pdd p = (a + b)*(c + 3*d) + 2; + std::cout << p << "\n"; + for (auto const& m : p) { + std::cout << m << "\n"; + } + } static void order() { std::cout << "order\n"; pdd_manager m(4); @@ -158,8 +174,8 @@ public : unsigned i = 0; for (auto const& m : p) { std::cout << m << "\n"; - SASSERT(i != 0 ||( m.vars.size() == 1 && m.vars[0] == vb)); - SASSERT(i != 1 ||( m.vars.size() == 1 && m.vars[0] == va)); + VERIFY(i != 0 ||( m.vars.size() == 1 && m.vars[0] == vb)); + VERIFY(i != 1 ||( m.vars.size() == 1 && m.vars[0] == va)); i++; } pdd ccbbaa = c*c*b*b*a*a; @@ -170,8 +186,8 @@ public : std::cout << "p = " << p << "\n"; for (auto const& m : p) { std::cout << m << "\n"; - SASSERT(i != 0 ||( m.vars.size() == 6 && m.vars[4] == vb)); // the first one has to be ccbbba - SASSERT(i != 1 ||( m.vars.size() == 6 && m.vars[4] == va)); // the second one has to be ccbbaa + VERIFY(i != 0 ||( m.vars.size() == 6 && m.vars[4] == vb)); // the first one has to be ccbbba + VERIFY(i != 1 ||( m.vars.size() == 6 && m.vars[4] == va)); // the second one has to be ccbbaa i++; } pdd dcbba = d*c*b*b*a; @@ -184,7 +200,7 @@ public : k = m.first_leading(k); do { if (m.is_val(k)) { - SASSERT(m.val(k).is_one()); + VERIFY(m.val(k).is_one()); break; } v.push_back(m.var(k)); @@ -198,8 +214,8 @@ public : for( ; it != v.end(); it ++) { std::cout << "*v" << *it; } - SASSERT(v.size() == 6); - SASSERT(v[0] == vc); + VERIFY(v.size() == 6); + VERIFY(v[0] == vc); std::cout << "\n"; v.reset(); p = d*c*c*d + b*c*c*b + d*d*d; @@ -208,7 +224,7 @@ public : k = m.first_leading(k); do { if (m.is_val(k)) { - SASSERT(m.val(k).is_one()); + VERIFY(m.val(k).is_one()); break; } v.push_back(m.var(k)); @@ -222,8 +238,8 @@ public : for( ; it != v.end(); it ++) { std::cout << "*v" << *it; } - SASSERT(v.size() == 4); - SASSERT(v[0] == vd); + VERIFY(v.size() == 4); + VERIFY(v[0] == vd); std::cout << "\n"; } @@ -243,8 +259,8 @@ public : unsigned i = 0; for (auto const& m : p) { std::cout << m << "\n"; - SASSERT(i != 0 ||( m.vars.size() == 1 && m.vars[0] == vb)); - SASSERT(i != 1 ||( m.vars.size() == 1 && m.vars[0] == va)); + VERIFY(i != 0 ||( m.vars.size() == 1 && m.vars[0] == vb)); + VERIFY(i != 1 ||( m.vars.size() == 1 && m.vars[0] == va)); i++; } pdd ccbbaa = c*c*b*b*a*a; @@ -255,29 +271,30 @@ public : std::cout << "p = " << p << "\n"; for (auto const& m : p) { std::cout << m << "\n"; - SASSERT(i != 0 ||( m.vars.size() == 6 && m.vars[4] == vb)); // the first one has to be ccbbba - SASSERT(i != 1 ||( m.vars.size() == 6 && m.vars[4] == va)); // the second one has to be ccbbaa + VERIFY(i != 0 ||( m.vars.size() == 6 && m.vars[4] == vb)); // the first one has to be ccbbba + VERIFY(i != 1 ||( m.vars.size() == 6 && m.vars[4] == va)); // the second one has to be ccbbaa i++; } pdd dcbba = d*c*b*b*a; pdd dd = d * d; pdd p0 = p + dd; std::cout << p0 << " > " << p << "\n"; - SASSERT(m.lm_lt(p, p0)); - SASSERT(m.lm_lt(p0 + a * b, p0 + b * b)); + VERIFY(m.lm_lt(p, p0)); + VERIFY(m.lm_lt(p0 + a * b, p0 + b * b)); } + }; } void tst_pdd() { - dd::test1(); - dd::test2(); - dd::test3(); - dd::test5(); - dd::test_reset(); - dd::test_iterator(); + dd::test::hello_world(); + dd::test::reduce(); + dd::test::large_product(); + dd::test::canonize(); + dd::test::reset(); + dd::test::iterator(); dd::test::order(); dd::test::order_lm(); }