adding pdd-grobner

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/dd/dd_pdd.cpp

@@ -338,6 +338,56 @@ namespace dd {
         }
     }
 
+    /*
+     * Compare leading monomials.
+     * The pdd format makes lexicographic comparison easy: compare based on
+     * the top variable and descend depending on whether hi(x) == hi(y)
+     */
+    bool pdd_manager::lt(pdd const& a, pdd const& b) {
+        PDD x = a.root;
+        PDD y = b.root;
+        if (x == y) return false;
+        while (true) {
+            SASSERT(x != y);
+            if (is_val(x)) 
+                return !is_val(y) || val(x) < val(y);
+            if (is_val(y)) 
+                return false;
+            if (level(x) == level(y)) {
+                if (hi(x) == hi(y)) {
+                    x = lo(x);
+                    y = lo(y);
+                }
+                else {
+                    x = hi(x);
+                    y = hi(y);
+                }
+            }
+            else {
+                return level(x) > level(y);
+            }
+        }
+    }
+
+    /**
+       Compare leading terms of pdds
+     */
+    bool pdd_manager::different_leading_term(pdd const& a, pdd const& b) {
+        PDD x = a.root;
+        PDD y = b.root;
+        while (true) {
+            if (x == y) return false;
+            if (is_val(x) || is_val(y)) return true;
+            if (level(x) == level(y)) {
+                x = hi(x);
+                y = hi(y);
+            }
+            else {
+                return false;
+            }
+        }
+    }
+
     void pdd_manager::push(PDD b) {
         m_pdd_stack.push_back(b);
     }

src/math/dd/dd_pdd.h

@@ -141,7 +141,6 @@ namespace dd {
         small_object_allocator     m_alloc;
         mutable svector<unsigned>  m_mark;
         mutable unsigned           m_mark_level;
-        mutable svector<double>    m_count;
         mutable svector<PDD>       m_todo;
         bool                       m_disable_gc;
         bool                       m_is_new_node;
@@ -230,6 +229,8 @@ namespace dd {
         pdd mul(rational const& c, pdd const& b);
         pdd reduce(pdd const& a, pdd const& b);
         bool try_spoly(pdd const& a, pdd const& b, pdd& r);
+        bool lt(pdd const& a, pdd const& b);
+        bool different_leading_term(pdd const& a, pdd const& b);
 
         std::ostream& display(std::ostream& out);
         std::ostream& display(std::ostream& out, pdd const& b);
@@ -243,7 +244,8 @@ namespace dd {
         pdd_manager* m;
         pdd(unsigned root, pdd_manager* m): root(root), m(m) { m->inc_ref(root); }
     public:
-        pdd(pdd & other): root(other.root), m(other.m) { m->inc_ref(root); }
+        pdd(pdd_manager& pm): root(0), m(&pm) { SASSERT(is_zero()); }
+        pdd(pdd const& other): root(other.root), m(other.m) { m->inc_ref(root); }
         pdd(pdd && other): root(0), m(other.m) { std::swap(root, other.root); }
         pdd& operator=(pdd const& other);
         ~pdd() { m->dec_ref(root); }
@@ -252,6 +254,7 @@ namespace dd {
         unsigned var() const { return m->var(root); }
         rational const& val() const { SASSERT(is_val()); return m->val(root); }
         bool is_val() const { return m->is_val(root); }
+        bool is_zero() const { return m->is_zero(root); }
 
         pdd operator+(pdd const& other) const { return m->add(*this, other); }
         pdd operator-(pdd const& other) const { return m->sub(*this, other); }
@@ -259,10 +262,13 @@ 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 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); }
 
         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; }
+        bool operator<(pdd const& b) const { return m->lt(*this, b); }
+
         unsigned size() const { return m->pdd_size(*this); }
     };
 
@@ -276,6 +282,7 @@ namespace dd {
 
     inline pdd operator-(pdd const& b, int x) { return b + (-rational(x)); }
 
+
     std::ostream& operator<<(std::ostream& out, pdd const& b);
 
 }

src/math/grobner/CMakeLists.txt

@@ -1,6 +1,7 @@
 z3_add_component(grobner
   SOURCES
     grobner.cpp
+    pdd_grobner.cpp
   COMPONENT_DEPENDENCIES
     ast
 )

src/math/grobner/pdd_grobner.cpp

