Reorganizing the code. Moved nlsat to its own directory.

Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>

src/nlsat/nlsat_assignment.h

@@ -0,0 +1,78 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_assignment.h
+
+Abstract:
+
+    Assignment: Var -> Algebraic Number 
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-08.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_ASSIGNMENT_H_
+#define _NLSAT_ASSIGNMENT_H_
+
+#include"nlsat_types.h"
+#include"algebraic_numbers.h"
+
+namespace nlsat {
+
+    /**
+       \brief A mapping from variables to values. 
+       This mapping is used to encode the current partial interpretation in nlsat.
+    */
+    class assignment : public polynomial::var2anum {
+        scoped_anum_vector m_values;
+        svector<bool>      m_assigned;
+    public:
+        assignment(anum_manager & _m):m_values(_m) {}
+        virtual ~assignment() {}
+        anum_manager & am() const { return m_values.m(); }
+        void swap(assignment & other) {
+            m_values.swap(other.m_values);
+            m_assigned.swap(other.m_assigned);
+        }
+        void set_core(var x, anum & v) {
+            m_values.reserve(x+1, anum());
+            m_assigned.reserve(x+1, false); 
+            m_assigned[x] = true;
+            am().swap(m_values[x], v); 
+        }
+        void set(var x, anum const & v) {
+            m_values.reserve(x+1, anum());
+            m_assigned.reserve(x+1, false); 
+            m_assigned[x] = true;
+            am().set(m_values[x], v); 
+        }
+        void reset(var x) { if (x < m_assigned.size()) m_assigned[x] = false; }
+        bool is_assigned(var x) const { return m_assigned.get(x, false); }
+        anum const & value(var x) const { return m_values[x]; }
+        virtual anum_manager & m() const { return am(); }
+        virtual bool contains(var x) const { return is_assigned(x); }
+        virtual anum const & operator()(var x) const { SASSERT(is_assigned(x)); return value(x); }
+    };
+    
+    /**
+       \brief Wrapper for temporarily unassigning a given variable y.
+       That is, given an assignment M, M' = undef_var_assignment(M, y) is identical
+       to M, but y is unassigned in M'
+    */
+    class undef_var_assignment : public polynomial::var2anum {
+        assignment const & m_assignment;
+        var                m_y;
+    public:
+        undef_var_assignment(assignment const & a, var y):m_assignment(a), m_y(y) {}
+        virtual anum_manager & m() const { return m_assignment.am(); }
+        virtual bool contains(var x) const { return x != m_y && m_assignment.is_assigned(x); }
+        virtual anum const & operator()(var x) const { return m_assignment.value(x); }
+    };
+};
+
+#endif

src/nlsat/nlsat_clause.cpp

@@ -0,0 +1,49 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_clause.cpp
+
+Abstract:
+
+    Clauses used in the Nonlinear arithmetic satisfiability procedure. 
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-02.
+
+Revision History:
+
+--*/
+#include"nlsat_clause.h"
+
+namespace nlsat {
+
+    clause::clause(unsigned id, unsigned sz, literal const * lits, bool learned, assumption_set as):
+        m_id(id),
+        m_size(sz),
+        m_capacity(sz),
+        m_learned(learned),
+        m_activity(0),
+        m_assumptions(as) {
+        for (unsigned i = 0; i < sz; i++) {
+            m_lits[i] = lits[i];
+        }
+    }
+
+    bool clause::contains(literal l) const {
+        for (unsigned i = 0; i < m_size; i++)
+            if (m_lits[i] == l)
+                return true;
+        return false;
+    }
+
+    bool clause::contains(bool_var v) const {
+        for (unsigned i = 0; i < m_size; i++)
+            if (m_lits[i].var() == v)
+                return true;
+        return false;
+    }
+
+};

src/nlsat/nlsat_clause.h

@@ -0,0 +1,61 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_clause.h
+
+Abstract:
+
+    Clauses used in the Nonlinear arithmetic satisfiability procedure. 
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-02.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_CLAUSE_H_
+#define _NLSAT_CLAUSE_H_
+
+#include"nlsat_types.h"
+#include"vector.h"
+
+namespace nlsat {
+
+    class clause {
+        friend class solver;
+        unsigned         m_id;
+        unsigned         m_size;
+        unsigned         m_capacity:31;
+        unsigned         m_learned:1;
+        unsigned         m_activity;
+        assumption_set   m_assumptions;
+        literal          m_lits[0];
+        static size_t get_obj_size(unsigned num_lits) { return sizeof(clause) + num_lits * sizeof(literal); }
+        size_t get_size() const { return get_obj_size(m_capacity); }
+        clause(unsigned id, unsigned sz, literal const * lits, bool learned, assumption_set as);
+    public:
+        unsigned size() const { return m_size; }
+        unsigned id() const { return m_id; }
+        literal & operator[](unsigned idx) { SASSERT(idx < m_size); return m_lits[idx]; }
+        literal const & operator[](unsigned idx) const { SASSERT(idx < m_size); return m_lits[idx]; }
+        bool is_learned() const { return m_learned; }
+        literal * begin() { return m_lits; }
+        literal * end() { return m_lits + m_size; }
+        literal const * c_ptr() const { return m_lits; }
+        void inc_activity() { m_activity++; }
+        void set_activity(unsigned v) { m_activity = v; }
+        unsigned get_activity() const { return m_activity; }
+        bool contains(literal l) const;
+        bool contains(bool_var v) const;
+        void shrink(unsigned num_lits) { SASSERT(num_lits <= m_size); if (num_lits < m_size) { m_size = num_lits; } }
+        assumption_set assumptions() const { return m_assumptions; }
+    };
+
+    typedef ptr_vector<clause> clause_vector;
+
+};
+
+#endif

src/nlsat/nlsat_evaluator.cpp

