#include "nlsat/levelwise.h"
#include "nlsat/nlsat_types.h"
#include "nlsat/nlsat_assignment.h"
#include <vector>
#include <unordered_map>
#include <map>
#include <set>
#include <algorithm>
#include <memory>
#include "math/polynomial/algebraic_numbers.h"
#include "nlsat_common.h"
#include <queue>
#include "math/polynomial/polynomial_cache.h"
#include "math/polynomial/polynomial.h"
#include "util/z3_exception.h"

bool is_set(unsigned k) { return static_cast<int>(k) != -1; }

namespace nlsat {
    enum relation_mode {
        biggest_cell,
        chain,
        lowest_degree
    };
    enum prop_enum {
        ir_ord,
        an_del,
        non_null,
        ord_inv,
        sgn_inv,
        connected,
        an_sub,
        sample_holds,
        repr,
        _count 
    };
    struct nullified_poly_exception {};
    struct levelwise::impl {
        // Utility: call fn for each distinct irreducible factor of poly
        template<typename Func>
        void for_each_distinct_factor(polynomial_ref& poly, Func&& fn) {
            polynomial::polynomial_ref_vector factors(m_pm);
            ::nlsat::factor(poly, m_cache, factors);
            for (unsigned i = 0; i < factors.size(); i++) {
                polynomial_ref pr(factors.get(i), m_pm);
                fn(pr);
            }
        }
  
        // todo: consider to key polynomials in a set by using m_pm.eq
        struct property {
            prop_enum m_prop_tag;
            polynomial_ref   m_poly;
            unsigned         m_root_index = -1; // index of the root function, if applicable; -1 means unspecified
            property(prop_enum pr, polynomial_ref const & pp, int si = -1) : m_prop_tag(pr), m_poly(pp), m_root_index(si) {}
            property(prop_enum pr, polynomial_ref const & pp, polynomial::manager& pm) : m_prop_tag(pr), m_poly(pp), m_root_index(-1) {}
            // have to pass polynomial::manager& to create a polynomial_ref even with a null object
            property(prop_enum pr, polynomial::manager& pm) : m_prop_tag(pr), m_poly(polynomial_ref(pm)), m_root_index(-1) {}
            
        };
        struct compare_prop_tags {
            bool operator()(const property& a, const property& b) const {
                return (int)a.m_prop_tag > (int)b.m_prop_tag; // ir_ord dequed first
            }
        };
        typedef std::priority_queue<property, std::vector<property>, compare_prop_tags> property_queue;

        // p_q_plus: Enhanced property queue that allows iteration without modification
        class p_q_plus {
        private:
            property_queue m_queue;
            std::vector<property> m_elements; // Maintains all enqueued elements for iteration
            
        public:
            p_q_plus() = default;
            
            // Push a property to the queue and maintain the elements vector
            void push(const property& p) {
                m_queue.push(p);
                m_elements.push_back(p);
            }

            size_t size() const { return m_queue.size(); }
            
            // Pop and return the top element
            property pop() {
                property p = m_queue.top();
                m_queue.pop();
                // Remove from elements vector
                auto it = std::find_if(m_elements.begin(), m_elements.end(),
                                       [&p](const property& elem) {
                                           return elem.m_prop_tag == p.m_prop_tag && 
                                               elem.m_poly == p.m_poly &&
                                               elem.m_root_index == p.m_root_index;
                                       });
                if (it != m_elements.end()) {
                    m_elements.erase(it);
                }
                return p;
            }
            
            // Check if queue is empty
            bool empty() const {
                return m_queue.empty();
            }
            
            // Get reference to underlying elements for iteration
            const std::vector<property>& elements() const {
                return m_elements;
            }
            
            // Iterator interface
            auto begin() const { return m_elements.begin(); }
            auto end() const { return m_elements.end(); }
            
            // Clear the queue and elements
            void clear() {
                while (!m_queue.empty()) {
                    m_queue.pop();
                }
                m_elements.clear();
            }
        };

        struct root_function {
            scoped_anum val;
            indexed_root_expr ire;
            root_function(anum_manager& am, polynomial_ref p, unsigned idx, anum const& v)
                : val(am), ire{ p, idx } { am.set(val, v); }
            root_function(root_function&& other) noexcept : val(other.val.m()), ire(other.ire) { val = other.val; }
            root_function(root_function const&) = delete;
            root_function& operator=(root_function const&) = delete;
            // allow move-assignment so we can sort a vector<root_function>
            root_function& operator=(root_function&& other) noexcept { val = other.val; ire = other.ire; return *this; }
        };
        solver&                      m_solver;
        polynomial_ref_vector const& m_P;
        unsigned                     m_n; // the max of max_var(p) of all the polynomials, the level of the conflict
        unsigned                     m_level;// the current level while refining the properties
        pmanager&                    m_pm;
        anum_manager&                m_am;
        std::vector<p_q_plus>        m_Q; // the set of properties to prove
        std::vector<property>        m_to_refine; 
        std::vector<root_function_interval> m_I; // intervals per level (indexed by variable/level)
        bool m_fail = false;
        /*
          Let us suppose that l = m_level, and j is such that m_rfunc[j](x_l) <= sample(x_l) and
          m_rfunc[j](x_l) reaches the maximum of such numbern. Let k is the symmetrical definition for 
          the upper bounds, and for simplicity assume that both j and k are defined.
          The paper does not address it formally, but the idea is that the relation is legal 
          if and only if for every p < j, there is a path (p, p1), (p1, p2), ...(p_n, j) inside of the relation,
          and the corresponding statement for the upper bound.
        */
        struct relation_E {
            std::vector<root_function> m_rfunc; // the root functions on a level
            std::vector<std::pair<unsigned, unsigned>> m_pairs; // of the relation
            bool empty() const { return m_rfunc.size() == 0 && m_pairs.size() == 0; }
            void clear() {
                m_pairs.clear();
                m_rfunc.clear();
                m_l_start = m_l_end = m_u_start = m_u_end = -1;
            }
            // the indices point te the m_rfunc vector
            unsigned m_l_start = -1;
            unsigned m_l_end = -1;
            unsigned m_u_start = -1;
            unsigned m_u_end = -1;
            void add_pair(unsigned j, unsigned k) { m_pairs.emplace_back(j, k);}
        };
        relation_E m_rel;
        relation_mode m_relation_mode = biggest_cell; // there are other choices as well
        assignment const & sample() const { return m_solver.sample();}
        assignment & sample() { return m_solver.sample(); }
        polynomial::cache &          m_cache; 
        polynomial_ref_vector        m_psc_tmp;
        // max_x plays the role of n in algorith 1 of the levelwise paper.

