From 179988e161ede103db1ed749eb90b82e880e69db Mon Sep 17 00:00:00 2001
From: Lev Nachmanson <levnach@hotmail.com>
Date: Wed, 5 May 2021 12:53:20 -0700
Subject: [PATCH] support recursive terms (#5246)

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
---
 src/math/lp/lar_solver.cpp | 24 +++++++++++++++++++++++-
 src/math/lp/lar_solver.h   |  2 ++
 src/math/lp/lar_term.h     | 14 ++++++++++++--
 3 files changed, 37 insertions(+), 3 deletions(-)

diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index fc4253e41..0837a7b5d 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -162,7 +162,7 @@ namespace lp {
                 continue;
             if (!m_terms[k]->contains(basis_j))
                 continue;
-            m_terms[k]->subst(basis_j, m_mpq_lar_core_solver.m_r_solver.m_pivot_row);
+            m_terms[k]->subst_in_row(basis_j, m_mpq_lar_core_solver.m_r_solver.m_pivot_row);
         }
     }
 
@@ -1717,6 +1717,27 @@ namespace lp {
         return true;
     }
 
+    void lar_solver::subst_known_terms(lar_term* t) {
+        std::set<unsigned> seen_terms;
+        for (const auto&p : *t) {
+            auto j = p.column();
+            if (this->column_corresponds_to_term(j)) {
+                seen_terms.insert(j);
+            }
+        }
+        while (!seen_terms.empty()) {
+            unsigned j = *seen_terms.begin();
+            seen_terms.erase(j);
+            auto tj = this->m_var_register.local_to_external(j);
+            auto& ot = this->get_term(tj);
+            for(const auto&  p : ot){
+              if (this->column_corresponds_to_term(p.column())) {
+                 seen_terms.insert(p.column());
+              }   
+            }
+            t->subst_by_term(ot, j);
+        }
+    }
     // do not register in m_var_register this term if ext_i == UINT_MAX
     var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>>& coeffs, unsigned ext_i) {
         TRACE("lar_solver_terms", print_linear_combination_of_column_indices_only(coeffs, tout) << ", ext_i =" << ext_i << "\n";);
@@ -1726,6 +1747,7 @@ namespace lp {
         if (strategy_is_undecided())
             return add_term_undecided(coeffs);
         lar_term* t = new lar_term(coeffs);
+        subst_known_terms(t);
         push_term(t);
         SASSERT(m_terms.size() == m_term_register.size());
         unsigned adjusted_term_index = m_terms.size() - 1;
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 452614caa..6e61354ff 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -584,6 +584,8 @@ public:
         }
         return out;
     }
+    
+    void subst_known_terms(lar_term*);
 
     inline std::ostream& print_column_bound_info(unsigned j, std::ostream& out) const {
         return m_mpq_lar_core_solver.m_r_solver.print_column_bound_info(j, out);
diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h
index bea2d24d3..a8ab37473 100644
--- a/src/math/lp/lar_term.h
+++ b/src/math/lp/lar_term.h
@@ -53,10 +53,20 @@ public:
     unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
 
     template <typename T>
-    const T & coeffs() const {
+    const T &   coeffs() const {
         return m_coeffs;
     }
     
+    void subst_by_term(const lar_term& t, unsigned term_column) {
+        auto it = this->m_coeffs.find_core(term_column);
+        if (it == nullptr) return;
+        mpq a = it->get_data().m_value;
+        this->m_coeffs.erase(term_column);
+        for (const auto & p : t) {
+            this->add_monomial(a * p.coeff(), p.column());
+        }
+    }
+
     lar_term(const vector<std::pair<mpq, unsigned>>& coeffs) {
         for (const auto & p : coeffs) {
             add_monomial(p.first, p.second);
@@ -94,7 +104,7 @@ public:
     }
 
     // j is the basic variable to substitute
-    void subst(unsigned j, indexed_vector<mpq> & li) {
+    void subst_in_row(unsigned j, indexed_vector<mpq> & li) {
         auto* it = m_coeffs.find_core(j);
         if (it == nullptr) return;
         const mpq & b = it->get_data().m_value;