@@ -0,0 +1,679 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_evaluator.cpp
+
+Abstract:
+
+    Helper class for computing the infeasible intervals of an
+    arithmetic literal.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-12.
+
+Revision History:
+
+--*/
+#include"nlsat_evaluator.h"
+
+namespace nlsat {
+
+    struct evaluator::imp {
+        assignment const &       m_assignment;
+        pmanager &               m_pm;
+        small_object_allocator & m_allocator;
+        anum_manager &           m_am;
+        interval_set_manager     m_ism;
+        scoped_anum_vector       m_tmp_values;
+        scoped_anum_vector       m_add_roots_tmp;
+        scoped_anum_vector       m_inf_tmp;
+        
+        // sign tables: light version
+        struct sign_table {
+            anum_manager &     m_am;
+            struct section {
+                anum     m_root;
+                unsigned m_pos;
+            };
+            svector<section>   m_sections;
+            unsigned_vector    m_sorted_sections; // refs to m_sections
+            unsigned_vector    m_poly_sections;   // refs to m_sections
+            svector<int>       m_poly_signs;
+            struct poly_info {
+                unsigned       m_num_roots;
+                unsigned       m_first_section;   // idx in m_poly_sections;
+                unsigned       m_first_sign;      // idx in m_poly_signs;
+                poly_info(unsigned num, unsigned first_section, unsigned first_sign):
+                    m_num_roots(num),
+                    m_first_section(first_section),
+                    m_first_sign(first_sign) {
+                }
+            };
+            svector<poly_info> m_info;
+            
+            sign_table(anum_manager & am):m_am(am) {}
+
+            ~sign_table() {
+                reset();
+            }
+
+            void reset() {
+                unsigned sz = m_sections.size();
+                for (unsigned i = 0; i < sz; i++)
+                    m_am.del(m_sections[i].m_root);
+                m_sections.reset();
+                m_sorted_sections.reset();
+                m_poly_sections.reset();
+                m_poly_signs.reset();
+                m_info.reset();
+            }
+
+            unsigned mk_section(anum & v, unsigned pos) {
+                unsigned new_id = m_sections.size();
+                m_sections.push_back(section());
+                section & new_s = m_sections.back();
+                m_am.set(new_s.m_root, v);
+                new_s.m_pos = pos;
+                return new_id;
+            }
+
+            // Merge the new roots of a polynomial p into m_sections & m_sorted_sections.
+            // Store the section ids for the new roots in p_section_ids
+            unsigned_vector new_sorted_sections;
+            void merge(anum_vector & roots, unsigned_vector & p_section_ids) {
+                new_sorted_sections.reset();  // new m_sorted_sections
+                unsigned i1  = 0;
+                unsigned sz1 = m_sorted_sections.size();
+                unsigned i2  = 0;
+                unsigned sz2 = roots.size();
+                unsigned j   = 0;
+                while (i1 < sz1 && i2 < sz2) {
+                    unsigned s1_id = m_sorted_sections[i1];
+                    section & s1   = m_sections[s1_id];
+                    anum &    r2   = roots[i2];
+                    int c = m_am.compare(s1.m_root, r2);
+                    if (c == 0) {
+                        new_sorted_sections.push_back(s1_id);
+                        p_section_ids.push_back(s1_id);
+                        s1.m_pos = j;
+                        i1++;
+                        i2++;
+                    }
+                    else if (c < 0) {
+                        new_sorted_sections.push_back(s1_id);
+                        s1.m_pos = j;
+                        i1++;
+                    }
+                    else {
+                        // create new section
+                        unsigned new_id = mk_section(r2, j);
+
+                        //
+                        new_sorted_sections.push_back(new_id);
+                        p_section_ids.push_back(new_id);
+                        i2++;
+                    }
+                    j++;
+                }
+                SASSERT(i1 == sz1 || i2 == sz2);
+                while (i1 < sz1) {
+                    unsigned s1_id = m_sorted_sections[i1];
+                    section & s1   = m_sections[s1_id];
+                    
+                    new_sorted_sections.push_back(s1_id);
+                    s1.m_pos = j;
+                    
+                    i1++;
+                    j++;
+                }
+                while (i2 < sz2) {
+                    anum &    r2   = roots[i2];
+                    // create new section
+                    unsigned new_id = mk_section(r2, j);
+                    
+                    //
+                    new_sorted_sections.push_back(new_id);
+                    p_section_ids.push_back(new_id);
+                    i2++;
+                    j++;
+                }
+                m_sorted_sections.swap(new_sorted_sections);
+            }
+
+            /**
+               \brief Add polynomial with the given roots and signs.
+            */
+            unsigned_vector p_section_ids;
+            void add(anum_vector & roots, svector<int> & signs) {
+                p_section_ids.reset();
+                if (!roots.empty())
+                    merge(roots, p_section_ids);
+                unsigned first_sign    = m_poly_signs.size();
+                unsigned first_section = m_poly_sections.size();
+                unsigned num_signs = signs.size();
+                // Must normalize signs since we use arithmetic operations such as *
+                // during evaluation.
+                // Without normalization, overflows may happen, and wrong results may be produced.
+                for (unsigned i = 0; i < num_signs; i++)
+                    m_poly_signs.push_back(normalize_sign(signs[i]));
+                m_poly_sections.append(p_section_ids);
+                m_info.push_back(poly_info(roots.size(), first_section, first_sign));
+                SASSERT(check_invariant());
+            }
+
+            /**
+               \brief Add constant polynomial 
+            */
+            void add_const(int sign) {
+                unsigned first_sign    = m_poly_signs.size();
+                unsigned first_section = m_poly_sections.size();
+                m_poly_signs.push_back(normalize_sign(sign));
+                m_info.push_back(poly_info(0, first_section, first_sign));
+            }
+
+            unsigned num_cells() const {
+                return m_sections.size() * 2 + 1;
+            }
+            
+            /**
+               \brief Return true if the given cell is a section (i.e., root)
+            */
+            bool is_section(unsigned c) const {
+                SASSERT(c < num_cells());
+                return c % 2 == 1;
+            }
+
+            /**
+               \brief Return true if the given cell is a sector (i.e., an interval between roots, or (-oo, min-root), or (max-root, +oo)).
+            */
+            bool is_sector(unsigned c) const {
+                SASSERT(c < num_cells());
+                return c % 2 == 0;
+            }
+
+            /**
+               \brief Return the root id associated with the given section.
+               
+               \pre is_section(c)
+            */
+            unsigned get_root_id(unsigned c) const {
+                SASSERT(is_section(c));
+                return m_sorted_sections[c/2];
+            }
+            
+            // Return value of root idx.
+            // If idx == UINT_MAX, it return zero (this is a hack to simplify the infeasible_intervals method
+            anum const & get_root(unsigned idx) const {
+                static anum zero;
+                if (idx == UINT_MAX)
+                    return zero;
+                SASSERT(idx < m_sections.size());
+                return m_sections[idx].m_root;
+            }
+
+            static unsigned section_id_to_cell_id(unsigned section_id) {
+                return section_id*2 + 1;
+            }
+            
+            // Return the cell_id of the root i of pinfo
+            unsigned cell_id(poly_info const & pinfo, unsigned i) const {
+                return section_id_to_cell_id(m_sections[m_poly_sections[pinfo.m_first_section + i]].m_pos);
+            }
+            
+            // Return the sign idx of pinfo
+            int sign(poly_info const & pinfo, unsigned i) const {
+                return m_poly_signs[pinfo.m_first_sign + i];
+            }
+            
+#define LINEAR_SEARCH_THRESHOLD 8
+            int sign_at(unsigned info_id, unsigned c) const {
+                poly_info const & pinfo  = m_info[info_id];
+                unsigned num_roots = pinfo.m_num_roots;
+                if (num_roots < LINEAR_SEARCH_THRESHOLD) {
+                    unsigned i = 0;
+                    for (; i < num_roots; i++) {
+                        unsigned section_cell_id = cell_id(pinfo, i);
+                        if (section_cell_id == c)
+                            return 0;
+                        else if (section_cell_id > c)
+                            break;
+                    }
+                    return sign(pinfo, i);
+                }
+                else {
+                    if (num_roots == 0)
+                        return sign(pinfo, 0);
+                    unsigned root_1_cell_id = cell_id(pinfo, 0);
+                    unsigned root_n_cell_id = cell_id(pinfo, num_roots - 1);
+                    if (c < root_1_cell_id)
+                        return sign(pinfo, 0);
+                    else if (c == root_1_cell_id || c == root_n_cell_id)
+                        return 0;
+                    else if (c > root_n_cell_id)
+                        return sign(pinfo, num_roots);
+                    int low  = 0;
+                    int high = num_roots-1;
+                    while (true) {
+                        SASSERT(0 <= low && high < static_cast<int>(num_roots));
+                        SASSERT(cell_id(pinfo, low) < c);
+                        SASSERT(c < cell_id(pinfo, high));
+                        if (high == low + 1) {
+                            SASSERT(cell_id(pinfo, low) < c);
+                            SASSERT(c < cell_id(pinfo, low+1));
+                            return sign(pinfo, low+1);
+                        }
+                        SASSERT(high > low + 1);
+                        int mid   = low + ((high - low)/2);
+                        SASSERT(low < mid && mid < high);
+                        unsigned mid_cell_id = cell_id(pinfo, mid);
+                        if (mid_cell_id == c) {
+                            return 0;
+                        }
+                        if (c < mid_cell_id) {
+                            high = mid;
+                        }
+                        else {
+                            SASSERT(mid_cell_id < c);
+                            low  = mid;
+                        }
+                    }
+                }
+            }
+            
+            bool check_invariant() const {
+                SASSERT(m_sections.size() == m_sorted_sections.size());
+                for (unsigned i = 0; i < m_sorted_sections.size(); i++) {
+                    SASSERT(m_sorted_sections[i] < m_sections.size());
+                    SASSERT(m_sections[m_sorted_sections[i]].m_pos == i);
+                }
+                unsigned total_num_sections = 0;
+                unsigned total_num_signs = 0;
+                for (unsigned i = 0; i < m_info.size(); i++) {
+                    SASSERT(m_info[i].m_first_section <= m_poly_sections.size());
+                    SASSERT(m_info[i].m_num_roots == 0 || m_info[i].m_first_section < m_poly_sections.size());
+                    SASSERT(m_info[i].m_first_sign < m_poly_signs.size());
+                    total_num_sections += m_info[i].m_num_roots;
+                    total_num_signs += m_info[i].m_num_roots + 1;
+                }
+                SASSERT(total_num_sections == m_poly_sections.size());
+                SASSERT(total_num_signs == m_poly_signs.size());
+                return true;
+            }
+
+            // Display sign table for the given variable
+            void display(std::ostream & out) const {
+                out << "sections:\n  ";
+                for (unsigned i = 0; i < m_sections.size(); i++) {
+                    if (i > 0) out << " < ";
+                    m_am.display_decimal(out, m_sections[m_sorted_sections[i]].m_root);
+                }
+                out << "\n";
+                out << "sign variations:\n";
+                for (unsigned i = 0; i < m_info.size(); i++) {
+                    out << "  ";
+                    for (unsigned j = 0; j < num_cells(); j++) {
+                        if (j > 0)
+                            out << " ";
+                        int s = sign_at(i, j);
+                        if (s < 0) out << "-";
+                        else if (s == 0) out << "0";
+                        else out << "+";
+                    }
+                    out << "\n";
+                }
+            }
+
+            // Display sign table for the given variable
+            void display_raw(std::ostream & out) const {
+                out << "sections:\n  ";
+                for (unsigned i = 0; i < m_sections.size(); i++) {
+                    if (i > 0) out << " < ";
+                    m_am.display_decimal(out, m_sections[m_sorted_sections[i]].m_root);
+                }
+                out << "\n";
+                out << "poly_info:\n";
+                for (unsigned i = 0; i < m_info.size(); i++) {
+                    out << "  roots:";
+                    poly_info const & info = m_info[i];
+                    for (unsigned j = 0; j < info.m_num_roots; j++) {
+                        out << " ";
+                        out << m_poly_sections[info.m_first_section+j];
+                    }
+                    out << ", signs:";
+                    for (unsigned j = 0; j < info.m_num_roots+1; j++) {
+                        out << " ";
+                        int s = m_poly_signs[info.m_first_sign+j];
+                        if (s < 0) out << "-";
+                        else if (s == 0) out << "0";
+                        else out << "+";
+                    }
+                    out << "\n";
+                }
+            }
+        };
+
+        sign_table m_sign_table_tmp;
+
+        imp(assignment const & x2v, pmanager & pm, small_object_allocator & allocator):
+            m_assignment(x2v),
+            m_pm(pm),
+            m_allocator(allocator),
+            m_am(m_assignment.am()),
+            m_ism(m_am, allocator),
+            m_tmp_values(m_am),
+            m_add_roots_tmp(m_am),
+            m_inf_tmp(m_am),
+            m_sign_table_tmp(m_am) {
+        }
+
+        var max_var(poly const * p) const {
+            return m_pm.max_var(p);
+        }
+
+        /**
+           \brief Return the sign of the polynomial in the current interpration.
+           
+           \pre All variables of p are assigned in the current interpretation.
+        */
+        int eval_sign(poly * p) {
+            // TODO: check if it is useful to cache results
+            SASSERT(m_assignment.is_assigned(max_var(p)));
+            return m_am.eval_sign_at(polynomial_ref(p, m_pm), m_assignment);
+        }
+        
+        bool satisfied(int sign, atom::kind k) {
+            return 
+                (sign == 0 && (k == atom::EQ || k == atom::ROOT_EQ || k == atom::ROOT_LE || k == atom::ROOT_GE)) ||
+                (sign <  0 && (k == atom::LT || k == atom::ROOT_LT || k == atom::ROOT_LE)) ||
+                (sign >  0 && (k == atom::GT || k == atom::ROOT_GT || k == atom::ROOT_GE));
+        }
+
+        bool satisfied(int sign, atom::kind k, bool neg) {
+            bool r = satisfied(sign, k);
+            return neg ? !r : r;
+        }
+
+        bool eval_ineq(ineq_atom * a, bool neg) {
+            SASSERT(m_assignment.is_assigned(a->max_var()));
+            // all variables of a were already assigned... 
+            atom::kind k = a->get_kind();
+            unsigned sz  = a->size();
+            int sign = 1;
+            for (unsigned i = 0; i < sz; i++) {
+                int curr_sign = eval_sign(a->p(i));
+                if (a->is_even(i) && curr_sign < 0)
+                    curr_sign = 1;
+                sign = sign * curr_sign;
+                if (sign == 0)
+                    break;
+            }
+            return satisfied(sign, k, neg);
+        }
+
+        bool eval_root(root_atom * a, bool neg) {
+            SASSERT(m_assignment.is_assigned(a->max_var()));
+            // all variables of a were already assigned... 
+            atom::kind k = a->get_kind();
+            scoped_anum_vector & roots = m_tmp_values;
+            roots.reset();
+            m_am.isolate_roots(polynomial_ref(a->p(), m_pm), undef_var_assignment(m_assignment, a->x()), roots);
+            SASSERT(a->i() > 0);
+            if (a->i() > roots.size())
+                return false; // p does have sufficient roots
+            int sign = m_am.compare(m_assignment.value(a->x()), roots[a->i() - 1]);
+            return satisfied(sign, k, neg);
+        }
+        
+        bool eval(atom * a, bool neg) {
+            return a->is_ineq_atom() ? eval_ineq(to_ineq_atom(a), neg) : eval_root(to_root_atom(a), neg);
+        }
+
+        svector<int> m_add_signs_tmp;
+        void add(poly * p, var x, sign_table & t) {
+            SASSERT(m_pm.max_var(p) <= x);
+            if (m_pm.max_var(p) < x) {
+                t.add_const(eval_sign(p));
+            }
+            else {
+                // isolate roots of p
+                scoped_anum_vector & roots = m_add_roots_tmp;
+                svector<int>       & signs = m_add_signs_tmp;
+                roots.reset();
+                signs.reset();
+                TRACE("nlsat_evaluator", tout << "x: " << x << " max_var(p): " << m_pm.max_var(p) << "\n";);
+                // Note: I added undef_var_assignment in the following statement, to allow us to obtain the infeasible interval sets
+                // even when the maximal variable is assigned. I need this feature to minimize conflict cores.
+                m_am.isolate_roots(polynomial_ref(p, m_pm), undef_var_assignment(m_assignment, x), roots, signs);
+                t.add(roots, signs);
+            }
+        }
+
+        // Evalute the sign of p1^e1*...*pn^en (of atom a) in cell c of table t. 
+        int sign_at(ineq_atom * a, sign_table const & t, unsigned c) const {
+            int sign = 1;
+            unsigned num_ps = a->size();
+            for (unsigned i = 0; i < num_ps; i++) {
+                int curr_sign = t.sign_at(i, c);
+                TRACE("nlsat_evaluator_bug", tout << "sign of i: " << i << " at cell " << c << "\n"; 
+                      m_pm.display(tout, a->p(i)); 
+                      tout << "\nsign: " << curr_sign << "\n";);
+                if (a->is_even(i) && curr_sign < 0)
+                    curr_sign = 1;
+                sign = sign * curr_sign;
+                if (sign == 0)
+                    break;
+            }
+            return sign;
+        }
+        
+        interval_set_ref infeasible_intervals(ineq_atom * a, bool neg) {
+            sign_table & table = m_sign_table_tmp;
+            table.reset();
+            unsigned num_ps = a->size();
+            var x = a->max_var();
+            for (unsigned i = 0; i < num_ps; i++) {
+                add(a->p(i), x, table);
+                TRACE("nlsat_evaluator_bug", tout << "table after:\n"; m_pm.display(tout, a->p(i)); tout << "\n"; table.display_raw(tout);); 
+                
+            }
+            TRACE("nlsat_evaluator", 
+                  tout << "sign table for:\n"; 
+                  for (unsigned i = 0; i < num_ps; i++) { m_pm.display(tout, a->p(i)); tout << "\n"; }
+                  table.display(tout););
+
+            interval_set_ref result(m_ism);
+            interval_set_ref set(m_ism);
+            literal jst(a->bvar(), neg);
+            atom::kind k = a->get_kind();
+            
+            anum dummy;
+            bool prev_sat         = true;
+            bool prev_inf         = true;
+            bool prev_open        = true;
+            unsigned prev_root_id = UINT_MAX;
+            
+            unsigned num_cells = table.num_cells();
+            for (unsigned c = 0; c < num_cells; c++) {
+                TRACE("nlsat_evaluator",
+                      tout << "cell: " << c << "\n";
+                      tout << "prev_sat: " << prev_sat << "\n";
+                      tout << "prev_inf: " << prev_inf << "\n";
+                      tout << "prev_open: " << prev_open << "\n";
+                      tout << "prev_root_id: " << prev_root_id << "\n";
+                      tout << "processing cell: " << c << "\n";
+                      tout << "interval_set so far:\n" << result << "\n";);
+                int sign = sign_at(a, table, c);
+                TRACE("nlsat_evaluator", tout << "sign: " << sign << "\n";);
+                if (satisfied(sign, k, neg)) {
+                    // current cell is satisfied
+                    if (!prev_sat) {
+                        SASSERT(c > 0);
+                        // add interval
+                        bool curr_open;
+                        unsigned curr_root_id;
+                        if (table.is_section(c)) {
+                            curr_open    = true;
+                            curr_root_id = table.get_root_id(c);
+                        }
+                        else {
+                            SASSERT(table.is_section(c-1));
+                            curr_open    = false;
+                            curr_root_id = table.get_root_id(c-1);
+                        }
+                        set = m_ism.mk(prev_open, prev_inf, table.get_root(prev_root_id),
+                                       curr_open, false,    table.get_root(curr_root_id), jst);
+                        result = m_ism.mk_union(result, set);   
+                        prev_sat = true;
+                    }
+                }
+                else {
+                    // current cell is not satisfied
+                    if (prev_sat) {
+                        if (c == 0) {
+                            if (num_cells == 1) {
+                                // (-oo, oo)
+                                result = m_ism.mk(true, true, dummy, true, true, dummy, jst); 
+                            }
+                            else {
+                                // save -oo as begining of infeasible interval
+                                prev_open    = true;
+                                prev_inf     = true;
+                                prev_root_id = UINT_MAX;
+                            }
+                        }
+                        else { 
+                            SASSERT(c > 0);
+                            prev_inf     = false;
+                            if (table.is_section(c)) {
+                                prev_open    = false;
+                                prev_root_id = table.get_root_id(c); 
+                                TRACE("nlsat_evaluator", tout << "updated prev_root_id: " << prev_root_id << " using cell: " << c << "\n";);
+                            }
+                            else {
+                                SASSERT(table.is_section(c-1));
+                                prev_open    = true;
+                                prev_root_id = table.get_root_id(c-1);
+                                TRACE("nlsat_evaluator", tout << "updated prev_root_id: " << prev_root_id << " using cell: " << (c - 1) << "\n";);
+                            }
+                        }
+                        prev_sat = false;
+                    }
+                    if (c == num_cells - 1) {
+                        // last cell add interval with  (prev, oo)
+                        set = m_ism.mk(prev_open, prev_inf, table.get_root(prev_root_id),
+                                       true, true, dummy, jst);
+                        result = m_ism.mk_union(result, set);
+                    } 
+                }
+            }
+            TRACE("nlsat_evaluator", tout << "interval_set: " << result << "\n";);
+            return result;
+        }
+
+        interval_set_ref infeasible_intervals(root_atom * a, bool neg) {
+            atom::kind k = a->get_kind();
+            unsigned i = a->i();
+            SASSERT(i > 0);
+            literal jst(a->bvar(), neg);
+            anum dummy;
+            scoped_anum_vector & roots = m_tmp_values;
+            roots.reset();
+            var x = a->max_var();
+            // Note: I added undef_var_assignment in the following statement, to allow us to obtain the infeasible interval sets
+            // even when the maximal variable is assigned. I need this feature to minimize conflict cores.
+            m_am.isolate_roots(polynomial_ref(a->p(), m_pm), undef_var_assignment(m_assignment, x), roots);
+            interval_set_ref result(m_ism);
+
+            if (i > roots.size()) {
+                // p does have sufficient roots
+                // atom is false by definition
+                if (neg) {
+                    result = m_ism.mk_empty(); 
+                }
+                else {
+                    result = m_ism.mk(true, true, dummy, true, true, dummy, jst); // (-oo, oo)
+                }
+            }
+            else {
+                anum const & r_i = roots[i-1];
+                switch (k) {
+                case atom::ROOT_EQ:
+                    if (neg) {
+                        result = m_ism.mk(false, false, r_i, false, false, r_i, jst); // [r_i, r_i]
+                    }
+                    else {
+                        interval_set_ref s1(m_ism), s2(m_ism);
+                        s1 = m_ism.mk(true, true, dummy, true, false, r_i, jst); // (-oo, r_i)
+                        s2 = m_ism.mk(true, false, r_i, true, true, dummy, jst); // (r_i, oo)
+                        result = m_ism.mk_union(s1, s2);
+                    }
+                    break;
+                case atom::ROOT_LT:
+                    if (neg)
+                        result = m_ism.mk(true, true, dummy, true, false, r_i, jst); // (-oo, r_i)
+                    else
+                        result = m_ism.mk(false, false, r_i, true, true, dummy, jst); // [r_i, oo)
+                    break;
+                case atom::ROOT_GT:
+                    if (neg) 
+                        result = m_ism.mk(true, false, r_i, true, true, dummy, jst); // (r_i, oo)
+                    else
+                        result = m_ism.mk(true, true, dummy, false, false, r_i, jst); // (-oo, r_i]
+                    break;
+                case atom::ROOT_LE:
+                    if (neg)
+                        result = m_ism.mk(true, true, dummy, false, false, r_i, jst); // (-oo, r_i]
+                    else
+                        result = m_ism.mk(true, false, r_i, true, true, dummy, jst); // (r_i, oo)
+                    break;
+                case atom::ROOT_GE:
+                    if (neg) 
+                        result = m_ism.mk(false, false, r_i, true, true, dummy, jst); // [r_i, oo)
+                    else
+                        result = m_ism.mk(true, true, dummy, true, false, r_i, jst); // (-oo, r_i)
+                    break;
+                default:
+                    UNREACHABLE();
+                    break;
+                }
+            }
+            TRACE("nlsat_evaluator", tout << "interval_set: " << result << "\n";);
+            return result;
+        }
+        
+        interval_set_ref infeasible_intervals(atom * a, bool neg) {
+            return a->is_ineq_atom() ? infeasible_intervals(to_ineq_atom(a), neg) : infeasible_intervals(to_root_atom(a), neg); 
+        }
+    };
+    
+    evaluator::evaluator(assignment const & x2v, pmanager & pm, small_object_allocator & allocator) {
+        m_imp = alloc(imp, x2v, pm, allocator);
+    }
+
+    evaluator::~evaluator() {
+        dealloc(m_imp);
+    }
+
+    interval_set_manager & evaluator::ism() const {
+        return m_imp->m_ism;
+    }
+
+    bool evaluator::eval(atom * a, bool neg) {
+        return m_imp->eval(a, neg);
+    }
+        
+    interval_set_ref evaluator::infeasible_intervals(atom * a, bool neg) {
+        return m_imp->infeasible_intervals(a, neg);
+    }
+
+    void evaluator::push() {
+        // do nothing
+    }
+
+    void evaluator::pop(unsigned num_scopes) {
+        // do nothing
+    }
+};