        impl(solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache & cache)
            : m_solver(solver), m_P(ps),  m_n(max_x),  m_pm(pm), m_am(am), m_cache(cache),  m_psc_tmp(m_pm) {
            TRACE(lws, tout  << "m_n:" << m_n << "\n";);
            m_I.reserve(m_n); // cannot just resize bcs of the absence of the default constructor of root_function_interval
            for (unsigned i = 0; i < m_n; ++i)
                m_I.emplace_back(m_pm);
            m_Q.resize(m_n + 1);
        }
        // end constructor


        // helper overload so callers can pass either raw poly* or polynomial_ref
        unsigned  max_var(poly* p) { return m_pm.max_var(p); }
        unsigned  max_var(polynomial_ref const & p) { return m_pm.max_var(p); }

        // Wrapper to use cached PSC chain computation
        void psc_chain(polynomial_ref & p, polynomial_ref & q, unsigned x, polynomial_ref_vector & result) {
            m_cache.psc_chain(p, q, x, result);
        }

        // Use the PSC chain of p and its derivative w.r.t x to obtain a discriminant candidate.
        polynomial_ref psc_discriminant(polynomial_ref const & p_in, unsigned x) {
            polynomial_ref p(p_in);
            polynomial_ref d(m_pm);
            d = derivative(p, x);
            polynomial_ref_vector & chain = m_psc_tmp;
            chain.reset();
            psc_chain(p, d, x, chain);
            polynomial_ref disc(m_pm);
            for (unsigned i = 0; i < chain.size(); ++i) {
                disc = polynomial_ref(chain.get(i), m_pm);
                if (!is_const(disc)) 
                     break;
            }
            return disc;
        }

        // Helper to print out m_Q
        std::ostream& display(std::ostream& out) {
            out << "{\n";
            unsigned level = 0;    
            for (auto &q: m_Q) {
                if (q.empty()) { level++; continue; }
                out << "level:" << level << "[\n";
                for (const auto& pr : q) {
                    display(out, pr) << "\n";
                }
                out << "]\n";   
                level ++; 
            }
            out << "}\n";
            return out;
        }

        bool is_square_free(poly* p) {
            return m_pm.is_square_free(p);            
        }

        struct level_t {
            unsigned val;
            explicit level_t(unsigned lvl) { val = lvl; }
        };

        bool find_in_Q(polynomial_ref & f, unsigned level) {
            if (level >= m_Q.size())
                return false;
            poly* fp = f.get();
            if (!fp)
                return false;
            unsigned const fid = m_pm.id(fp);
            for (auto const& pr : m_Q[level]) {
                poly* qp = pr.m_poly.get();
                if (qp && m_pm.id(qp) == fid)
                    return true;
            }
            return false;
        }

        /*
          prepare the initial properties:         
          scan input polynomials, add sgn_inv / an_del properties and collect polynomials at level m_n
        */ 
        void collect_top_level_properties(polynomial_ref_vector & ps_of_n_level) {
            for (unsigned i = 0; i < m_P.size(); ++i) {
                polynomial_ref pi(m_P[i], m_pm);
                for_each_distinct_factor(pi, [&](polynomial_ref& f) {
                    unsigned level = max_var(f);
                    TRACE(lws, tout << "before normalize f:"  << f << "\n";);
                    normalize(f);
                    TRACE(lws, tout << "after normalize f:"  << f << "\n";);
                    if (find_in_Q(f, level)) {
                        TRACE(lws, tout << "found in Q\n";);
                        return;
                    }
                    if (level < m_n)
                        mk_prop(prop_enum::sgn_inv, f, level_t(level));
                    else if (level == m_n) {
                        mk_prop(prop_enum::an_del, f, level_t(level));
                        ps_of_n_level.push_back(f);
                    }
                    else 
                        UNREACHABLE();
                });
            }
        }

        // isolate and collect algebraic roots for the given polynomials
        void collect_roots_for_ps( polynomial_ref_vector const & ps_of_n_level, std::vector<std::pair<scoped_anum, poly*>> & root_vals) {
            for (unsigned i = 0; i < ps_of_n_level.size(); i++)  {               
                scoped_anum_vector roots(m_am);
                polynomial_ref p(ps_of_n_level[i], m_pm);
                m_am.isolate_roots(p, undef_var_assignment(sample(), m_n), roots);
                TRACE(lws,  
                      ::nlsat::display(tout << "roots of ", m_solver, p) << ": ";
                      ::nlsat::display(tout, roots);
                    );
                unsigned num_roots = roots.size();
                for (unsigned k = 0; k < num_roots; ++k) {
                    scoped_anum v(m_am);
                    m_am.set(v, roots[k]);
                    root_vals.emplace_back(std::move(v), p);
                }                
            }
        }

