Merge branch 'master' into polysat

src/math/dd/dd_pdd.h

@@ -413,7 +413,8 @@ namespace dd {
         bool is_var() const { return m.is_var(root); }
         /** Polynomial is of the form a * x + b for numerals a, b. */
         bool is_unilinear() const { return !is_val() && lo().is_val() && hi().is_val(); }
-        bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); } 
+        bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); }
+        bool is_offset() const { return !is_val() && lo().is_val() && hi().is_one(); }
         bool is_binary() const { return m.is_binary(root); }
         bool is_monomial() const { return m.is_monomial(root); }
         bool is_univariate() const { return m.is_univariate(root); }

src/math/grobner/grobner.cpp

@@ -116,10 +116,13 @@ void grobner::reset() {
 }
 
 void grobner::display_var(std::ostream & out, expr * var) const {
+    out << "#" << var->get_id();
+#if 0
     if (is_app(var) && to_app(var)->get_num_args() > 0)
         out << mk_bounded_pp(var, m_manager);
     else
         out << mk_pp(var, m_manager);
+#endif    
 }
 
 void grobner::display_vars(std::ostream & out, unsigned num_vars, expr * const * vars) const {

src/math/grobner/pdd_solver.cpp

@@ -160,6 +160,9 @@ namespace dd {
         }
         while (simplified && !eq.poly().is_val());
 
+        if (eq.poly().is_unary() && eq.poly().hi().val() < 0)
+            eq = -eq.poly();
+
         TRACE("dd.solver", display(tout << "simplification result: ", eq););
     }
 

src/math/lp/explanation.h

@@ -53,7 +53,8 @@ public:
         if (e.m_vector.empty()) {
             for (constraint_index j : e.m_set)
                 push_back(j);
-        } else {
+        } 
+        else {
             for (const auto & p : e.m_vector) {
                 add_pair(p.first, p.second);
             }
@@ -71,15 +72,17 @@ public:
         constraint_index ci() const { return m_var; }
         const mpq &coeff() const { return m_coeff; }
     };