src/nlsat/nlsat_evaluator.h

@@ -0,0 +1,61 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_evaluator.h
+
+Abstract:
+
+    Helper class for computing the infeasible intervals of an
+    arithmetic literal.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-12.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_EVALUATOR_H_
+#define _NLSAT_EVALUATOR_H_
+
+#include"nlsat_types.h"
+#include"nlsat_assignment.h"
+#include"nlsat_interval_set.h"
+
+namespace nlsat {
+
+    class evaluator {
+        struct imp;
+        imp *  m_imp;
+    public:
+        evaluator(assignment const & x2v, pmanager & pm, small_object_allocator & allocator);
+        ~evaluator();
+
+        interval_set_manager & ism() const;
+
+        /**
+           \brief Evaluate the given literal in the current model.
+           All variables in the atom must be assigned.
+           
+           The result is true if the literal is satisfied, and false otherwise.
+        */
+        bool eval(atom * a, bool neg);
+        
+        /**
+           \brief Return the infeasible interval set for the given literal.
+           All but the a->max_var() must be assigned in the current model.
+
+           Let x be a->max_var(). Then, the resultant set specifies which
+           values of x falsify the given literal.
+        */
+        interval_set_ref infeasible_intervals(atom * a, bool neg);
+
+        void push();
+        void pop(unsigned num_scopes);
+    };
+    
+};
+
+#endif