        // Given collected roots (possibly unsorted), sort them, and add ord_inv(resultant(...))
        // for adjacent roots coming from different polynomials. Avoid adding the same unordered pair twice.
        // Returns false on failure (e.g. when encountering an ambiguous zero resultant).
        void add_adjacent_resultants(std::vector<std::pair<scoped_anum, poly*>> & root_vals) {
            if (root_vals.size() < 2) return;
            std::sort(root_vals.begin(), root_vals.end(), [&](auto const & a, auto const & b){ return m_am.lt(a.first, b.first); });
            std::set<std::pair<unsigned,unsigned>> added_pairs;
            TRACE(lws, 
                  tout << "root_vals:";
                  for (const auto& rv : root_vals) {
                      tout << " [";
                      m_am.display(tout, rv.first);
                      if (rv.second) {
                          tout << ", poly: ";
                          ::nlsat::display(tout, m_solver, polynomial_ref(rv.second, m_pm));
                      }
                      tout << "]";
                  }
                  tout << std::endl;
                );
            for (unsigned j = 0; j + 1 < root_vals.size(); ++j) {
                poly* p1 = root_vals[j].second;
                poly* p2 = root_vals[j+1].second;
                if (p1 == p2) continue; // delineability of p1 handled by an_del
                unsigned id1 = polynomial::manager::id(polynomial_ref(p1, m_pm));
                unsigned id2 = polynomial::manager::id(polynomial_ref(p2, m_pm));
                std::pair<unsigned,unsigned> key = id1 < id2 ? std::make_pair(id1, id2) : std::make_pair(id2, id1);
                if (added_pairs.find(key) != added_pairs.end())
                    continue;
                added_pairs.insert(key);
                add_ord_inv_resultant(p1, p2, m_n);
            }
        }

        /*
          See comment 7.1 Generating explanation for MCSAT.
          Return Q = { sgn_inv(p) | level(p) < m_n }
          ∪ { an_del(p)  | level(p) == m_n }
          ∪ { ord_inv(resultant(p_j,p_{j+1})) for adjacent root functions }.
        */
        void init_properties() {
            polynomial_ref_vector ps_of_n_level(m_pm);
            collect_top_level_properties(ps_of_n_level);
            std::vector<std::pair<scoped_anum, poly*>> root_vals;
            collect_roots_for_ps(ps_of_n_level, root_vals);
            add_adjacent_resultants(root_vals);
        }
        

        // Compute root function interval from sorted roots.
        void compute_interval_from_sorted_roots() {
            root_function_interval & I = m_I[m_level];
            // default: whole line sector (-inf, +inf)
            I.section = false;
            I.l = nullptr; I.u = nullptr; 
            if (m_rel.empty()) return;
            if (!sample().is_assigned(m_level)) return;
            anum const& y_val = sample().value(m_level);
            TRACE(lws, tout << "sample val:";  m_am.display_decimal(tout, y_val); tout << "\n";);
                        
            // find first index where roots[idx].val >= y_val
            const auto & rfs = m_rel.m_rfunc;
            unsigned idx = 0;
            while (idx < rfs.size() && m_am.compare(rfs[idx].val, y_val) < 0) {
                TRACE(lws, tout << "idx:" << idx << ", val:"; m_am.display_decimal(tout, rfs[idx].val); tout << "\n";);
                ++idx;
            }
            if (idx < rfs.size() && m_am.compare(rfs[idx].val, y_val) == 0) {
                TRACE(lws, tout << "exact match at idx:" << idx << ", it's a section\n";);
                auto const& ire = rfs[idx].ire;
                I.section = true;
                I.l = ire.p; I.l_index = ire.i;
                I.u = nullptr; I.u_index = -1; // the section is defined by the I.l
                TRACE(lws, tout << "section bound -> p:"; if (I.l) m_pm.display(tout, I.l); tout << ", index:" << I.l_index << "\n";);
                m_rel.m_l_start = m_rel.m_l_end = idx;
                while (++idx < rfs.size() && m_am.compare(rfs[idx].val, y_val) == 0) {
                    m_rel.m_l_end = idx;
                    TRACE(lws, tout << "idx:" << idx << ", val:"; m_am.display_decimal(tout, rfs[idx].val); tout << "\n";);
                }
                TRACE(lws, display_relation(tout););
                return;
            }

            // sector: lower bound is last root with val < y, upper bound is first root with val > y
            if (idx > 0) {
                // find start,end of equal-valued group for lower bound
                unsigned start = idx - 1;
                m_rel.m_l_end = start;
                while (start > 0 && m_am.compare(rfs[start-1].val, rfs[start].val) == 0) {
                    --start;
                    TRACE(lws, tout << "start:" << start << ", val:"; m_am.display_decimal(tout, rfs[start].val); tout << "\n";);
                }
                m_rel.m_l_start = start;
                auto const& ire = rfs[start].ire;
                I.l = ire.p; I.l_index = ire.i;
            } 
            if (idx < rfs.size()) {
                // find start, end of equal-valued group for upper bound
                unsigned start = idx;
                m_rel.m_u_start = idx;
                while (start + 1 < rfs.size() && m_am.compare(rfs[start].val, rfs[start + 1].val) == 0) {
                    ++start;
                    TRACE(lws, tout << "start:" << start << ", val:"; m_am.display_decimal(tout, rfs[start].val); tout << "\n";);
                }
                auto const& ire = rfs[start].ire;
                m_rel.m_u_end = start;
                I.u = ire.p; I.u_index = ire.i;
            }
            TRACE(lws, display_relation(tout) << std::endl;);
        }

        property pop(p_q_plus& q) {
            return q.pop();
        }

        //works on m_level
        void apply_property_rules(prop_enum prop_to_avoid) {
            auto& q = m_Q[m_level];
            while(!q.empty()) {
                property p = pop(q); // there is a choice here of what property to pop
                if (p.m_prop_tag == prop_to_avoid) {
                    q.push(p);
                    break;
                }
                apply_pre(p);
            }
        }

        // Part B of construct_interval: build (I, E, ≼) representation for level i
        void build_representation() {
            collect_E();
            // todo: this order needs to be abstracted: it does not have to be linear.
            // We need a boolean function E_rel(a, b)
            std::sort(m_rel.m_rfunc.begin(), m_rel.m_rfunc.end(), [&](root_function const& a, root_function const& b){
                return m_am.lt(a.val, b.val);
            });
            TRACE(lws,
                  if (m_rel.empty()) tout << "E is empty\n";
                  else { tout << "E:\n";
                      tout << "roots:\n";
                      for (const auto& rf: m_rel.m_rfunc) {
                          display(tout, rf) << "\n";
                      } 
                      tout << "pairs:\n";
                      for (unsigned kk = 0; kk < m_rel.m_pairs.size(); ++kk) {
                          auto pair = m_rel.m_pairs[kk];
                          display(tout, m_rel.m_rfunc[pair.first]) << "<<<" ; display(tout, m_rel.m_rfunc[pair.second])<< "\n"; 
                      }
                  });
            compute_interval_from_sorted_roots();
            fill_relation_pairs();
            TRACE(lws, display(tout << "m_I[" << m_level << "]:", m_I[m_level]) << std::endl;);
        }