@@ -0,0 +1,253 @@
+/*++
+  Copyright (c) 2019 Microsoft Corporation
+
+  Abstract:
+
+    Grobner core based on pdd representation of polynomials
+
+  Author:
+     Nikolaj Bjorner (nbjorner)
+     Lev Nachmanson (levnach)
+
+  --*/
+
+#include "math/grobner/pdd_grobner.h"
+
+namespace dd {
+
+    /***
+        The main algorithm maintains two sets (S, A), 
+        where S is m_processed, and A is m_to_simplify.
+        Initially S is empty and A contains the initial equations.
+        
+        Each step proceeds as follows:
+        - pick a in A, and remove a from A
+        - simplify a using S
+        - simplify S using a
+        - for s in S:
+             b = superpose(a, s)
+             add b to A
+        - add a to S
+        - simplify A using a
+    */
+
+    grobner::grobner(reslimit& lim, pdd_manager& m) : 
+        m(m),
+        m_limit(lim), 
+        m_conflict(false)
+    {}
+
+
+    grobner::~grobner() {
+        del_equations(0);
+    }
+
+    void grobner::compute_basis() {
+        while (!done() && compute_basis_step()) {
+            TRACE("grobner", display(tout););
+        }
+    }
+
+    bool grobner::compute_basis_step() {
+        m_stats.m_compute_steps++;
+        equation* e = pick_next();
+        if (!e) return false;
+        equation& eq = *e;
+        SASSERT(!eq.is_processed());
+        if (!simplify_using(eq, m_processed)) return false;
+        if (eliminate(eq)) return true;
+        if (check_conflict(eq)) return false;
+        if (!simplify_using(m_processed, eq)) return false;
+        TRACE("grobner", tout << "eq = "; display_equation(tout, eq););
+        superpose(eq);
+        toggle_processed(eq);
+        return simplify_using(m_to_simplify, eq);
+    }
+    
+    grobner::equation* grobner::pick_next() {
+        equation* eq = nullptr;
+        ptr_buffer<equation> to_delete;
+        for (auto* curr : m_to_simplify) {
+            if (!eq|| is_simpler(*curr, *eq)) {
+                eq = curr;
+            }
+        }
+        for (auto* e : to_delete)
+            del_equation(e);
+        if (eq) 
+            m_to_simplify.erase(eq);
+        return eq;
+    }
+
+    void grobner::superpose(equation const & eq) {
+        for (equation * target : m_processed) {
+            superpose(eq, *target);
+        }
+    }
+
+    bool grobner::simplify_using(equation& eq, equation_set const& eqs) {
+        bool simplified, changed_leading_term;
+        TRACE("grobner", tout << "simplifying: "; display_equation(tout, eq); tout << "using equalities of size " << eqs.size() << "\n";);
+        do {
+            simplified = false;
+            for (equation* p : eqs) {
+                if (simplify_source_target(*p, eq, changed_leading_term)) {
+                    simplified = true;
+                }
+                if (canceled() || eq.poly().is_val()) {
+                    break;
+                }
+            }
+        }
+        while (simplified && !eq.poly().is_val());
+
+        TRACE("grobner", tout << "simplification result: ";  display_equation(tout, eq););
+
+        check_conflict(eq);
+        return !done();
+    }
+
+    /*
+      Use the given equation to simplify equations in set
+    */
+    bool grobner::simplify_using(equation_set& set, equation const& eq) {
+        TRACE("grobner", tout << "poly " << eq.poly() <<  "\n";);
+        ptr_buffer<equation> to_delete;
+        ptr_buffer<equation> to_move;
+        bool changed_leading_term = false;
+        for (equation* target : set) {
+            if (canceled())
+                return false;
+            if (simplify_source_target(eq, *target, changed_leading_term)) {
+                if (is_trivial(*target))
+                    to_delete.push_back(target);
+                else if (check_conflict(*target))
+                    return false;
+                else if (changed_leading_term && target->is_processed()) {
+                    to_move.push_back(target);
+                }
+            }            
+        }
+        for (equation* eq : to_delete)
+            del_equation(eq);
+        for (equation* eq : to_move) 
+            toggle_processed(*eq);
+        return true;
+    }    
+    
+    // return true iff simplified
+    bool grobner::simplify_source_target(equation const& source, equation& target, bool& changed_leading_term) {
+        TRACE("grobner", tout << "simplifying: " << target.poly() << "\nusing: " << source.poly() << "\n";);
+        m_stats.m_simplified++;
+        pdd r = source.poly().reduce(target.poly());
+        if (r == source.poly()) {
+            return false;
+        }
+        changed_leading_term = target.is_processed() && m.different_leading_term(r, source.poly());
+        target = r;
+        target = m_dep_manager.mk_join(target.dep(), source.dep());
+        update_stats_max_degree_and_size(target);
+
+        TRACE("grobner", tout << "simplified " << target.poly() << "\n";);
+        return true;
+    }
+
+    // let eq1: ab+q=0, and eq2: ac+e=0, then qc - eb = 0
+    void grobner::superpose(equation const& eq1, equation const& eq2) {
+        TRACE("grobner_d", tout << "eq1="; display_equation(tout, eq1) << "eq2="; display_equation(tout, eq2););
+        pdd r(m);
+        if (!m.try_spoly(eq1.poly(), eq2.poly(), r)) return;
+        m_stats.m_superposed++;
+        if (r.is_zero()) return;
+        unsigned idx = m_equations.size();
+        equation* eq = alloc(equation, r, m_dep_manager.mk_join(eq1.dep(), eq2.dep()), idx);
+        update_stats_max_degree_and_size(*eq);
+        if (r.is_val()) set_conflict();
+        m_to_simplify.insert(eq);
+    }
+
+    void grobner::toggle_processed(equation& eq) {
+        if (eq.is_processed()) {
+            m_to_simplify.insert(&eq);
+            m_processed.erase(&eq);
+        }
+        else {
+            m_processed.insert(&eq);
+            m_to_simplify.erase(&eq);            
+        }
+        eq.set_processed(!eq.is_processed());        
+    }
+
+    grobner::equation_set const& grobner::equations() {
+        m_all_eqs.reset();
+        for (equation* eq : m_equations) if (eq) m_all_eqs.insert(eq);
+        return m_all_eqs;
+    }
+
+    void grobner::reset() {
+        del_equations(0);
+        SASSERT(m_equations.empty());    
+        m_processed.reset();
+        m_to_simplify.reset();
+        m_stats.reset();
+        m_conflict = false;
+    }
+
+    void grobner::add(pdd const& p, u_dependency * dep) {
+        if (p.is_zero()) return;
+        if (p.is_val()) set_conflict();
+        equation * eq = alloc(equation, p, dep, m_equations.size());
+        m_equations.push_back(eq);
+        m_to_simplify.insert(eq);
+    }   
+
+    void grobner::del_equations(unsigned old_size) {
+        for (unsigned i = m_equations.size(); i-- > old_size; ) {
+            del_equation(m_equations[i]);
+        }
+        m_equations.shrink(old_size);    
+    }
+
+    void grobner::del_equation(equation* eq) {
+        if (!eq) return;
+        unsigned idx = eq->idx();
+        m_processed.erase(eq);
+        m_to_simplify.erase(eq);
+        dealloc(m_equations[idx]);
+        m_equations[idx] = nullptr;
+    }
+    
+    bool grobner::canceled() {
+        return m_limit.get_cancel_flag();
+    }
+    
+    bool grobner::done() {
+        return m_to_simplify.size() + m_processed.size() >= m_config.m_eqs_threshold || canceled() || m_conflict;
+    }
+    
+    void grobner::update_stats_max_degree_and_size(const equation& e) {
+        // TBD m_stats.m_max_expr_size = std::max(m_stats.m_max_expr_size, e.poly().size());
+        // TBD m_stats.m_max_expr_degree = std::max(m_stats.m_max_expr_degree, e.poly.degree());
+    }
+    
+    std::ostream& grobner::print_stats(std::ostream & out) const {
+        return out << "stats:\nsteps = " << m_stats.m_compute_steps << "\nsimplified: " <<
+            m_stats.m_simplified << "\nsuperposed: " <<
+            m_stats.m_superposed << "\nexpr degree: " << m_stats.m_max_expr_degree <<
+            "\nexpr size: " << m_stats.m_max_expr_size << "\n";
+    }
+            
+    std::ostream& grobner::display_equation(std::ostream & out, const equation & eq) const {
+        out << "expr = " << eq.poly() << "\n";
+        m_print_dep(eq.dep(), out);
+        return out;
+    }
+
+    std::ostream& grobner::display(std::ostream& out) const {
+        out << "processed\n"; for (auto e : m_processed) display_equation(out, *e);
+        out << "to_simplify\n"; for (auto e : m_to_simplify) display_equation(out, *e);
+        return out;
+    }
+
+}
+