src/nlsat/nlsat_explain.h

@@ -0,0 +1,64 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_explain.h
+
+Abstract:
+
+    Functor that implements the "explain" procedure defined in Dejan and Leo's paper.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-13.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_EXPLAIN_H_
+#define _NLSAT_EXPLAIN_H_
+
+#include"nlsat_solver.h"
+#include"nlsat_scoped_literal_vector.h"
+#include"polynomial_cache.h"
+
+namespace nlsat {
+    class evaluator;
+    
+    class explain {
+        struct imp;
+        imp * m_imp;
+    public:
+        explain(solver & s, assignment const & x2v, polynomial::cache & u, atom_vector const & atoms, atom_vector const & x2eq,
+                evaluator & ev);
+        ~explain();
+
+        void reset();
+        void set_simplify_cores(bool f);
+        void set_full_dimensional(bool f);
+        void set_minimize_cores(bool f);
+
+        /**
+           \brief Given a set of literals ls[0], ... ls[n-1] s.t.
+           - n > 0
+           - all of them are arithmetic literals.
+           - all of them have the same maximal variable.
+           - (ls[0] and ... and ls[n-1]) is infeasible in the current interpretation.
+
+           Let x be the maximal variable in {ls[0], ..., ls[n-1]}.
+
+           Remark: the current interpretation assigns all variables in ls[0], ..., ls[n-1] but x.
+
+           This procedure stores in result a set of literals: s_1, ..., s_m
+           s.t.
+                 - (s_1 or ... or s_m or ~ls[0] or ... or ~ls[n-1]) is a valid clause
+                 - s_1, ..., s_m do not contain variable x.
+                 - s_1, ..., s_m are false in the current interpretation
+        */
+        void operator()(unsigned n, literal const * ls, scoped_literal_vector & result);
+    };
+
+};
+
+#endif

src/nlsat/nlsat_interval_set.cpp