        void fill_relation_with_biggest_cell_heuristic() {
            unsigned l = m_rel.m_l_end; 
            if (is_set(l)) 
                for (unsigned j = 0; j < l; j++) 
                    m_rel.add_pair(j, l);
            
            unsigned u = m_rel.m_u_start;
            if (is_set(u)) 
                for (unsigned j = u + 1; j < m_rel.m_rfunc.size(); j++)
                    m_rel.add_pair(u, j);

            if (is_set(l) && is_set(u)) {
                SASSERT(l + 1 == u);
                m_rel.add_pair(l, u);        
            }
        }

        void fill_relation_pairs() {
            if (m_relation_mode == biggest_cell) 
                fill_relation_with_biggest_cell_heuristic();
            else
                NOT_IMPLEMENTED_YET();
        }

        // Step 1a: collect E the set of root functions on m_level
        void collect_E() {
            TRACE(lws, tout << "enter\n";);
            m_rel.clear();

            for (auto const& q : m_Q[m_level]) {
                if (q.m_prop_tag != prop_enum::sgn_inv)
                    continue;
                polynomial_ref p(m_pm);
                p = q.m_poly;

                SASSERT(max_var(p) == m_level);
                scoped_anum_vector roots(m_am);
                m_am.isolate_roots(polynomial_ref(p, m_pm), undef_var_assignment(sample(), m_level), roots);
                
                unsigned num_roots = roots.size();
                TRACE(lws, ::nlsat::display(tout << "p:", m_solver, p) << ",";
                      tout << "roots (" << num_roots << "):";
                      for (unsigned kk = 0; kk < num_roots; ++kk) {
                          tout << " "; m_am.display(tout, roots[kk]);
                      }
                      tout << std::endl;
                );
                for (unsigned k = 0; k < num_roots; ++k)
                    m_rel.m_rfunc.emplace_back(m_am, p, k + 1, roots[k]);
            }
            TRACE(lws, tout << "exit\n";);
        }

        void normalize(polynomial_ref & p) {
            if (!m_solver.lws_norm())
                return;
            unsigned level = max_var(p);
            if (find_in_Q(p, level)) return;
            TRACE(lws_norm, tout << "p:" << p << std::endl;);

            p = ::primitive(p, level);
            p = m_pm.flip_sign_if_lm_neg(p);
            p = m_cache.mk_unique(p); 
            TRACE(lws_norm,  tout << "normalized p:" << p << "\n";);
        }

        bool same_polynomial_up_to_constant(poly* a, poly* b) {
            if (a == b || m_pm.id(a) == m_pm.id(b))
                return true;
            polynomial_ref pa(a, m_pm);
            polynomial_ref pb(b, m_pm);
            if (is_zero(pa) || is_zero(pb) || is_const(pa) || is_const(pb))
                return false;
            unsigned v = max_var(pa);
            if (v != max_var(pb))
                return false;
            polynomial_ref ca(m_pm), cb(m_pm);
            m_pm.primitive(pa, v, ca);
            m_pm.primitive(pb, v, cb);
            ca = m_pm.flip_sign_if_lm_neg(ca);
            cb = m_pm.flip_sign_if_lm_neg(cb);
            return m_pm.eq(ca, cb);
        }

        void add_ord_inv_resultant(poly* a, poly* b, unsigned x) {
            polynomial_ref pa(a, m_pm);
            polynomial_ref pb(b, m_pm);
            polynomial_ref r(m_pm);
            polynomial_ref_vector & chain = m_psc_tmp;
            chain.reset();
            psc_chain(pa, pb, x, chain);

            for (unsigned i = 0; i < chain.size(); ++i) {
                r = polynomial_ref(chain.get(i), m_pm);
                if (is_const(r))
                    continue;
                break;
            }

            TRACE(lws,
                      tout << "psc resultant wrt x" << x << "\n";
                      ::nlsat::display(tout << "a: ", m_solver, pa) << "\n";
                      ::nlsat::display(tout << "b: ", m_solver, pb) << "\n";
                      ::nlsat::display(tout << "r: ", m_solver, r) << "\n";);
            
            for_each_distinct_factor(r, [&](polynomial_ref& f) {
                TRACE(lws, tout << "f:" << f << "\n";);
                mk_prop(prop_enum::ord_inv, f);
            });
        }


        // add a property to m_Q if an equivalent one is not already present.
        // Equivalence: same m_prop_tag and same level; polynomials checked for equality
        // (polynomials are already in primitive form, so constant multipliers are normalized away).
        void add_to_Q_if_new(property pr, unsigned level) {            
            if (pr.m_poly)
                normalize(pr.m_poly);
            SASSERT(pr.m_prop_tag != ord_inv || (pr.m_poly && !is_zero(pr.m_poly)));
            for (auto const & q : m_Q[level]) {
                if (q.m_prop_tag != pr.m_prop_tag) continue;
                if (q.m_root_index != pr.m_root_index) continue;
                if ((q.m_poly && !pr.m_poly) || (!q.m_poly && pr.m_poly)) continue;
                if (!q.m_poly && !pr.m_poly) return;
                if (m_pm.id(q.m_poly.get()) != m_pm.id(pr.m_poly.get())) continue;
                TRACE(lws, display(tout << "matched q:", q) << std::endl;);
                return;
            }
            SASSERT(!pr.m_poly || is_square_free(pr.m_poly));
            m_Q[level].push(pr);
        }
        