+
     class iterator {
         bool      m_run_on_vector;
+        mpq       m_one = one_of_type<mpq>();
         pair_vec::const_iterator    m_vi;
         ci_set::iterator            m_ci;
     public:
         cimpq operator*() const {
             return m_run_on_vector?
                 cimpq( m_vi->first, m_vi->second) :
-                cimpq( *m_ci, one_of_type<mpq>()); 
+                cimpq( *m_ci, m_one); 
         }        
         iterator operator++() {
             if (m_run_on_vector)

src/math/lp/gomory.cpp

@@ -82,7 +82,8 @@ class create_cut {
         TRACE("gomory_cut_detail", tout << "new_a = " << new_a << ", k = " << m_k << ", lcm_den = " << m_lcm_den << "\n";);
 #if SMALL_CUTS
         // if (numerator(new_a).is_big()) throw found_big(); 
-        if (numerator(new_a) > m_big_number) throw found_big(); 
+        if (numerator(new_a) > m_big_number)
+            throw found_big(); 
 #endif
     }
 
@@ -90,28 +91,24 @@ class create_cut {
         TRACE("gomory_cut_detail_real", tout << "j = " << j << ", a = " << a << ", m_k = " << m_k << "\n";);
         mpq new_a;
         if (at_lower(j)) {
-            if (a.is_pos()) {
+            if (a.is_pos()) 
                 // the delta is a (x - f) is positive it has to grow and fight m_one_minus_f
                 new_a = a / m_one_minus_f;
-            }
-            else {
+            else 
                 // the delta is negative and it works again m_f
                 new_a = - a / m_f;
-            }
             m_k.addmul(new_a, lower_bound(j).x); // is it a faster operation than
             // k += lower_bound(j).x * new_a;  
             m_ex->push_back(column_lower_bound_constraint(j));
         }
         else {
             lp_assert(at_upper(j));
-            if (a.is_pos()) {
+            if (a.is_pos()) 
                 // the delta is works again m_f
                 new_a =  - a / m_f; 
-            }
-            else {
+            else 
                 // the delta is positive works again m_one_minus_f
                 new_a =   a / m_one_minus_f; 
-            }
             m_k.addmul(new_a, upper_bound(j).x); //  k += upper_bound(j).x * new_a; 
             m_ex->push_back(column_upper_bound_constraint(j));
         }
@@ -121,7 +118,8 @@ class create_cut {
         
 #if SMALL_CUTS
         // if (numerator(new_a).is_big()) throw found_big(); 
-        if (numerator(new_a) > m_big_number) throw found_big(); 
+        if (numerator(new_a) > m_big_number)
+            throw found_big(); 
 #endif
     }
 
@@ -146,13 +144,15 @@ class create_cut {
                 if (!m_k.is_int())
                     m_k = ceil(m_k);
                 m_t.add_monomial(mpq(1), v);
-            } else {
+            }
+            else {
                 m_k /= -a;
                 if (!m_k.is_int())
                     m_k = ceil(m_k);
                 m_t.add_monomial(-mpq(1), v);
             }
-        } else {
+        }
+        else {
             m_lcm_den = lcm(m_lcm_den, denominator(m_k));
             lp_assert(m_lcm_den.is_pos());
             TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;);
@@ -176,14 +176,12 @@ class create_cut {
     }
 
     std::ostream& dump_coeff_val(std::ostream & out, const mpq & a) const {
-        if (a.is_int()) {
+        if (a.is_int()) 
             out << a;
-        } 
         else if ( a >= zero_of_type<mpq>())
             out << "(/ " << numerator(a) << " " << denominator(a) << ")";
-        else {
-            out << "(- ( / " <<   numerator(-a) << " " << denominator(-a) << "))";
-        }
+        else 
+            out << "(- (/ " <<   numerator(-a) << " " << denominator(-a) << "))";
         return out;
     }
     
@@ -197,12 +195,10 @@ class create_cut {
     
     std::ostream& dump_row_coefficients(std::ostream & out) const {
         mpq lc(1);
-        for (const auto& p : m_row) {
+        for (const auto& p : m_row) 
             lc = lcm(lc, denominator(p.coeff())); 
-        }        
-        for (const auto& p : m_row) {
+        for (const auto& p : m_row) 
             dump_coeff_val(out << " (* ", p.coeff()*lc) << " " << var_name(p.var()) << ")";
-        }
         return out;
     }
     
@@ -218,9 +214,8 @@ class create_cut {
     
     void dump_declarations(std::ostream& out) const {
         // for a column j the var name is vj
-        for (const auto & p : m_row) {
+        for (const auto & p : m_row) 
             dump_declaration(out, p.var());
-        }
         for (lar_term::ival p : m_t) {
             auto t = lia.lra.column2tv(p.column());
             if (t.is_term()) {
@@ -239,12 +234,11 @@ class create_cut {
     void dump_explanations(std::ostream& out) const {
         for (const auto & p : m_row) {            
             unsigned j = p.var();
-            if (j == m_inf_col || (!is_real(j) && p.coeff().is_int())) {
+            if (j == m_inf_col || (!is_real(j) && p.coeff().is_int())) 
                 continue;
-            }
-            else if (at_lower(j)) {
+            else if (at_lower(j)) 
                 dump_lower_bound_expl(out, j);
-            } else {
+            else {
                 lp_assert(at_upper(j));
                 dump_upper_bound_expl(out, j);
             }
@@ -252,9 +246,8 @@ class create_cut {
     }
 
     std::ostream& dump_term_coefficients(std::ostream & out) const {
-        for (lar_term::ival p : m_t) {
+        for (lar_term::ival p : m_t) 
             dump_coeff(out, p);
-        }
         return out;
     }
     
@@ -281,9 +274,8 @@ class create_cut {
 public:
     void dump(std::ostream& out) {
         out << "applying cut at:\n"; print_linear_combination_indices_only<row_strip<mpq>, mpq>(m_row, out); out << std::endl;
-        for (auto & p : m_row) {
+        for (auto & p : m_row) 
             lia.lra.print_column_info(p.var(), out);
-        }
         out << "inf_col = " << m_inf_col << std::endl;
     }
 
@@ -304,7 +296,8 @@ public:
         m_abs_max = 0;
         for (const auto & p : m_row) {
             mpq t = abs(ceil(p.coeff()));
-            if (t > m_abs_max) m_abs_max = t;
+            if (t > m_abs_max)
+                m_abs_max = t;
         }
         m_big_number = m_abs_max.expt(2);
 #endif
@@ -324,9 +317,8 @@ public:
                     m_ex->push_back(column_upper_bound_constraint(j));
                     continue;
                 }
-                if (is_real(j)) {  
+                if (is_real(j))   
                     real_case_in_gomory_cut(- p.coeff(), j);
-                } 
                 else if (!p.coeff().is_int()) {
                     some_int_columns = true;
                     m_fj = fractional_part(-p.coeff());
@@ -411,7 +403,8 @@ int gomory::find_basic_var() {
 lia_move gomory::operator()() {
     lra.move_non_basic_columns_to_bounds(true);
     int j = find_basic_var();
-    if (j == -1) return lia_move::undef;
+    if (j == -1)
+        return lia_move::undef;
     unsigned r = lia.row_of_basic_column(j);
     const row_strip<mpq>& row = lra.get_row(r);
     SASSERT(lra.row_is_correct(r));

src/math/lp/lp_settings.h

@@ -24,6 +24,7 @@ Revision History:
 #include <algorithm>
 #include <limits>
 #include <iomanip>
+#include <iostream>
 #include <cstring>
 #include "util/stopwatch.h"
 #include "util/statistics.h"

src/math/lp/lu.h

@@ -33,7 +33,7 @@ Revision History:
 #include "math/lp/static_matrix.h"
 #include <string>
 #include "math/lp/numeric_pair.h"
-#include <iostream>
+#include <ostream>
 #include <fstream>
 #include "math/lp/row_eta_matrix.h"
 #include "math/lp/square_dense_submatrix.h"

src/math/lp/mps_reader.h

@@ -26,7 +26,7 @@ Revision History:
 #include <string>
 #include "util/vector.h"
 #include <unordered_map>
-#include <iostream>
+#include <ostream>
 #include <fstream>
 #include <locale>
 #include "math/lp/lp_primal_simplex.h"

src/math/polynomial/algebraic_numbers.cpp

@@ -1466,7 +1466,10 @@ namespace algebraic_numbers {
                 qm().add(il, nbv, il);
                 qm().add(iu, nbv, iu);
                 // (il, iu) is an isolating refinable (rational) interval for the new polynomial.
-                upm().convert_q2bq_interval(m_add_tmp.size(), m_add_tmp.data(), il, iu, bqm(), l, u);
+                if (!upm().convert_q2bq_interval(m_add_tmp.size(), m_add_tmp.data(), il, iu, bqm(), l, u)) {
+                    TRACE("algebraic", tout << "conversion failed\n");
+                }
+                    
             }
             TRACE("algebraic",
                   upm().display(tout, m_add_tmp.size(), m_add_tmp.data());
@@ -1576,7 +1579,9 @@ namespace algebraic_numbers {
                 if (is_neg)
                     qm().swap(il, iu);
                 // (il, iu) is an isolating refinable (rational) interval for the new polynomial.
-                upm().convert_q2bq_interval(mulp.size(), mulp.data(), il, iu, bqm(), l, u);
+                if (!upm().convert_q2bq_interval(mulp.size(), mulp.data(), il, iu, bqm(), l, u)) {
+                    TRACE("algebraic", tout << "conversion failed\n");
+                }
             }
             TRACE("algebraic",
                   upm().display(tout, mulp.size(), mulp.data());
@@ -1690,7 +1695,10 @@ namespace algebraic_numbers {
                 qm().swap(inv_lower, inv_upper);
                 TRACE("algebraic_bug", tout << "inv new_bounds: " << qm().to_string(inv_lower) << ", " << qm().to_string(inv_upper) << "\n";);
                 // convert isolating interval back as a binary rational bound
-                upm().convert_q2bq_interval(cell_a->m_p_sz, cell_a->m_p, inv_lower, inv_upper, bqm(), lower(cell_a), upper(cell_a));
+                if (!upm().convert_q2bq_interval(cell_a->m_p_sz, cell_a->m_p, inv_lower, inv_upper, bqm(), lower(cell_a), upper(cell_a))) {
+                    TRACE("algebraic_bug", tout << "root isolation failed\n");                    
+                    throw algebraic_exception("inversion of algebraic number failed");
+                }
                 TRACE("algebraic_bug", tout << "after inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
                 update_sign_lower(cell_a);
                 SASSERT(acell_inv(*cell_a));       

src/math/realclosure/realclosure.cpp

@@ -29,6 +29,7 @@ Notes:
 #include "util/ref_vector.h"
 #include "util/ref_buffer.h"
 #include "util/common_msgs.h"
+#include <iostream>
 
 #ifndef REALCLOSURE_INI_BUFFER_SIZE
 #define REALCLOSURE_INI_BUFFER_SIZE 32

src/math/simplex/model_based_opt.cpp

@@ -485,6 +485,26 @@ namespace opt {
         }
     }
 
+    model_based_opt::row& model_based_opt::row::normalize() {
+#if 0
+        if (m_type == t_mod)
+            return *this;
+        rational D(denominator(abs(m_coeff)));
+        if (D == 0)
+            D = 1;
+        for (auto const& [id, coeff] : m_vars)
+            if (coeff != 0)
+                D = lcm(D, denominator(abs(coeff)));
+        if (D == 1)
+            return *this;
+        SASSERT(D > 0);
+        for (auto & [id, coeff] : m_vars)
+            coeff *= D;
+        m_coeff *= D;
+#endif
+        return *this;
+    }
+
     // 
     // Let
     //   row1: t1 + a1*x <= 0
@@ -923,9 +943,9 @@ namespace opt {
     }
 
     void model_based_opt::get_live_rows(vector<row>& rows) {
-        for (row const& r : m_rows) {
+        for (row & r : m_rows) {
             if (r.m_alive) {
-                rows.push_back(r);
+                rows.push_back(r.normalize());
             }
         }
     }

src/math/simplex/model_based_opt.h

@@ -59,6 +59,7 @@ namespace opt {
             bool        m_alive;        // rows can be marked dead if they have been processed.
             void reset() { m_vars.reset(); m_coeff.reset(); m_value.reset(); }
 
+            row& normalize();
             void neg() { for (var & v : m_vars) v.m_coeff.neg(); m_coeff.neg(); m_value.neg(); }
             rational get_coefficient(unsigned x) const;
         };

src/math/simplex/simplex.cpp

@@ -18,6 +18,7 @@ Notes:
 --*/
 
 #include "math/simplex/simplex.h"
+#include "math/simplex/sparse_matrix_ops.h"
 #include "math/simplex/sparse_matrix_def.h"
 #include "math/simplex/simplex_def.h"
 #include "util/rational.h"
@@ -36,6 +37,9 @@ namespace simplex {
         }
     }
 
+    void kernel(sparse_matrix<mpq_ext>& M, vector<vector<rational>>& K) {
+        sparse_matrix_ops::kernel(M, K);
+    }
 
     void ensure_rational_solution(simplex<mpq_ext>& S) {
         rational delta(1);

src/math/simplex/simplex.h

@@ -33,6 +33,7 @@ Notes:
 
 #include "math/simplex/sparse_matrix.h"
 #include "util/mpq_inf.h"
+#include "util/rational.h"
 #include "util/heap.h"
 #include "util/lbool.h"
 #include "util/uint_set.h"
@@ -200,5 +201,7 @@ namespace simplex {
     };
 
     void ensure_rational_solution(simplex<mpq_ext>& s);
+
+    void kernel(sparse_matrix<mpq_ext>& s, vector<vector<rational>>& K);
 };
 

src/math/simplex/sparse_matrix.h

@@ -143,6 +143,7 @@ namespace simplex {
         svector<int>            m_var_pos;          // temporary map from variables to positions in row
         unsigned_vector         m_var_pos_idx;      // indices in m_var_pos
         stats                   m_stats;
+        scoped_numeral          m_zero;
 
         bool well_formed_row(unsigned row_id) const;
         bool well_formed_column(unsigned column_id) const;
@@ -151,7 +152,7 @@ namespace simplex {
 
     public:
 
-        sparse_matrix(manager& _m): m(_m) {}
+        sparse_matrix(manager& _m): m(_m), m_zero(m) {}
         ~sparse_matrix();
         void reset();
         
@@ -208,41 +209,40 @@ namespace simplex {
         row_iterator row_begin(row const& r) { return row_iterator(m_rows[r.id()], true); }
         row_iterator row_end(row const& r) { return row_iterator(m_rows[r.id()], false); }
 
-        struct row_entries_t {
-            sparse_matrix& m;
-            row const& r;
-            row_entries_t(sparse_matrix & m, row const& r): m(m), r(r) {}
-            row_iterator begin() { return m.row_begin(r); }
-            row_iterator end() { return m.row_end(r); }
+        class row_vars {
+            friend class sparse_matrix;
+            sparse_matrix& s;
+            row r;
+            row_vars(sparse_matrix& s, row r): s(s), r(r) {}
+        public:
+            row_iterator begin() { return s.row_begin(r); }
+            row_iterator end() { return s.row_end(r); }
         };
 
-        row_entries_t row_entries(row const& r) { return row_entries_t(*this, r); }
+        row_vars get_row(row r) { return row_vars(*this, r); }
 
         unsigned column_size(var_t v) const { return m_columns[v].size(); }
 
+        unsigned num_vars() const { return m_columns.size(); }
+        unsigned num_rows() const { return m_rows.size(); }
+
         class col_iterator {
             friend class sparse_matrix;
             unsigned             m_curr;
-            int                  m_var;
-            sparse_matrix const& m_sm;
-
-            column const& col() const { 
-                return m_sm.m_columns[m_var]; 
-            }
-
+            column const&        m_col;
+            vector<_row>&        m_rows;
             void move_to_used() {
                 while (m_curr < col().num_entries() && col().m_entries[m_curr].is_dead()) {
                     ++m_curr;
                 }
             }
-
-            col_iterator(int var, sparse_matrix const& sm, bool begin): 
-                m_curr(0), m_var(var), m_sm(sm) {
-                ++col().m_refs;
-                if (begin) 
+            col_iterator(column const& c, vector<_row>& r, bool begin): 
+                m_curr(0), m_col(c), m_rows(r) {
+                ++m_col.m_refs;
+                if (begin)
                     move_to_used();
                 else 
-                    m_curr = col().num_entries();
+                    m_curr = m_col.num_entries();
             }
 
         public:
@@ -255,18 +255,20 @@ namespace simplex {
                 m_curr(other.m_curr),
                 m_var(other.m_var),
                 m_sm(other.m_sm) {
-                ++col().m_refs;
+                ++m_col.m_refs;
             }
 
             row get_row() const { 
-                return row(col().m_entries[m_curr].m_row_id); 
+                return row(m_col.m_entries[m_curr].m_row_id);
             }
-            row_entry const& get_row_entry() const {
-                col_entry const& c = col().m_entries[m_curr];
+
+            row_entry& get_row_entry() {
+                col_entry const& c = m_col.m_entries[m_curr];
                 int row_id = c.m_row_id;
                 return m_sm.m_rows[row_id].m_entries[c.m_row_idx];
             }
-            
+
+            std::pair<row, row_entry*> operator*() { return std::make_pair(get_row(), &get_row_entry()); }
             col_iterator & operator++() { ++m_curr; move_to_used(); return *this; }
             col_iterator operator++(int) { col_iterator tmp = *this; ++*this; return tmp; }
             bool operator==(col_iterator const & it) const { return m_curr == it.m_curr; }
@@ -274,24 +276,67 @@ namespace simplex {
             col_iterator& operator*() { return *this; }
         };
 
-        col_iterator col_begin(int v) const { return col_iterator(v, *this, true); }
-        col_iterator col_end(int v) const { return col_iterator(v, *this, false); }
+        col_iterator col_begin(int v) { return col_iterator(m_columns[v], m_rows, true); }
+        col_iterator col_end(int v) { return col_iterator(m_columns[v], m_rows, false); }
 
-        class col_entries_t {
-            sparse_matrix const& m;
+        class var_rows {
+            friend class sparse_matrix;
+            sparse_matrix& s;
             int v;
+            var_rows(sparse_matrix& s, int v):s(s), v(v) {}
         public:
-            col_entries_t(sparse_matrix const& m, int v): m(m), v(v) {}
-            col_iterator begin() { return m.col_begin(v); }
-            col_iterator end() { return m.col_end(v); }
+            col_iterator begin() { return s.col_begin(v); }
+            col_iterator end() { return s.col_end(v); }
+        };
+        
+        var_rows get_rows(int v) { return var_rows(*this, v); }
+
+        class all_row_iterator {
+            friend class sparse_matrix;
+            unsigned m_curr;
+            vector<_row> const& m_rows;
+            void move_to_next() {
+                while (m_curr < m_rows.size() && m_rows[m_curr].size() == 0) {
+                    //std::cout << "size is 0 for " << m_curr << "\n";
+                    ++m_curr;
+                }
+            }
+        public:
+            all_row_iterator(unsigned curr, vector<_row> const& rows): m_curr(curr), m_rows(rows) {
+                move_to_next();
+            }
+            row operator*() { return row(m_curr); }
+            all_row_iterator & operator++() { m_curr++; move_to_next(); return *this; }
+            all_row_iterator operator++(int) { all_row_iterator tmp = *this; ++*this; return tmp; }
+            bool operator==(all_row_iterator const& it) const { return m_curr == it.m_curr; }
+            bool operator!=(all_row_iterator const& it) const { return m_curr != it.m_curr; }
+        };
+        
+        class all_rows {
+            friend class sparse_matrix;
+            sparse_matrix& s;
+            all_rows(sparse_matrix& s): s(s) {}
+        public:
+            all_row_iterator begin() { return all_row_iterator(0, s.m_rows); }
+            all_row_iterator end() { return all_row_iterator(s.m_rows.size(), s.m_rows); }
         };
 
-        col_entries_t col_entries(int v) { return col_entries_t(*this, v); }
+
+        all_rows get_rows() { return all_rows(*this); }
+
+        numeral const& get_coeff(row r, unsigned v) {
+            for (auto & row : get_row(r)) 
+                if (row.m_var == v) 
+                    return row.m_coeff;
+            return m_zero;
+        }
 
         void display(std::ostream& out);
         void display_row(std::ostream& out, row const& r) const;
         bool well_formed() const;
 
+        manager& get_manager() { return m; }
+
         void collect_statistics(::statistics & st) const;
 
     };
@@ -312,4 +357,4 @@ namespace simplex {
         typedef unsynch_mpq_inf_manager eps_manager;
     };
 
-};
+}

src/math/simplex/sparse_matrix_def.h

@@ -298,6 +298,8 @@ namespace simplex {
 
     template<typename Ext>
     void sparse_matrix<Ext>::add_var(row dst, numeral const& n, var_t v) {
+        if (m.is_zero(n))
+            return;
         _row& r   = m_rows[dst.id()];
         column& c = m_columns[v];
         unsigned r_idx;
@@ -316,6 +318,9 @@ namespace simplex {
     */
     template<typename Ext>
     void sparse_matrix<Ext>::add(row row1, numeral const& n, row row2) {
+
+        if (m.is_zero(n))
+            return;
         m_stats.m_add_rows++;
         _row & r1 = m_rows[row1.id()];
         

src/math/simplex/sparse_matrix_ops.h

@@ -0,0 +1,85 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    sparse_matrix_ops.h
+
+Abstract:
+
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-01-15
+
+Notes:
+
+--*/
+
+#pragma once
+
+#include "math/simplex/sparse_matrix.h"
+#include "util/rational.h"
+
+namespace simplex {
+
+class sparse_matrix_ops {
+  public:
+    template <typename Ext>
+    static void kernel(sparse_matrix<Ext> &M, vector<vector<rational>> &K) {
+        using scoped_numeral = typename Ext::scoped_numeral;
+
+        vector<unsigned> d, c;
+        unsigned n_vars = M.num_vars(), n_rows = M.num_rows();
+        c.resize(n_rows, 0u);
+        d.resize(n_vars, 0u);
+
+        auto &m = M.get_manager();
+        scoped_numeral m_ik(m);
+        scoped_numeral D(m);
+
+        for (unsigned k = 0; k < n_vars; ++k) {
+            d[k] = 0;
+            for (auto [row, row_entry] : M.get_rows(k)) {
+                if (c[row.id()] != 0) continue;
+                auto &m_jk = row_entry->m_coeff;
+                if (mpq_manager<false>::is_zero(m_jk)) continue;
+
+                // D = rational(-1) / m_jk;
+                m.set(D, m_jk);
+                m.inv(D);
+                m.neg(D);
+
+                M.mul(row, D);
+                for (auto [row_i, row_i_entry] : M.get_rows(k)) {
+                    if (row_i.id() == row.id()) continue;
+                    m.set(m_ik, row_i_entry->m_coeff);
+                    // row_i += m_ik * row
+                    M.add(row_i, m_ik, row);
+                }
+                c[row.id()] = k + 1;
+                d[k] = row.id() + 1;
+                break;
+            }
+        }
+
+        for (unsigned k = 0; k < n_vars; ++k) {
+            if (d[k] != 0) continue;
+            K.push_back(vector<rational>());
+            for (unsigned i = 0; i < n_vars; ++i) {
+                if (d[i] > 0) {
+                    auto r = sparse_matrix<mpq_ext>::row(d[i] - 1);
+                    K.back().push_back(rational(M.get_coeff(r, k)));
+                } else if (i == k)
+                    K.back().push_back(rational(1));
+                else
+                    K.back().push_back(rational(0));
+            }
+        }
+    }
+
+    static void kernel(sparse_matrix<mpq_ext> &M, vector<vector<rational>> &K) {
+        kernel<mpq_ext>(M, K);
+    }
+};
+} // namespace simplex

src/math/subpaving/subpaving_t.h

@@ -18,7 +18,7 @@ Revision History:
 --*/
 #pragma once
 
-#include<iostream>
+#include<ostream>
 #include "util/tptr.h"
 #include "util/small_object_allocator.h"
 #include "util/chashtable.h"

src/math/subpaving/tactic/subpaving_tactic.cpp

@@ -26,6 +26,7 @@ Revision History:
 #include "util/mpff.h"
 #include "util/mpfx.h"
 #include "util/f2n.h"
+#include <iostream>
 
 class subpaving_tactic : public tactic {