@@ -0,0 +1,762 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_interval_set.cpp
+
+Abstract:
+
+    Sets of disjoint infeasible intervals.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-11.
+
+Revision History:
+
+--*/
+#include"nlsat_interval_set.h"
+#include"algebraic_numbers.h"
+#include"buffer.h"
+
+namespace nlsat {
+
+    struct interval {
+        unsigned  m_lower_open:1;
+        unsigned  m_upper_open:1;
+        unsigned  m_lower_inf:1;
+        unsigned  m_upper_inf:1;
+        literal   m_justification;
+        anum      m_lower;
+        anum      m_upper;
+    };
+
+    class interval_set {
+    public:
+        static unsigned get_obj_size(unsigned num) { return sizeof(interval_set) + num*sizeof(interval); }
+        unsigned  m_num_intervals;
+        unsigned  m_ref_count:31;
+        unsigned  m_full:1;
+        interval  m_intervals[0];
+    };
+
+    void display(std::ostream & out, anum_manager & am, interval const & curr) {
+        if (curr.m_lower_inf) {
+            out << "(-oo, ";
+        }
+        else {
+            if (curr.m_lower_open)
+                    out << "(";
+            else
+                out << "[";
+            am.display_decimal(out, curr.m_lower);
+            out << ", ";
+        }
+        if (curr.m_justification.sign())
+            out << "~";
+        out << "p";
+        out << curr.m_justification.var() << ", ";
+        if (curr.m_upper_inf) {
+            out << "oo)";
+        }
+        else {
+            am.display_decimal(out, curr.m_upper);
+            if (curr.m_upper_open)
+                out << ")";
+            else
+                out << "]";
+        }
+    }
+
+    bool check_interval(anum_manager & am, interval const & i) {
+        if (i.m_lower_inf) {
+            SASSERT(i.m_lower_open);
+        }
+        if (i.m_upper_inf) {
+            SASSERT(i.m_upper_open);
+        }
+        if (!i.m_lower_inf && !i.m_upper_inf) {
+            int s = am.compare(i.m_lower, i.m_upper);
+            TRACE("nlsat_interval", tout << "lower: "; am.display_decimal(tout, i.m_lower); tout << ", upper: "; am.display_decimal(tout, i.m_upper);
+                  tout << "\ns: " << s << "\n";);
+            SASSERT(s <= 0);
+            if (s == 0) {
+                SASSERT(!i.m_lower_open && !i.m_upper_open);
+            }
+        }
+        return true;
+    }
+
+    bool check_no_overlap(anum_manager & am, interval const & curr, interval const & next) {
+        SASSERT(!curr.m_upper_inf);
+        SASSERT(!next.m_lower_inf);
+        int sign = am.compare(curr.m_upper, next.m_lower);
+        CTRACE("nlsat", sign > 0, display(tout, am, curr); tout << "  "; display(tout, am, next); tout << "\n";);
+        SASSERT(sign <= 0);
+        if (sign == 0) {
+            SASSERT(curr.m_upper_open || next.m_lower_open);
+        }
+        return true;
+    }
+    
+    // Check if the intervals are valid, ordered, and are disjoint.
+    bool check_interval_set(anum_manager & am, unsigned sz, interval const * ints) {
+        for (unsigned i = 0; i < sz; i++) {
+            interval const & curr = ints[i];
+            SASSERT(check_interval(am, curr));
+            if (i < sz - 1) {
+                interval const & next = ints[i+1];
+                SASSERT(check_no_overlap(am, curr, next));
+            }
+        }
+        return true;
+    }
+
+    interval_set_manager::interval_set_manager(anum_manager & m, small_object_allocator & a):
+        m_am(m),
+        m_allocator(a) {
+    }
+     
+    interval_set_manager::~interval_set_manager() {
+    }
+    
+    void interval_set_manager::del(interval_set * s) {
+        if (s == 0)
+            return;
+        unsigned num    = s->m_num_intervals;
+        unsigned obj_sz = interval_set::get_obj_size(num);
+        for (unsigned i = 0; i < num; i++) {
+            m_am.del(s->m_intervals[i].m_lower);
+            m_am.del(s->m_intervals[i].m_upper);
+        }
+        s->~interval_set();
+        m_allocator.deallocate(obj_sz, s);
+    }
+
+    void interval_set_manager::dec_ref(interval_set * s) {
+        SASSERT(s->m_ref_count > 0);
+        s->m_ref_count--;
+        if (s->m_ref_count == 0)
+            del(s);
+    }
+
+    void interval_set_manager::inc_ref(interval_set * s) {
+        s->m_ref_count++;
+    }
+    
+    interval_set * interval_set_manager::mk(bool lower_open, bool lower_inf, anum const & lower, 
+                                            bool upper_open, bool upper_inf, anum const & upper,
+                                            literal justification) {
+        void * mem = m_allocator.allocate(interval_set::get_obj_size(1));
+        interval_set * new_set = new (mem) interval_set();
+        new_set->m_num_intervals = 1;
+        new_set->m_ref_count  = 0;
+        new_set->m_full       = lower_inf && upper_inf;
+        interval * i = new (new_set->m_intervals) interval();
+        i->m_lower_open = lower_open;
+        i->m_lower_inf  = lower_inf;
+        i->m_upper_open = upper_open;
+        i->m_upper_inf  = upper_inf;
+        i->m_justification = justification;
+        if (!lower_inf)
+            m_am.set(i->m_lower, lower);
+        if (!upper_inf)
+            m_am.set(i->m_upper, upper);
+        SASSERT(check_interval_set(m_am, 1, new_set->m_intervals));
+        return new_set;
+    }
+
+    inline int compare_lower_lower(anum_manager & am, interval const & i1, interval const & i2) {
+        if (i1.m_lower_inf && i2.m_lower_inf)
+            return 0;
+        if (i1.m_lower_inf)
+            return -1;
+        if (i2.m_lower_inf)
+            return 1;
+        SASSERT(!i1.m_lower_inf && !i2.m_lower_inf);
+        int s = am.compare(i1.m_lower, i2.m_lower);
+        if (s != 0)
+            return s;
+        if (i1.m_lower_open == i2.m_lower_open)
+            return 0;
+        if (i1.m_lower_open)
+            return 1;
+        else
+            return -1;
+    }
+
+    inline int compare_upper_upper(anum_manager & am, interval const & i1, interval const & i2) {
+        if (i1.m_upper_inf && i2.m_upper_inf)
+            return 0;
+        if (i1.m_upper_inf)
+            return 1;
+        if (i2.m_upper_inf)
+            return -1;
+        SASSERT(!i1.m_upper_inf && !i2.m_upper_inf);
+        int s = am.compare(i1.m_upper, i2.m_upper);
+        if (s != 0)
+            return s;
+        if (i1.m_upper_open == i2.m_upper_open)
+            return 0;
+        if (i1.m_upper_open)
+            return -1;
+        else
+            return 1;
+    }
+
+    inline int compare_upper_lower(anum_manager & am, interval const & i1, interval const & i2) {
+        if (i1.m_upper_inf || i2.m_lower_inf)
+            return 1;
+        SASSERT(!i1.m_upper_inf && !i2.m_lower_inf);
+        int s = am.compare(i1.m_upper, i2.m_lower);
+        if (s != 0)
+            return s;
+        if (!i1.m_upper_open && !i2.m_lower_open)
+            return 0;
+        return -1;
+    }
+    
+    typedef sbuffer<interval, 128> interval_buffer;
+
+    // Given two interval in an interval set s.t. curr occurs before next.
+    // We say curr and next are "adjacent" iff
+    //      there is no "space" between them.
+    bool adjacent(anum_manager & am, interval const & curr, interval const & next) {
+        SASSERT(!curr.m_upper_inf);
+        SASSERT(!next.m_lower_inf);
+        int sign = am.compare(curr.m_upper, next.m_lower);
+        SASSERT(sign <= 0);
+        if (sign == 0) {
+            SASSERT(curr.m_upper_open || next.m_lower_open);
+            return !curr.m_upper_open || !next.m_lower_open;
+        }
+        return false;
+    }
+
+    inline void push_back(anum_manager & am, interval_buffer & buf, 
+                          bool lower_open, bool lower_inf, anum const & lower, 
+                          bool upper_open, bool upper_inf, anum const & upper,
+                          literal justification) {
+        buf.push_back(interval());
+        interval & i = buf.back();
+        i.m_lower_open = lower_open;
+        i.m_lower_inf  = lower_inf;
+        am.set(i.m_lower, lower);
+        i.m_upper_open = upper_open;
+        i.m_upper_inf  = upper_inf;
+        am.set(i.m_upper, upper);
+        i.m_justification = justification;
+        SASSERT(check_interval(am, i));
+    }
+
+    inline void push_back(anum_manager & am, interval_buffer & buf, interval const & i) {
+        push_back(am, buf,
+                  i.m_lower_open, i.m_lower_inf, i.m_lower,
+                  i.m_upper_open, i.m_upper_inf, i.m_upper,
+                  i.m_justification);
+    }
+
+    inline interval_set * mk_interval(small_object_allocator & allocator, interval_buffer & buf, bool full) {
+        unsigned sz = buf.size();
+        void * mem = allocator.allocate(interval_set::get_obj_size(sz));
+        interval_set * new_set = new (mem) interval_set();
+        new_set->m_full = full;
+        new_set->m_ref_count  = 0;
+        new_set->m_num_intervals = sz;
+        memcpy(new_set->m_intervals, buf.c_ptr(), sizeof(interval)*sz);
+        return new_set;
+    }
+
+    interval_set * interval_set_manager::mk_union(interval_set const * s1, interval_set const * s2) {
+        TRACE("nlsat_interval", tout << "mk_union\ns1: "; display(tout, s1); tout << "\ns2: "; display(tout, s2); tout << "\n";);
+        if (s1 == 0 || s1 == s2)
+            return const_cast<interval_set*>(s2);
+        if (s2 == 0)
+            return const_cast<interval_set*>(s1);
+        if (s1->m_full)
+            return const_cast<interval_set*>(s1);
+        if (s2->m_full)
+            return const_cast<interval_set*>(s2);
+        interval_buffer result;
+        unsigned sz1 = s1->m_num_intervals;
+        unsigned sz2 = s2->m_num_intervals;
+        unsigned i1  = 0;
+        unsigned i2  = 0;
+        while (true) {
+            if (i1 >= sz1) {
+                while (i2 < sz2) {
+                    TRACE("nlsat_interval", tout << "adding remaining intervals from s2: "; nlsat::display(tout, m_am, s2->m_intervals[i2]); tout << "\n";);
+                    push_back(m_am, result, s2->m_intervals[i2]);
+                    i2++;
+                }
+                break;
+            }
+            if (i2 >= sz2) {
+                while (i1 < sz1) {
+                    TRACE("nlsat_interval", tout << "adding remaining intervals from s1: "; nlsat::display(tout, m_am, s1->m_intervals[i1]); tout << "\n";);
+                    push_back(m_am, result, s1->m_intervals[i1]);
+                    i1++;
+                }
+                break;
+            }
+            interval const & int1 = s1->m_intervals[i1];
+            interval const & int2 = s2->m_intervals[i2];
+            int l1_l2_sign = compare_lower_lower(m_am, int1, int2);
+            int u1_u2_sign = compare_upper_upper(m_am, int1, int2);
+            TRACE("nlsat_interval", 
+                  tout << "i1: " << i1 << ", i2: " << i2 << "\n";
+                  tout << "int1: "; nlsat::display(tout, m_am, int1); tout << "\n";
+                  tout << "int2: "; nlsat::display(tout, m_am, int2); tout << "\n";);
+            if (l1_l2_sign <= 0) {
+                if (u1_u2_sign == 0) {
+                    // Cases:
+                    // 1)  [     ]
+                    //     [     ]
+                    //
+                    // 2) [     ]
+                    //      [   ]
+                    //
+                    TRACE("nlsat_interval", tout << "l1_l2_sign <= 0, u1_u2_sign == 0\n";);
+                    push_back(m_am, result, int1);
+                    i1++;
+                    i2++;
+                }
+                else if (u1_u2_sign > 0) {
+                    // Cases:
+                    //
+                    // 1) [        ]
+                    //    [     ]
+                    //
+                    // 2) [        ]
+                    //      [   ]
+                    i2++;
+                    TRACE("nlsat_interval", tout << "l1_l2_sign <= 0, u1_u2_sign > 0\n";);
+                    // i1 may consume other intervals of s2
+                }
+                else {
+                    SASSERT(u1_u2_sign < 0);
+                    int u1_l2_sign = compare_upper_lower(m_am, int1, int2);
+                    if (u1_l2_sign < 0) {
+                        SASSERT(l1_l2_sign < 0);
+                        // Cases:
+                        // 1)   [      ]
+                        //                [     ]
+                        TRACE("nlsat_interval", tout << "l1_l2_sign <= 0, u1_u2_sign < 0, u1_l2_sign < 0\n";);
+                        push_back(m_am, result, int1);
+                        i1++;
+                    }
+                    else if (u1_l2_sign == 0) {
+                        SASSERT(l1_l2_sign <= 0);
+                        SASSERT(!int1.m_upper_open && !int2.m_lower_open);
+                        SASSERT(!int2.m_lower_inf);
+                        TRACE("nlsat_interval", tout << "l1_l2_sign <= 0, u1_u2_sign < 0, u1_l2_sign == 0\n";);
+                        // Cases:
+                        if (l1_l2_sign != 0) {
+                            SASSERT(l1_l2_sign < 0);
+                            // 1)   [   ]
+                            //          [    ]
+                            SASSERT(!int2.m_lower_open);
+                            push_back(m_am, result, 
+                                      int1.m_lower_open, int1.m_lower_inf,  int1.m_lower,
+                                      true /* open */, false /* not +oo */, int1.m_upper, 
+                                      int1.m_justification); 
+                            i1++;
+                        }
+                        else {
+                            SASSERT(l1_l2_sign == 0);
+                            // 2)   u          <<< int1 is a singleton
+                            //      [     ]
+                            // just consume int1
+                            i1++;
+                        }
+                    }
+                    else {
+                        SASSERT(l1_l2_sign <= 0);
+                        SASSERT(u1_u2_sign < 0);
+                        SASSERT(u1_l2_sign > 0);
+                        TRACE("nlsat_interval", tout << "l1_l2_sign <= 0, u1_u2_sign < 0, u1_l2_sign > 0\n";);
+                        if (l1_l2_sign == 0) {
+                            // Case:
+                            // 1)   [      ]
+                            //      [          ]
+                            // just consume int1
+                            i1++;
+                        }
+                        else {
+                            SASSERT(l1_l2_sign < 0);
+                            SASSERT(u1_u2_sign < 0);
+                            SASSERT(u1_l2_sign > 0);
+                            // 2) [        ]
+                            //         [        ]
+                            push_back(m_am, result, 
+                                      int1.m_lower_open,     int1.m_lower_inf, int1.m_lower,
+                                      !int2.m_lower_open, false /* not +oo */, int2.m_lower,
+                                      int1.m_justification); 
+                            i1++;
+                        }
+                    }
+                }
+            }
+            else {
+                SASSERT(l1_l2_sign > 0);
+                if (u1_u2_sign == 0) {
+                    TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign == 0\n";);
+                    // Case:
+                    // 1)    [  ]
+                    //    [     ]
+                    //
+                    push_back(m_am, result, int2);
+                    i1++;
+                    i2++;
+                }
+                else if (u1_u2_sign < 0) {
+                    TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0\n";);
+                    // Case:
+                    // 1)   [   ]
+                    //    [       ]
+                    i1++;
+                    // i2 may consume other intervals of s1
+                }
+                else {
+                    int u2_l1_sign = compare_upper_lower(m_am, int2, int1);
+                    if (u2_l1_sign < 0) {
+                        TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0, u2_l1_sign < 0\n";);
+                        // Case:
+                        // 1)           [      ]
+                        //     [     ]
+                        push_back(m_am, result, int2);
+                        i2++;
+                    }
+                    else if (u2_l1_sign == 0) {
+                        TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0, u2_l1_sign == 0\n";);
+                        SASSERT(!int1.m_lower_open && !int2.m_upper_open);
+                        SASSERT(!int1.m_lower_inf);
+                        // Case:
+                        //      [    ]   
+                        //  [   ]
+                        push_back(m_am, result, 
+                                  int2.m_lower_open,     int2.m_lower_inf, int2.m_lower,
+                                  true /* open */,    false /* not +oo */, int2.m_upper, 
+                                  int2.m_justification); 
+                        i2++;
+                    }
+                    else {
+                        TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0, u2_l1_sign > 0\n";);
+                        SASSERT(l1_l2_sign > 0);
+                        SASSERT(u1_u2_sign > 0);
+                        SASSERT(u2_l1_sign > 0);
+                        // Case:
+                        //     [        ]
+                        // [        ]
+                        push_back(m_am, result, 
+                                  int2.m_lower_open,     int2.m_lower_inf, int2.m_lower,
+                                  !int1.m_lower_open, false /* not +oo */, int1.m_lower,
+                                  int2.m_justification); 
+                        i2++;
+                    }
+                }
+            }
+            SASSERT(result.size() <= 1 ||
+                    check_no_overlap(m_am, result[result.size() - 2], result[result.size() - 1]));
+        }
+
+        SASSERT(!result.empty());
+        SASSERT(check_interval_set(m_am, result.size(), result.c_ptr()));
+        // Compress
+        // Remark: we only combine adjacent intervals when they have the same justification
+        unsigned j  = 0;
+        unsigned sz = result.size(); 
+        for (unsigned i = 1; i < sz; i++) {
+            interval & curr = result[j];
+            interval & next = result[i];
+            if (curr.m_justification == next.m_justification && 
+                adjacent(m_am, curr, next)) {
+                // merge them
+                curr.m_upper_inf  = next.m_upper_inf;
+                curr.m_upper_open = next.m_upper_open;
+                m_am.swap(curr.m_upper, next.m_upper);
+            }
+            else {
+                j++;
+                if (i != j) {
+                    interval & next_curr = result[j];
+                    next_curr.m_lower_inf = next.m_lower_inf;
+                    next_curr.m_lower_open = next.m_lower_open;
+                    m_am.swap(next_curr.m_lower, next.m_lower);
+                    next_curr.m_upper_inf = next.m_upper_inf;
+                    next_curr.m_upper_open = next.m_upper_open;
+                    m_am.swap(next_curr.m_upper, next.m_upper);
+                    next_curr.m_justification = next.m_justification;
+                }
+            }
+        }
+        j++;
+        for (unsigned i = j; i < sz; i++) {
+            interval & curr = result[i];
+            m_am.del(curr.m_lower);
+            m_am.del(curr.m_upper);
+        }
+        result.shrink(j);
+        SASSERT(check_interval_set(m_am, result.size(), result.c_ptr()));
+        sz = j;
+        SASSERT(sz >= 1);
+        bool found_slack  = !result[0].m_lower_inf || !result[sz-1].m_upper_inf;
+        // Check if full
+        for (unsigned i = 0; i < sz - 1 && !found_slack; i++) {
+            if (!adjacent(m_am, result[i], result[i+1]))
+                found_slack = true;
+        }
+        // Create new interval set
+        interval_set * new_set = mk_interval(m_allocator, result, !found_slack);
+        SASSERT(check_interval_set(m_am, sz, new_set->m_intervals));
+        return new_set;
+    }
+
+    bool interval_set_manager::is_full(interval_set const * s) {
+        if (s == 0)
+            return false;
+        return s->m_full == 1;
+    }
+
+    unsigned interval_set_manager::num_intervals(interval_set const * s) const {
+        if (s == 0) return 0;
+        return s->m_num_intervals;
+    }
+    
+    bool interval_set_manager::subset(interval_set const * s1, interval_set const * s2) {
+        if (s1 == s2)
+            return true;
+        if (s1 == 0)
+            return true;
+        if (s2 == 0)
+            return false;
+        if (s2->m_full)
+            return true;
+        if (s1->m_full)
+            return false;
+        unsigned sz1 = s1->m_num_intervals;
+        unsigned sz2 = s2->m_num_intervals;
+        SASSERT(sz1 > 0 && sz2 > 0);
+        unsigned i1  = 0;
+        unsigned i2  = 0;
+        while (i1 < sz1 && i2 < sz2) {
+            interval const & int1 = s1->m_intervals[i1];
+            interval const & int2 = s2->m_intervals[i2];
+            TRACE("nlsat_interval", tout << "subset main loop, i1: " << i1 << ", i2: " << i2 << "\n";
+                  tout << "int1: "; nlsat::display(tout, m_am, int1); tout << "\n";
+                  tout << "int2: "; nlsat::display(tout, m_am, int2); tout << "\n";);
+            if (compare_lower_lower(m_am, int1, int2) < 0) {
+                TRACE("nlsat_interval", tout << "done\n";);
+                // interval [int1.lower1, int2.lower2] is not in s2
+                // s1: [ ...
+                // s2:    [ ...
+                return false;
+            }
+            while (i2 < sz2) {
+                interval const & int2 = s2->m_intervals[i2];
+                TRACE("nlsat_interval", tout << "inner loop, i2: " << i2 << "\n";
+                      tout << "int2: "; nlsat::display(tout, m_am, int2); tout << "\n";);
+                int u1_u2_sign = compare_upper_upper(m_am, int1, int2);
+                if (u1_u2_sign == 0) {
+                    TRACE("nlsat_interval", tout << "case 1, break\n";);
+                    // consume both
+                    // s1: ... ]
+                    // s2: ... ]
+                    i1++;
+                    i2++;
+                    break;
+                }
+                else if (u1_u2_sign < 0) {
+                    TRACE("nlsat_interval", tout << "case 2, break\n";);
+                    // consume only int1, int2 may cover other intervals of s1 
+                    // s1:  ... ]
+                    // s2:    ... ]
+                    i1++;
+                    break;
+                }
+                else {
+                    SASSERT(u1_u2_sign > 0);
+                    int u2_l1_sign = compare_upper_lower(m_am, int2, int1);
+                    TRACE("nlsat_interval", tout << "subset, u2_l1_sign: " << u2_l1_sign << "\n";);
+                    if (u2_l1_sign < 0) {
+                        TRACE("nlsat_interval", tout << "case 3, break\n";);
+                        // s1:           [ ...
+                        // s2: [ ... ]  ...
+                        i2++;
+                        break;
+                    }
+                    SASSERT(u2_l1_sign >= 0);
+                    // s1:      [ ...  ]
+                    // s2: [ ...    ]
+                    if (i2 == sz2 - 1) {
+                        TRACE("nlsat_interval", tout << "case 4, done\n";);
+                        // s1:   ... ]
+                        // s2: ...]
+                        // the interval [int2.upper, int1.upper] is not in s2
+                        return false; // last interval of s2
+                    }
+                    interval const & next2 = s2->m_intervals[i2+1];
+                    if (!adjacent(m_am, int2, next2)) {
+                        TRACE("nlsat_interval", tout << "not adjacent, done\n";);
+                        // s1:   ... ]
+                        // s2: ... ]   [
+                        // the interval [int2.upper, min(int1.upper, next2.lower)] is not in s2
+                        return false;
+                    }
+                    TRACE("nlsat_interval", tout << "continue..\n";);
+                    // continue with adjacent interval of s2
+                    // s1:    ...  ]
+                    // s2:  ..][ ...
+                    i2++;
+                }
+            }
+        }
+        return i1 == sz1;
+    }
+
+    bool interval_set_manager::set_eq(interval_set const * s1, interval_set const * s2) {
+        if (s1 == 0 || s2 == 0)
+            return s1 == s2;
+        if (s1->m_full || s2->m_full)
+            return s1->m_full == s2->m_full;
+        // TODO: check if bottleneck, then replace simple implementation
+        return subset(s1, s2) && subset(s2, s1);
+    }
+        
+    bool interval_set_manager::eq(interval_set const * s1, interval_set const * s2) {
+        if (s1 == 0 || s2 == 0)
+            return s1 == s2;
+        if (s1->m_num_intervals != s2->m_num_intervals)
+            return false;
+        for (unsigned i = 0; i < s1->m_num_intervals; i++) {
+            interval const & int1 = s1->m_intervals[i];
+            interval const & int2 = s2->m_intervals[i]; 
+            if (int1.m_lower_inf  != int2.m_lower_inf ||
+                int1.m_lower_open != int2.m_lower_open ||
+                int1.m_upper_inf  != int2.m_upper_inf ||
+                int1.m_upper_open != int2.m_upper_open ||
+                int1.m_justification != int2.m_justification ||
+                !m_am.eq(int1.m_lower, int2.m_lower) ||
+                !m_am.eq(int1.m_upper, int2.m_upper))
+                return false;
+        }
+        return true;
+    }
+
+    void interval_set_manager::get_justifications(interval_set const * s, literal_vector & js) {
+        js.reset();
+        unsigned num = num_intervals(s);
+        for (unsigned i = 0; i < num; i++) {
+            literal l     = s->m_intervals[i].m_justification;
+            unsigned lidx = l.index();
+            if (m_already_visited.get(lidx, false))
+                continue;
+            m_already_visited.setx(lidx, true, false);
+            js.push_back(l);
+        }
+        for (unsigned i = 0; i < num; i++) {
+            literal l     = s->m_intervals[i].m_justification;
+            unsigned lidx = l.index();
+            m_already_visited[lidx] = false;
+        }
+    }
+    
+    interval_set * interval_set_manager::get_interval(interval_set const * s, unsigned idx) const {
+        SASSERT(idx < num_intervals(s));
+        interval_buffer result;
+        push_back(m_am, result, s->m_intervals[idx]);
+        bool found_slack  = !result[0].m_lower_inf || !result[0].m_upper_inf;
+        interval_set * new_set = mk_interval(m_allocator, result, !found_slack);
+        SASSERT(check_interval_set(m_am, result.size(), new_set->m_intervals));
+        return new_set;
+    }
+
+    void interval_set_manager::peek_in_complement(interval_set const * s, anum & w, bool randomize) {
+        SASSERT(!is_full(s));
+        if (s == 0) {
+            if (randomize) {
+                int num = m_rand() % 2 == 0 ? 1 : -1;
+#define MAX_RANDOM_DEN_K 4
+                int den_k = (m_rand() % MAX_RANDOM_DEN_K);
+                int den   = 1 << den_k;
+                scoped_mpq _w(m_am.qm());
+                m_am.qm().set(_w, num, den);
+                m_am.set(w, _w);
+                return;
+            }
+            else {
+                m_am.set(w, 0);
+                return;
+            }
+        }
+        
+        unsigned n = 0;
+        
+        unsigned num = num_intervals(s);
+        if (!s->m_intervals[0].m_lower_inf) {
+            // lower is not -oo
+            n++;
+            m_am.int_lt(s->m_intervals[0].m_lower, w);
+            if (!randomize)
+                return;
+        }
+        if (!s->m_intervals[num-1].m_upper_inf) {
+            // upper is not oo
+            n++;
+            if (n == 1 || m_rand()%n == 0)
+                m_am.int_gt(s->m_intervals[num-1].m_upper, w);
+            if (!randomize)
+                return;
+        }
+        
+        // Try to find a gap that is not an unit.
+        for (unsigned i = 1; i < num; i++) {
+            if (m_am.lt(s->m_intervals[i-1].m_upper, s->m_intervals[i].m_lower)) {
+                n++;
+                if (n == 1 || m_rand()%n == 0)
+                    m_am.select(s->m_intervals[i-1].m_upper, s->m_intervals[i].m_lower, w);
+                if (!randomize)
+                    return;
+            }
+        }
+        
+        if (n > 0)
+            return;
+        
+        // Try to find a rational
+        unsigned irrational_i = UINT_MAX;
+        for (unsigned i = 1; i < num; i++) {
+            if (s->m_intervals[i-1].m_upper_open && s->m_intervals[i].m_lower_open) {
+                SASSERT(m_am.eq(s->m_intervals[i-1].m_upper, s->m_intervals[i].m_lower)); // otherwise we would have found it in the previous step
+                if (m_am.is_rational(s->m_intervals[i-1].m_upper)) {
+                    m_am.set(w, s->m_intervals[i-1].m_upper);
+                    return;
+                }
+                if (irrational_i == UINT_MAX)
+                    irrational_i = i-1;
+            }
+        }
+        SASSERT(irrational_i != UINT_MAX);
+        // Last option: peek irrational witness :-(
+        SASSERT(s->m_intervals[irrational_i].m_upper_open && s->m_intervals[irrational_i+1].m_lower_open);
+        m_am.set(w, s->m_intervals[irrational_i].m_upper);
+    }
+
+    void interval_set_manager::display(std::ostream & out, interval_set const * s) const {
+        if (s == 0) {
+            out << "{}";
+            return;
+        }
+        out << "{";
+        for (unsigned i = 0; i < s->m_num_intervals; i++) {
+            if (i > 0)
+                out << ", ";
+            nlsat::display(out, m_am, s->m_intervals[i]);
+        }
+        out << "}";
+        if (s->m_full)
+            out << "*";
+    }
+
+};