        // construct_interval: compute representation for level i and apply post rules.
        // Returns false on failure.
        // works on m_level
        void construct_interval() {
            m_rel.clear();
            apply_property_rules(prop_enum::sgn_inv);

            TRACE(lws, display(tout << "m_Q:") << std::endl;);

            build_representation();
            SASSERT(invariant());
            apply_property_rules(prop_enum::_count);
        }
        // handle ord_inv(discriminant_{x_{i+1}}(p)) for an_del pre-processing
        void add_ord_inv_discriminant_for(const property& p) {
            unsigned x = max_var(p.m_poly);
            if (degree(p.m_poly,  x) < 2) return;
            polynomial::polynomial_ref disc(m_pm);
            disc = psc_discriminant(p.m_poly, x);
            TRACE(lws, display(tout << "p:", p) << "\n"; ::nlsat::display(tout << "discriminant by x" << max_var(p.m_poly)<< ": ", m_solver, disc) << "\n";);
            if (!is_const(disc)) {
                for_each_distinct_factor(disc, [&](polynomial::polynomial_ref f) {
                    mk_prop(prop_enum::ord_inv, f);
                });
            }
        }
        
        
        //  handle sgn_inv(leading_coefficient_{x_{i+1}}(p)) for an_del pre-processing
        void add_sgn_inv_leading_coeff_for(const property& p) {
            TRACE(lws, tout << "p:"; display(tout, p) << "\n";);
            poly * pp = p.m_poly.get();
            unsigned lvl = max_var(p.m_poly);
            unsigned deg = m_pm.degree(pp, lvl);
            polynomial_ref coeff(m_pm);
            coeff = m_pm.coeff(pp, lvl, deg);
            if (!is_const(coeff)) {
                for_each_distinct_factor(coeff, [&](polynomial_ref& f) {
                    normalize(f);
                    mk_prop(prop_enum::sgn_inv, f);
                });
            }
        }

        void apply_pre_an_del(const property& p) {
            unsigned p_lvl = max_var(p.m_poly);
            if (p_lvl > 0) {
                mk_prop(prop_enum::an_sub, level_t(p_lvl - 1));
                mk_prop(prop_enum::connected, level_t(p_lvl - 1));
            }
            mk_prop(prop_enum::non_null, p.m_poly);
            add_ord_inv_discriminant_for(p);
            add_sgn_inv_leading_coeff_for(p);
        }

        // Pre-processing for connected(i) (Rule 4.11)
        void apply_pre_connected(const property & p) {
            // Rule 4.11 special-case: if the connected property refers to level 0 there's nothing to refine
            // further; just remove the property from Q and return.
            if (m_level == 0) {
                TRACE(lws, tout << "apply_pre_connected: level 0 -> erasing connected property and returning" << std::endl;);
                return;
            }

            // p.level > 0
            // Rule 4.11 precondition: when processing connected(i) we must ensure the next lower level
            // has connected(i-1) and repr(I,s) available. 
            /*
              Let Q := connected(i − 1)(R) ∧ repr(I, s)(R).
              Q, I = (sector, l, u), l ̸= −∞, u ̸= ∞, l (≼t) u, add ir_ord(≼, s)
              Q, I = (sector, l, u), l= −∞ ∨ u = ∞ 
              Q, I = (section, b) ⊢ connected(i)(R)
            */
            if (m_level) {
                mk_prop(prop_enum::connected, level_t(m_level - 1));
                mk_prop(prop_enum::repr, level_t(m_level - 1));            
            }
            if (!have_representation()) 
                return; // no change since the cell representation is not available            

            const auto& I = m_I[m_level];
            TRACE(lws, display(tout << "interval m_I[" << m_level << "]\n", I) << "\n";);
            if (I.is_section()) return;
            SASSERT(I.is_sector());
            if (!I.l_inf() && !I.u_inf()) {
                if (m_level)
                    mk_prop(ir_ord, level_t(m_level - 1));
            }
        }

        // If a polynomial has a coefficient which is non-zero constant then non_null holds
        bool has_const_coeff(const polynomial_ref& p) {
            unsigned level = max_var(p);            
            unsigned deg = m_pm.degree(p, level);
            for (unsigned j = 0; j <= deg; ++j) {
                polynomial_ref coeff(m_pm);
                coeff = m_pm.coeff(p, level, j);
                TRACE(lws, tout << "coeff:" << coeff << "\n";);
                if(m_pm.nonzero_const_coeff(p, level, j)) 
                    return true;
            }
            return false;
        }
        
        void apply_pre_non_null(const property& p) {
            TRACE(lws, tout << "p:"; display(tout, p) << std::endl;);
            if (has_const_coeff(p.m_poly))
                return;
            if (!try_non_null_via_coeffs(p))
                fail();
            return;
            // fallback: apply the first subrule
            // TODO: choose between try_non_null_via_coeffs and apply_pre_non_null_fallback
            apply_pre_non_null_fallback(p);
        }


        // If some coefficient c_j of p is constant non-zero at the sample, or
        // if c_j evaluates non-zero at the sample and we already have sgn_inv(c_j) in m_Q,
        // then non_null(p) holds on the region represented by 'rs' (if provided).
        // Returns true if non_null was established and the property p was removed.
        // TODO: use cached non-null properties
        bool try_non_null_via_coeffs(const property& p) {
            unsigned level = max_var(p.m_poly);
            poly* pp = p.m_poly.get();
            unsigned deg = m_pm.degree(pp, level);
            for (unsigned j = 0; j <= deg; ++j) {
                polynomial_ref coeff(m_pm);
                coeff = m_pm.coeff(pp, level, j);
                if (sign(coeff, sample(), m_am) == 0)
                    continue;

                for_each_distinct_factor(coeff, [&](const polynomial::polynomial_ref & f) {
                    mk_prop(prop_enum::sgn_inv, f);
                });
                return true;
            }
            TRACE(lws, tout << "false\n";);
            return false;
        }