src/math/grobner/pdd_grobner.h

@@ -0,0 +1,126 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+#pragma once
+
+#include "math/dd/dd_pdd.h"
+#include "util/dependency.h"
+#include "util/obj_hashtable.h"
+#include "util/region.h"
+#include "util/rlimit.h"
+
+namespace dd {
+
+class grobner {
+public:
+    struct stats {
+        unsigned m_simplified;
+        unsigned m_max_expr_size;
+        unsigned m_max_expr_degree;
+        unsigned m_superposed;
+        unsigned m_compute_steps;
+        void reset() { memset(this, 0, sizeof(*this)); }
+    };
+
+    struct config {
+        unsigned m_eqs_threshold;
+    };
+
+private:
+    class equation {
+        bool                       m_processed;  //!< state
+        unsigned                   m_idx;        //!< position at m_equations
+        pdd                        m_poly;       //   simplified expressionted monomials
+        u_dependency *             m_dep;        //!< justification for the equality
+    public:
+        equation(pdd const& p, u_dependency* d, unsigned idx): 
+            m_processed(false),
+            m_idx(idx),
+            m_poly(p),
+            m_dep(d)
+        {}
+
+        const pdd& poly() const { return m_poly; }        
+        u_dependency * dep() const { return m_dep; }
+        unsigned hash() const { return m_idx; }
+        unsigned idx() const { return m_idx; }
+        void operator=(pdd& p) { m_poly = p; }
+        void operator=(u_dependency* d) { m_dep = d; }
+        bool is_processed() const { return m_processed; }
+        void set_processed(bool p) { m_processed = p; }
+    };
+
+    typedef obj_hashtable<equation> equation_set;
+    typedef ptr_vector<equation> equation_vector;
+    typedef std::function<void (u_dependency* d, std::ostream& out)> print_dep_t;
+
+    pdd_manager&                                 m;
+    reslimit&                                    m_limit;
+    stats                                        m_stats;
+    config                                       m_config;
+    print_dep_t                                  m_print_dep;
+    equation_vector                              m_equations;
+    equation_set                                 m_processed;
+    equation_set                                 m_to_simplify;
+    mutable u_dependency_manager                 m_dep_manager;
+    equation_set                                 m_all_eqs;
+    bool                                         m_conflict;
+public:
+    grobner(reslimit& lim, pdd_manager& m);
+    ~grobner();
+
+    void operator=(print_dep_t& pd) { m_print_dep = pd; }
+    void operator=(config const& c) { m_config = c; }
+
+    void reset();
+    void compute_basis();
+    void add(pdd const&, u_dependency * dep);
+    equation_set const& equations();
+    u_dependency_manager& dep() const { return m_dep_manager;  }
+    std::ostream& display_equation(std::ostream& out, const equation& eq) const;
+    std::ostream& display(std::ostream& out) const;
+
+private:
+    bool compute_basis_step();
+    equation* pick_next();
+    bool canceled();
+    bool done();
+    void superpose(equation const& eq1, equation const& eq2);
+    void superpose(equation const& eq);
+    bool simplify_source_target(equation const& source, equation& target, bool& changed_leading_term);
+    bool simplify_using(equation_set& set, equation const& eq);
+    bool simplify_using(equation& eq, equation_set const& eqs);
+    void toggle_processed(equation& eq);
+
+    bool eliminate(equation& eq) { return is_trivial(eq) && (del_equation(&eq), true); }
+    bool is_trivial(equation const& eq) const { return eq.poly().is_zero(); }    
+    bool is_simpler(equation const& eq1, equation const& eq2) { return eq1.poly() < eq2.poly(); }
+    bool check_conflict(equation const& eq) { return eq.poly().is_val() && !is_trivial(eq) && (set_conflict(), true); }    
+    void set_conflict() { m_conflict = true; }
+
+    void del_equations(unsigned old_size);
+    void del_equation(equation* eq);    
+
+    void update_stats_max_degree_and_size(const equation& e);
+
+    std::ostream& print_stats(std::ostream&) const;
+
+};
+
+}

src/util/dependency.h

@@ -321,6 +321,8 @@ public:
 // Implement old dependency manager used by interval and groebner 
 typedef scoped_dependency_manager<void*>             v_dependency_manager;
 typedef scoped_dependency_manager<void*>::dependency v_dependency;
+typedef scoped_dependency_manager<unsigned>             u_dependency_manager;
+typedef scoped_dependency_manager<unsigned>::dependency u_dependency;
 
 #endif /* DEPENDENCY_H_ */