src/nlsat/nlsat_interval_set.h

@@ -0,0 +1,123 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_interval_set.h
+
+Abstract:
+
+    Sets of disjoint infeasible intervals.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-11.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_INTERVAL_SET_H_
+#define _NLSAT_INTERVAL_SET_H_
+
+#include"nlsat_types.h"
+
+namespace nlsat {
+
+    class interval_set;
+
+    class interval_set_manager {
+        anum_manager &           m_am;
+        small_object_allocator & m_allocator;
+        svector<char>            m_already_visited;
+        random_gen               m_rand;
+        void del(interval_set * s);
+    public:
+        interval_set_manager(anum_manager & m, small_object_allocator & a);
+        ~interval_set_manager();
+        
+        void set_seed(unsigned s) { m_rand.set_seed(s); }
+
+        /**
+           \brief Return the empty set.
+        */
+        interval_set * mk_empty() { return 0; }
+        
+        /**
+           \brief Return a set of composed of a single interval.
+        */
+        interval_set * mk(bool lower_open, bool lower_inf, anum const & lower, 
+                          bool upper_open, bool upper_inf, anum const & upper,
+                          literal justification);
+        
+        /**
+           \brief Return the union of two sets.
+        */
+        interval_set * mk_union(interval_set const * s1, interval_set const * s2);
+        
+        /**
+           \brief Reference counting
+        */
+        void dec_ref(interval_set * s);
+        void inc_ref(interval_set * s);
+        
+        /**
+           \brief Return true if s is the empty set.
+        */
+        bool is_empty(interval_set const * s) {
+            return s == 0;
+        }
+        
+        /**
+           \brief Return true if the set contains all real numbers.
+        */
+        bool is_full(interval_set const * s);
+
+        /**
+           `\brief Return true if s1 is a subset of s2.
+        */
+        bool subset(interval_set const * s1, interval_set const * s2);
+
+        /**
+           \brief Return true if s1 and s2 cover the same subset of R.
+           The justifications are ignored
+        */
+        bool set_eq(interval_set const * s1, interval_set const * s2);
+        
+        /**
+           \brief Return true if s1 and s2 are the same (the justifications are taking into account).
+        */
+        bool eq(interval_set const * s1, interval_set const * s2);
+
+        /**
+           \brief Return a set of literals that justify s.
+        */
+        void get_justifications(interval_set const * s, literal_vector & js);
+        
+        void display(std::ostream & out, interval_set const * s) const;
+        
+        unsigned num_intervals(interval_set const * s) const;
+        
+        /**
+           \brief (For debugging purposes) Return one of the intervals in s.
+           \pre idx < num_intervals()
+        */
+        interval_set * get_interval(interval_set const * s, unsigned idx) const;
+        
+        /**
+           \brief Select a witness w in the complement of s.
+           
+           \pre !is_full(s)
+        */
+        void peek_in_complement(interval_set const * s, anum & w, bool randomize);
+    };
+
+    typedef obj_ref<interval_set, interval_set_manager> interval_set_ref;
+
+    inline std::ostream & operator<<(std::ostream & out, interval_set_ref const & s) {
+        s.m().display(out, s);
+        return out;
+    }
+
+};
+
+#endif