        // Helper for Rule 4.2, subrule 1: fallback when subrule 2 does not apply.
        // sample(s)(R), degx_{i+1} (p) > 1, disc(x_{i+1} (p)(s)) ̸= 0, sgn_inv(disc(x_{i+1} (p))(R) 
        void apply_pre_non_null_fallback(const property& p) {
            unsigned level = max_var(p.m_poly);

            poly * pp = p.m_poly.get();
            unsigned deg = m_pm.degree(pp, level);
            // fallback applies only for degree > 1
            if (deg <= 1) {
                fail();
                return;
            }
            // compute discriminant w.r.t. the variable at p.level
            polynomial_ref disc(m_pm);
            disc = psc_discriminant(p.m_poly, level);
            SASSERT (!is_zero(disc)); // p.m_poly is supposed to be square free
            TRACE(lws, ::nlsat::display(tout << "discriminant: ", m_solver, disc) << "\n";);

            // If discriminant is non-constant, add sign-invariance requirement for it
            if (!is_const(disc)) {
                for_each_distinct_factor(disc, [&](const polynomial::polynomial_ref & f) {
                    mk_prop(prop_enum::sgn_inv, f);
                });
            }

            // non_null is established by the discriminant being non-zero at the sample
        }

        // an_sub(R) iff R is an analitical manifold
        // Rule 4.7
        void apply_pre_an_sub(const property& p) {
            mk_prop(prop_enum::repr, level_t(m_level)) ;
            if (m_level > 0) 
                mk_prop(prop_enum::an_sub, level_t(m_level -1));
            // if level == 0 then an_sub holds - bcs an empty set is an analytical submanifold 
        }

        /*
          Rule 4.13 : the precondition holds by construction
          The property repr(I, s) holds on R if and only if I.l ∈ irExpr(I.l.p, s) (if I.l ̸= −∞), I.u ∈ irExpr(I.u.p, s)
          (if I.u ̸= ∞) respectively I.b ∈ irExpr(I.b.p, s) and one of the following holds:
          • I = (sector, l, u), dom(θl,s ) ∩ dom(θu,s ) ⊇ R ↓[i−1] and R = {(r, r′) | r ∈ R ↓[i−1], r′ ∈ (θl,s (r), θu,s (r))};
          or
          • I = (sector, −∞, u), dom(θu,s ) ⊇ R ↓[i−1] and R = {(r, r′) | r ∈ R ↓[i−1], r′ ∈ (−∞, θu,s (r))}; or
          • I = (sector, l, ∞), dom(θl,s ) ⊇ R ↓[i−1] and R = {(r, r′) | r ∈ R ↓[i−1], r′ ∈ (θl,s (r), ∞)}; or
          • I = (sector, −∞, ∞) and R = {(r, r′) | r ∈ R ↓[i−1], r′ ∈ R}; or
          • I = (section, b), dom(θb,s ) ⊇ R ↓[i−1] and R = {(r, r′) | r ∈ R ↓[i−1], r′
          = θb,s (r)},
        */
        void apply_pre_repr(const property& p) {
            const auto& I = m_I[m_level];
            TRACE(lws, display(tout << "interval m_I[" << m_level << "]\n", I) << "\n";);
            if (m_level) 
                mk_prop(sample_holds, level_t(m_level - 1));
            if (I.is_section()) {
                /*sample(s)(R), holds(I)(R), I = (section, b), an_del(b.p)(R) */
                mk_prop(an_del, I.l);
            } else {
                /* sample(s)(R), holds(I)(R), I = (sector, l, u), l= −∞ ∨ an_del(l.p)(R), u = ∞ ∨ an_del(u.p)(R*/
                SASSERT(I.is_sector());
                
                if (!I.l_inf())
                    mk_prop(an_del, I.l);
                if (!I.u_inf())
                    mk_prop(an_del, I.u);
            }                
        }

        void apply_pre_sample(const property& p) {
            if (m_level == 0)
                return;
            mk_prop(sample_holds, level_t(m_level - 1));
            mk_prop(prop_enum::repr, level_t(m_level - 1));
        }

        void mk_prop(prop_enum pe, level_t level) {
            SASSERT(is_set(level.val));
            add_to_Q_if_new(property(pe, m_pm), level.val);
        }

        void mk_prop(prop_enum pe, polynomial_ref poly) {            
            add_to_Q_if_new(property(pe, poly), max_var(poly));
        }

        void mk_prop(prop_enum pe, polynomial_ref poly, level_t level) {            
            add_to_Q_if_new(property(pe, poly), level.val);
        }
        

        void apply_pre_sgn_inv(const property& p) {
            SASSERT(is_square_free(p.m_poly));
            scoped_anum_vector roots(m_am);
            SASSERT(max_var(p.m_poly) == m_level);
            m_am.isolate_roots(p.m_poly, undef_var_assignment(sample(), m_level), roots);
            if (roots.size() == 0) {
                /*  Rule 4.6. Let i ∈ N, R ⊆ Ri, s ∈ R_{i−1}, and realRoots(p(s, xi )) = ∅. 
                    sample(s)(R), an_del(p)(R) ⊢ sgn_inv(p)(R) */
                if (m_level)
                    mk_prop(sample_holds, level_t(m_level - 1));
                mk_prop(prop_enum::an_del, p.m_poly);
                return;
            }
            // now we have some roots at s
            const auto &I = m_I[m_level];
            TRACE(lws, display(tout << "I:", I) << std::endl;);
            if (I.section) {
                /* Rule 4.8. Let i ∈ N>0 , R ⊆ Ri
                   , s ∈ R_{i−1}
                   , p ∈ Q[x1, . . . , xi ], level(p) = i, and I be a root function interval
                   of level i. Assume that p is irreducible, and I = (section, b).
                   Let Q := an_sub(i − 1)(R) ∧ connected(i − 1)(R) ∧ repr(I, s)(R) ∧ an_del(b.p)(R).
                   Q, b.p= p ⊢ sgn_inv(p)(R)
                   Q, b.p ̸= p, ord_inv(resxi (b.p, p))(R) ⊢ sgn_inv(p)(R)
                */ 
                if (m_level) { 
                    mk_prop(prop_enum::an_sub, level_t(m_level - 1));
                    mk_prop(prop_enum::connected, level_t(m_level - 1));
                    mk_prop(prop_enum::repr, level_t(m_level - 1)); 
                }
                mk_prop(prop_enum::an_del, m_I[m_level].l);
                if (I.l == p.m_poly.get()) {
                    // nothing is added
                } else {
                        add_ord_inv_resultant(I.l, p.m_poly.get(), m_level);
                }
            } else {
                /*
                  Rule 4.10. Let i ∈ N>0 , R ⊆ Ri
                  , s ∈ Ri−1
                  , p ∈ Q[x1, . . . , xi ], level(p) = i, I be a root function interval of
                  level i, ≼ be an indexed root ordering of level i, and ≼t be the reflexive and transitive closure of ≼.
                  We choose l, u such that either I = (sector, l, u) or (I = (section, b) for b = l = u).
                  Assume that p is irreducible, irExpr(p, s) ̸= ∅, ξ.p is irreducible for all ξ ∈ dom(≼), ≼ matches s,
                  and for all ξ ∈ irExpr(p, s) it holds either ξ ≼t l or u ≼t ξ.
                  sample(s)(R), repr(I, s)(R), ir_ord(≼, s)(R), an_del(p)(R) ⊢ sgn_inv(p)(R)
                */
                // todo - read the preconditions on p it needs to be diff
                SASSERT(precondition_on_sign_inv(p));
                if (m_level) {
                    mk_prop(sample_holds, level_t(m_level - 1)); 
                    mk_prop(repr, level_t(m_level - 1));
                }
                mk_prop(ir_ord, level_t(m_level));
                mk_prop(an_del, p.m_poly);
            }
        }


        /*Assume that p is irreducible, irExpr(p, s) ̸= ∅, ξ.p is irreducible for all ξ ∈ dom(≼), ≼ matches s,
          and for all ξ ∈ irExpr(p, s) it holds either ξ ≼t l or u ≼t ξ.*/
        bool precondition_on_sign_inv(const property &p) {
            SASSERT(is_square_free(p.m_poly));
            SASSERT(max_var(p.m_poly) == m_level);
            return true;
        }

        /*
          Rule 4.5. Let i ∈ N>0 , R ⊆ Ri
          , s ∈ Ri
          , and p ∈ Q[x1, . . . , xi ], level(p) = i. Assume that p is irreducible.
          p(s) ̸= 0, sample(s)(R), sgn_inv(p)(R) ⊢ ord_inv(p)(R)
          p(s)= 0, sample(s)(R), an_sub(i− 1)(R), connected(i)(R), sgn_inv(p)(R), an_del(p)(R) ⊢ ord_inv(p)(R)
        */
        void apply_pre_ord_inv(const property& p) {
            SASSERT(p.m_prop_tag == prop_enum::ord_inv && is_square_free(p.m_poly));
            unsigned level = max_var(p.m_poly);
            auto sign_on_sample = sign(p.m_poly, sample(), m_am);
            mk_prop(sample_holds, level_t(level));
            if (sign_on_sample) {
                mk_prop(prop_enum::sgn_inv, p.m_poly);
            } else { // sign is zero
                if (level)
                    mk_prop(prop_enum::an_sub, level_t(level - 1));
                mk_prop(prop_enum::connected, level_t(level));                
                mk_prop(prop_enum::sgn_inv, p.m_poly);
                mk_prop(prop_enum::an_del, p.m_poly);
            }
        }
        
        void apply_pre(const property& p) {
            TRACE(lws, tout << "apply_pre BEGIN m_Q:"; display(tout) << std::endl; 
                  display(tout << "pre p:", p) << std::endl;);
            switch (p.m_prop_tag) {
            case prop_enum::an_del:
                apply_pre_an_del(p);
                break;
            case prop_enum::connected:
                apply_pre_connected(p);
                break;
            case prop_enum::non_null:
                apply_pre_non_null(p);
                break;
            case prop_enum::an_sub:
                apply_pre_an_sub(p);
                break;
            case prop_enum::repr:
                apply_pre_repr(p);
                break;
            case sample_holds:
                apply_pre_sample(p);
                break;
            case prop_enum::sgn_inv:
                apply_pre_sgn_inv(p);
                break;
            case prop_enum::ord_inv:
                apply_pre_ord_inv(p);
                break;
            case prop_enum::ir_ord:
                apply_pre_ir_ord(p);
                break;
            default:
                NOT_IMPLEMENTED_YET();
                break;
            }
            TRACE(lws,  tout << "apply_pre END m_Q:"; display(tout) << std::endl;);
            SASSERT (invariant());
        }

        bool have_representation() const { return m_rel.empty() == false; }

        void fail() {
            throw nullified_poly_exception();;
        }

        void apply_pre_ir_ord(const property& p) {
            /*Rule 4.9. Let i ∈ N, R ⊆ Ri, s ∈ Ri, and ≼ be an indexed root ordering of level i + 1.
              Assume that ξ.p is irreducible for all ξ ∈ dom(≼), and that ≼ matches s.
              sample(s)(R), an_sub(i)(R), connected(i)(R), ∀ξ ∈ dom(≼). an_del(ξ.p)(R), ∀(ξ,ξ′) ∈≼. ord_inv(resx_{i+1} (ξ.p, ξ′.p))(R) ⊢ ir_ord(≼, s)(R)
            */
            if (m_level > 0) {
                mk_prop(sample_holds, level_t(m_level -1 ));
                mk_prop(an_sub, level_t(m_level - 1));
                mk_prop(connected, level_t(m_level - 1));
            }
            for (const auto & pair: m_rel.m_pairs) {
                poly *a = m_rel.m_rfunc[pair.first].ire.p;
                poly *b = m_rel.m_rfunc[pair.second].ire.p;
                SASSERT(max_var(a) == max_var(b) && max_var(b) == m_level) ;
                if (m_pm.id(a) != m_pm.id(b))
                    add_ord_inv_resultant(a, b, m_level);
            }
        }
        