src/nlsat/nlsat_justification.h

@@ -0,0 +1,90 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_justification.h
+
+Abstract:
+
+    An explanation for a (Boolean) assignment in the 
+    nlsat procedure
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-10.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_JUSTIFICATION_H_
+#define _NLSAT_JUSTIFICATION_H_
+
+#include"nlsat_types.h"
+#include"tptr.h"
+
+namespace nlsat {
+
+    // There are two kinds of justifications in nlsat:
+    //
+    // - clause
+    //
+    // - lazy_justification: it is a set of arithmetic literals s.t.
+    // the maximal variable in each literal is the same. 
+    // The set is inconsistent in the current model.
+    // Thus, our nonlinear procedure may be applied to it
+    // to produce a clause.
+    //
+
+    class lazy_justification {
+        unsigned m_num_literals;
+        literal  m_literals[0];
+    public:
+        static unsigned get_obj_size(unsigned num) { return sizeof(lazy_justification) + sizeof(literal)*num; }
+        lazy_justification(unsigned num, literal const * lits):
+            m_num_literals(num) {
+            memcpy(m_literals, lits, sizeof(literal)*num);
+        }
+        unsigned size() const { return m_num_literals; }
+        literal operator[](unsigned i) const { SASSERT(i < size()); return m_literals[i]; }
+        literal const * lits() const { return m_literals; }
+    };
+
+    class justification {
+        void * m_data;
+    public:
+        enum kind { NULL_JST = 0, DECISION, CLAUSE, LAZY };
+        justification():m_data(TAG(void *, static_cast<void*>(0), NULL_JST)) { SASSERT(is_null()); }
+        justification(bool):m_data(TAG(void *, static_cast<void*>(0), DECISION)) { SASSERT(is_decision()); }
+        justification(clause * c):m_data(TAG(void *, c, CLAUSE)) { SASSERT(is_clause()); }
+        justification(lazy_justification * j):m_data(TAG(void *, j, LAZY)) { SASSERT(is_lazy()); }
+        kind get_kind() const { return static_cast<kind>(GET_TAG(m_data)); }
+        bool is_null() const { return get_kind() == NULL_JST; }
+        bool is_decision() const { return get_kind() == DECISION; }
+        bool is_clause() const { return get_kind() == CLAUSE; }
+        bool is_lazy() const { return get_kind() == LAZY; }
+        clause * get_clause() const { return UNTAG(clause*, m_data); }
+        lazy_justification * get_lazy() const { return UNTAG(lazy_justification*, m_data); }
+        bool operator==(justification other) const { return m_data == other.m_data; }
+        bool operator!=(justification other) const { return m_data != other.m_data; }
+    };
+    
+    const justification null_justification;
+    const justification decided_justification(true);
+
+    inline justification mk_clause_jst(clause const * c) { return justification(const_cast<clause*>(c)); }
+    inline justification mk_lazy_jst(small_object_allocator & a, unsigned num, literal const * lits) {
+        void * mem = a.allocate(lazy_justification::get_obj_size(num));
+        return justification(new (mem) lazy_justification(num, lits));
+    }
+
+    inline void del_jst(small_object_allocator & a, justification jst) {
+        if (jst.is_lazy()) {
+            lazy_justification * ptr = jst.get_lazy();
+            unsigned obj_sz = lazy_justification::get_obj_size(ptr->size());
+            a.deallocate(obj_sz, ptr);
+        }
+    }
+};
+
+#endif

src/nlsat/nlsat_scoped_literal_vector.h

@@ -0,0 +1,89 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_scoped_literal_vector.h
+
+Abstract:
+
+    Scoped vector for nlsat literals.
+    Just to be "cancel" safe.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-13.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_SCOPED_LITERAL_VECTOR_H_
+#define _NLSAT_SCOPED_LITERAL_VECTOR_H_
+
+#include"nlsat_solver.h"
+
+namespace nlsat {
+
+    class scoped_literal_vector {
+        solver &       m_solver;
+        literal_vector m_lits;
+    public:
+        scoped_literal_vector(solver & s):m_solver(s) {}
+        ~scoped_literal_vector() { reset(); }
+        unsigned size() const { return m_lits.size(); }
+        bool empty() const { return m_lits.empty(); }
+        literal operator[](unsigned i) const { return m_lits[i]; }
+        void reset() {
+            unsigned sz = m_lits.size();
+            for (unsigned i = 0; i < sz; i++) {
+                m_solver.dec_ref(m_lits[i]);
+            }
+            m_lits.reset();
+        }
+        void push_back(literal l) {
+            m_solver.inc_ref(l);
+            m_lits.push_back(l);
+        }
+        void set(unsigned i, literal l) {
+            m_solver.inc_ref(l);
+            m_solver.dec_ref(m_lits[i]);
+            m_lits[i] = l;
+        }
+        literal const * c_ptr() const { return m_lits.c_ptr(); }
+        void shrink(unsigned new_sz) {
+            SASSERT(new_sz <= m_lits.size());
+            unsigned sz = m_lits.size();
+            if (new_sz == sz)
+                return;
+            for (unsigned i = new_sz; i < sz; i++) {
+                m_solver.dec_ref(m_lits[i]);
+            }
+            m_lits.shrink(new_sz);
+        }
+        void append(unsigned sz, literal const * ls) {
+            for (unsigned i = 0; i < sz; i++)
+                push_back(ls[i]);
+        }
+    };
+
+
+    class scoped_literal {
+        solver &  m_solver;
+        literal   m_lit;
+    public:
+        scoped_literal(solver & s):m_solver(s), m_lit(null_literal) {}
+        ~scoped_literal() { m_solver.dec_ref(m_lit); }
+        scoped_literal & operator=(literal l) { 
+            m_solver.inc_ref(l);
+            m_solver.dec_ref(m_lit);
+            m_lit = l;
+            return *this;
+        }
+        scoped_literal & operator=(scoped_literal const & l) { return operator=(l.m_lit); }
+        operator literal&() { return m_lit; }
+        operator literal const &() const { return m_lit; }
+        void neg() { m_lit.neg(); }
+    };
+};
+
+#endif

src/nlsat/nlsat_solver.h