        bool invariant() {
            for (unsigned i = 0; i < m_Q.size(); i++) {                
                bool level_is_ok = std::all_of(m_Q[i].begin(), m_Q[i].end(), [this, i](const property& p){
                    if (!p.m_poly) return true;
                    if (max_var(p.m_poly) != i) {
                        enable_trace("lws");
                        display(std::cout << "max_var != i in p:", p) << std::endl;
                        return false;
                    }
                    if (! m_cache.contains(p.m_poly)) {
                        enable_trace("lws");
                        display(std::cout << "p.m_poly is not normalized:", p) << std::endl;
                        return false;
                    }
                    return true;
                });

                if (! level_is_ok)
                    return false;
            }
            return true;
        }

        
        // return an empty vector on failure, otherwise returns the cell representations with intervals
        std::vector<root_function_interval> single_cell_work() {
            TRACE(lws, m_solver.display_assignment(tout << "sample():\n") << std::endl; ::nlsat::display(tout << "m_P:\n", m_solver, m_P) << "\n";);
            if (m_n == 0) return m_I; // we have an empty sample
            m_level = m_n;
            
            init_properties(); // initializes m_Q as a queue of properties on levels <= m_n
            SASSERT(m_rel.empty());
            apply_property_rules(prop_enum::_count); // reduce the level from m_n to m_n - 1 to be consumed by construct_interval
            SASSERT(m_Q[m_n].size() == 0);
            SASSERT(m_level == m_n); 
            do {  // m_level changes from m_n - 1 to 0
                m_level--;
                construct_interval();
            } while (m_level != 0);
            return m_I;
        }
        std::vector<root_function_interval> single_cell() {
            try {
                single_cell_work();
            } catch (nullified_poly_exception &) {
                m_fail = true;
            }
            return m_I;
        }
        
        // Pretty-print helpers
        static const char* prop_name(prop_enum p) {
            switch (p) {
            case prop_enum::ir_ord:               return "ir_ord";
            case prop_enum::an_del:               return "an_del";
            case prop_enum::non_null:             return "non_null";
            case prop_enum::ord_inv:              return "ord_inv";
            case prop_enum::sgn_inv:              return "sgn_inv";
            case prop_enum::connected:            return "connected";
            case prop_enum::an_sub:               return "an_sub";
            case sample_holds:                    return "sample";
            case prop_enum::repr:                 return "repr";
            case prop_enum::_count:               return "_count";
            }
            return "?";
        }

        std::ostream& display(std::ostream& out, std::vector<property> & prs) const {
            for (const auto &pr : prs) {
                display(out, pr) << std::endl;
            }
            return out;
        }
 
        std::ostream& display_relation(std::ostream& out) const {
            out << "Relation_E:\n";
            if (m_rel.empty()) {
                out << "  (empty)\n";
                return out;
            }
            out << "  Root functions:\n";
            for (size_t i = 0; i < m_rel.m_rfunc.size(); ++i) {
                out << "    [" << i << "]: ";
                display(out, m_rel.m_rfunc[i], true) << "\n";
            }
            out << "  Pairs:\n";
            for (const auto& pair : m_rel.m_pairs) {
                out << "    (" << pair.first << ", " << pair.second << ")\n";
            }
            out << "  Indices:\n";
            out << "    m_l_start:" << (int)m_rel.m_l_start << ", m_l_end:" << (int)m_rel.m_l_end << "\n";
            out << "    m_u_start:" << (int)m_rel.m_u_start << ", m_u_end:" << (int)m_rel.m_u_end << "\n";
            return out;
        }

        std::ostream& display(std::ostream& out, const property & pr) const {
            out << "{prop:" << prop_name(pr.m_prop_tag);
            if (is_set(pr.m_root_index)) out   << ", m_root_index:" << pr.m_root_index;
            if (pr.m_poly) {
                out << ", poly:";
                ::nlsat::display(out, m_solver, pr.m_poly);
            }
            else {
                out << ", p:null";
            }
            out << "}";
            return out;
        }
       
        std::ostream& display(std::ostream& out, const indexed_root_expr& ire ) const {
            out << "RootExpr: p:";
            ::nlsat::display(out, m_solver, ire.p);
            out << ", root_index:" << ire.i;
            return out;
        }

        // Print the indexed root function's value. If print_approx is true print a decimal
        // approximation, otherwise print the full representation.
        std::ostream& display(std::ostream& out, const root_function& f, bool print_approx = true ) const {
            display(out << "indexed_root_function:", f.ire) << "\n" << "val:";
            if (print_approx)
                m_am.display_decimal(out, f.val);
            else
                m_am.display(out, f.val);
            return out;
        }
        
        std::ostream& display(std::ostream& out, const root_function_interval& I) const {
            return ::nlsat::display(out, m_solver, I);
        }
    };
// constructor
    levelwise::levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var n, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache& cache)
        : m_impl(new impl(solver, ps, n, s, pm, am, cache)) {}

    levelwise::~levelwise() { delete m_impl; }

    std::vector<levelwise::root_function_interval> levelwise::single_cell() {
        return m_impl->single_cell();
    }

    bool levelwise::failed() const { return m_impl->m_fail; }
} // namespace nlsat

// Free pretty-printer for symbolic_interval
std::ostream& nlsat::display(std::ostream& out, solver& s, levelwise::root_function_interval const& I) {
    if (I.section) {
        out << "Section: ";
        if (I.l == nullptr)
            out << "(undef)";
        else {
            ::nlsat::display(out, s, I.l);
            out << "[root_index:" << I.l_index << "]";
        }
    }
    else {
        out << "Sector: (";
        if (I.l_inf())
            out << "-oo";
        else {
            ::nlsat::display(out, s, I.l);
            out << "[root_index:" << I.l_index << "]";
        }
        out << ", ";
        if (I.u_inf())
            out << "+oo";
        else {
            ::nlsat::display(out, s, I.u);
            out << "[root_index:" << I.u_index << "]";
        }
        out << ")";
    }
    return out;
}