@@ -0,0 +1,187 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_solver.h
+
+Abstract:
+
+    Nonlinear arithmetic satisfiability procedure.  The procedure is
+    complete for nonlinear real arithmetic, but it also has limited
+    support for integers.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-02.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_SOLVER_H_
+#define _NLSAT_SOLVER_H_
+
+#include"nlsat_types.h"
+#include"params.h"
+#include"statistics.h"
+
+namespace nlsat {
+
+    class solver {
+        struct imp;
+        imp * m_imp;
+    public:
+        solver(params_ref const & p);
+        ~solver();
+
+        /**
+           \brief Return reference to rational manager.
+        */
+        unsynch_mpq_manager & qm();
+        
+        /**
+           \brief Return reference to algebraic number manager.
+        */
+        anum_manager & am();
+
+        /**
+           \brief Return a reference to the polynomial manager used by the solver.
+        */
+        pmanager & pm();
+
+        void set_display_var(display_var_proc const & proc);
+
+        // -----------------------
+        //
+        // Variable, Atoms, Clauses & Assumption creation
+        //
+        // -----------------------
+
+        /**
+           \brief Create a fresh boolean variable that is not associated with any 
+                  nonlinear arithmetic atom.
+        */
+        bool_var mk_bool_var();
+    
+        /**
+           \brief Create a real/integer variable.
+        */
+        var mk_var(bool is_int);
+        
+        /**
+           \brief Create an atom of the form: p=0, p<0, p>0
+           Where p = ps[0]^e[0]*...*ps[sz-1]^e[sz-1]
+
+           e[i] = 1 if is_even[i] is false
+           e[i] = 2 if is_even[i] is true
+
+           \pre sz > 0
+        */
+        bool_var mk_ineq_atom(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even);
+
+        /**
+           \brief Create a literal for the p=0, p<0, p>0.
+           Where p = ps[0]^e[0]*...*ps[sz-1]^e[sz-1]  for sz > 0
+                 p = 1                                for sz = 0
+
+           e[i] = 1 if is_even[i] is false
+           e[i] = 2 if is_even[i] is true
+        */
+        literal mk_ineq_literal(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even);
+
+        /**
+           \brief Create an atom of the form: x=root[i](p), x<root[i](p), x>root[i](p)
+        */
+        bool_var mk_root_atom(atom::kind k, var x, unsigned i, poly * p);
+
+        void inc_ref(bool_var b);
+        void inc_ref(literal l) { inc_ref(l.var()); }
+        void dec_ref(bool_var b);
+        void dec_ref(literal l) { dec_ref(l.var()); }
+        
+        /**
+           \brief Create a new clause.
+        */
+        void mk_clause(unsigned num_lits, literal * lits, assumption a = 0);
+
+        // -----------------------
+        //
+        // Basic
+        //
+        // -----------------------
+
+        /**
+           \brief Return the number of Boolean variables.
+        */
+        unsigned num_bool_vars() const;
+        
+        /**
+           \brief Get atom associated with Boolean variable. Return 0 if there is none.
+        */
+        atom * bool_var2atom(bool_var b);
+
+        /**
+           \brief Return number of integer/real variables
+        */
+        unsigned num_vars() const;
+
+        bool is_int(var x) const;
+
+        // -----------------------
+        //
+        // Search
+        //
+        // -----------------------
+        lbool check();
+
+        // -----------------------
+        //
+        // Model
+        //
+        // -----------------------
+        anum const & value(var x) const;
+        lbool bvalue(bool_var b) const;
+        bool is_interpreted(bool_var b) const;
+
+        lbool value(literal l) const;
+
+        // -----------------------
+        //
+        // Misc
+        //
+        // -----------------------
+        void updt_params(params_ref const & p);
+        static void collect_param_descrs(param_descrs & d);
+
+        void set_cancel(bool f);
+        void collect_statistics(statistics & st);
+        void reset_statistics();
+        void display_status(std::ostream & out) const;
+
+        // -----------------------
+        //
+        // Pretty printing
+        //
+        // -----------------------
+   
+        /**
+           \brief Display solver's state.
+        */
+        void display(std::ostream & out) const;
+
+        /**
+           \brief Display literal
+        */
+        void display(std::ostream & out, literal l) const;
+
+        /**
+           \brief Display variable
+        */
+        void display(std::ostream & out, var x) const;
+        
+        display_var_proc const & display_proc() const;
+    };
+
+};
+
+#endif

src/nlsat/nlsat_types.cpp

@@ -0,0 +1,70 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_types.cpp
+
+Abstract:
+
+    Basic types used in the nonlinear arithmetic satisfiability procedure.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-02.
+
+Revision History:
+
+--*/
+#include"nlsat_types.h"
+#include"debug.h"
+#include"hash.h"
+#include"polynomial.h"
+ 
+namespace nlsat {
+
+    ineq_atom::ineq_atom(kind k, unsigned sz, poly * const * ps, bool const * is_even, var max_var):
+        atom(k, max_var),
+        m_size(sz) {
+        for (unsigned i = 0; i < m_size; i++) {
+            m_ps[i] = TAG(poly *, ps[i], is_even[i] ? 1 : 0);
+        }
+        SASSERT(is_ineq_atom());
+    }
+
+    unsigned ineq_atom::hash_proc::operator()(ineq_atom const * a) const {
+        return get_composite_hash<ineq_atom const *, ineq_atom::khasher, ineq_atom::chasher>(a, a->m_size);
+    }
+
+    bool ineq_atom::eq_proc::operator()(ineq_atom const * a1, ineq_atom const * a2) const {
+        if (a1->m_size != a2->m_size || a1->m_kind != a2->m_kind)
+            return false;
+        unsigned sz = a1->m_size;
+        for (unsigned i = 0; i < sz; i++) {
+            if (a1->m_ps[i] != a2->m_ps[i])
+                return false;
+        }
+        return true;
+    }
+
+    root_atom::root_atom(kind k, var x, unsigned i, poly * p):
+        atom(k, x),
+        m_x(x),
+        m_i(i),
+        m_p(p) {
+        SASSERT(is_root_atom());
+    }
+
+    unsigned root_atom::hash_proc::operator()(root_atom const * a) const {
+        unsigned _a = a->m_x;
+        unsigned _b = ((a->m_i << 2) | (a->m_kind));
+        unsigned _c = polynomial::manager::id(a->m_p);
+        mix(_a, _b, _c);
+        return _c;
+    }
+
+    bool root_atom::eq_proc::operator()(root_atom const * a1, root_atom const * a2) const {
+        return a1->m_kind == a2->m_kind && a1->m_x == a2->m_x && a1->m_i == a2->m_i && a1->m_p == a2->m_p;
+    }
+
+};

src/nlsat/nlsat_types.h

@@ -0,0 +1,150 @@
+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    nlsat_types.h
+
+Abstract:
+
+    Basic types used in the nonlinear arithmetic satisfiability procedure.
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-01-02.
+
+Revision History:
+
+--*/
+#ifndef _NLSAT_TYPES_H_
+#define _NLSAT_TYPES_H_
+
+#include"polynomial.h"
+#include"buffer.h"
+#include"sat_types.h"
+#include"z3_exception.h"
+
+namespace algebraic_numbers {
+    class anum;
+    class manager;
+};
+
+namespace nlsat {
+#define NLSAT_VB_LVL 10
+    typedef void *                     assumption;
+    typedef void *                     assumption_set;
+
+    typedef sat::bool_var              bool_var;
+    typedef sat::bool_var_vector       bool_var_vector;
+    const bool_var null_bool_var = sat::null_bool_var;
+    typedef sat::literal               literal;
+    const literal null_literal = sat::null_literal;
+    typedef sat::literal_vector        literal_vector;
+
+    inline literal to_literal(unsigned x) { return sat::to_literal(x); }
+    
+    typedef polynomial::var            var;
+    typedef polynomial::var_vector     var_vector;
+    typedef polynomial::manager        pmanager;
+    typedef polynomial::polynomial     poly;
+    const var null_var = polynomial::null_var;
+    
+    const var true_bool_var = 0;
+    const literal true_literal(true_bool_var, false);
+    const literal false_literal(true_bool_var, true);
+
+    typedef polynomial::display_var_proc display_var_proc;
+
+    class atom;
+    class ineq_atom; // atoms of the form: p=0, p>0, p<0,  where p is a product of polynomials
+    class root_atom; // atoms of the form: x=root[i](p), x<root[i](p), x>root[i](p), where x is a variable and p a polynomial.
+    
+    class clause; 
+    class solver;
+
+    class atom {
+    public:
+        enum kind { EQ=0, LT, GT, ROOT_EQ=10, ROOT_LT, ROOT_GT, ROOT_LE, ROOT_GE };
+        static kind flip(kind k) {
+            SASSERT(k == EQ || k == LT || k == GT);
+            if (k == LT) return GT;
+            if (k == GT) return LT;
+            return EQ;
+        }
+    protected:
+        friend class solver;
+        kind     m_kind;
+        unsigned m_ref_count;
+        bool_var m_bool_var;
+        var      m_max_var;
+    public:
+        atom(kind k, var max_var):m_kind(k), m_ref_count(0), m_bool_var(null_bool_var), m_max_var(max_var) {}
+        bool is_eq() const { return m_kind == EQ || m_kind == ROOT_EQ; }
+
+        bool is_ineq_atom() const { return m_kind <= GT; }
+        bool is_root_atom() const { return m_kind >= ROOT_EQ; }
+        kind get_kind() const { return m_kind; }
+        bool_var bvar() const { return m_bool_var; }
+        var max_var() const { return m_max_var; }
+        unsigned ref_count() const { return m_ref_count; }
+        void inc_ref() { m_ref_count++; }
+        void dec_ref() { SASSERT(m_ref_count > 0); m_ref_count--; }
+    };
+
+    typedef ptr_vector<atom> atom_vector;
+
+    class ineq_atom : public atom {
+        friend class solver;
+        unsigned     m_size;
+        poly *       m_ps[0];
+        ineq_atom(kind k, unsigned sz, poly * const * ps, bool const * is_even, var max_var);
+        static unsigned get_obj_size(unsigned sz) { return sizeof(ineq_atom) + sizeof(poly*)*sz; }
+    public:
+        unsigned size() const { return m_size; }
+        poly * p(unsigned i) const { SASSERT(i < size()); return UNTAG(poly*, m_ps[i]); }
+        // Return true if i-th factor has odd degree
+        bool is_odd(unsigned i) const { SASSERT(i < size()); return GET_TAG(m_ps[i]) == 0; }
+        bool is_even(unsigned i) const { return !is_odd(i); }
+        struct khasher { unsigned operator()(ineq_atom const * a) const { return a->m_kind; } };
+        struct chasher { unsigned operator()(ineq_atom const * a, unsigned idx) const { return polynomial::manager::id(a->p(idx)); } };
+        struct hash_proc { unsigned operator()(ineq_atom const * a) const; };
+        struct eq_proc { bool operator()(ineq_atom const * a1, ineq_atom const * a2) const; };
+    };
+
+    class root_atom : public atom {
+        friend class solver;
+        var      m_x;
+        unsigned m_i;
+        poly *   m_p;
+        root_atom(kind k, var x, unsigned i, poly * p);
+    public:
+        var x() const { return m_x; }
+        unsigned i() const { return m_i; }
+        poly * p() const { return m_p; }
+        struct hash_proc { unsigned operator()(root_atom const * a) const; };
+        struct eq_proc { bool operator()(root_atom const * a1, root_atom const * a2) const; };
+    };
+
+    inline ineq_atom * to_ineq_atom(atom * a) { SASSERT(a->is_ineq_atom()); return static_cast<ineq_atom*>(a); }
+    inline root_atom * to_root_atom(atom * a) { SASSERT(a->is_root_atom()); return static_cast<root_atom*>(a); }
+    inline ineq_atom const * to_ineq_atom(atom const * a) { SASSERT(a->is_ineq_atom()); return static_cast<ineq_atom const *>(a); }
+    inline root_atom const * to_root_atom(atom const * a) { SASSERT(a->is_root_atom()); return static_cast<root_atom const *>(a); }
+
+    typedef algebraic_numbers::anum    anum;
+    typedef algebraic_numbers::manager anum_manager;
+
+    class solver_exception : public default_exception {
+    public:
+        solver_exception(char const * msg):default_exception(msg) {}
+    };
+
+    class assignment;
+    
+    inline int normalize_sign(int s) {
+        if (s < 0)  return -1;
+        if (s == 0) return 0;
+        return 1;
+    }
+};
+
+#endif