diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp
index 1fedb4221..93dffc76f 100644
--- a/src/math/lp/dioph_eq.cpp
+++ b/src/math/lp/dioph_eq.cpp
@@ -127,7 +127,7 @@ namespace lp {
            of variables in X and of integer constants showing the substitutions
         */
         u_map<term_o> m_sigma;
-
+        
     public:
         int_solver& lia;
         lar_solver& lra;
@@ -135,7 +135,7 @@ namespace lp {
 
         // we start assigning with UINT_MAX and go down, print it as l(UINT_MAX - m_last_fresh_x_var)
         unsigned  m_last_fresh_x_var = UINT_MAX;
-
+        bool m_report_branch = false;
 
         // set F
         std::list<unsigned> m_f;  // F = {λ(t):t in m_f}
@@ -144,7 +144,6 @@ namespace lp {
         vector<unsigned> m_k2s; // k is substituted by using equation in m_eprime[m_k2s[k]]
 
         unsigned            m_conflict_index = -1;  // m_eprime[m_conflict_index] gives the conflict
-
         imp(int_solver& lia, lar_solver& lra): lia(lia), lra(lra) {}
 
         term_o row_to_term(const row_strip<mpq>& row) const {
@@ -160,11 +159,13 @@ namespace lp {
         }
 
         void init() {
+            m_report_branch = false;
             unsigned n_of_rows = lra.A_r().row_count();
             m_k2s.clear();
             m_k2s.resize(lra.column_count(), -1);
             m_conflict_index = -1;
             m_infeas_explanation.clear();
+            lia.get_term().clear();
 
             auto all_vars_are_int = [this](const auto& row) {
                 for (const auto& p : row) {
@@ -222,6 +223,8 @@ namespace lp {
             return lra.print_expl(out, ex);
         }
         // returns true if no conflict is found and false otherwise
+        // this function devides all cooficients, and the free constant, of the ep.m_e  by the gcd of all coefficients, 
+        // it is needed by the next steps
         bool normalize_e_by_gcd(eprime_pair& ep) {
             TRACE("dioph_eq", print_term_o(ep.m_e, tout << "m_e:") << std::endl;
                   print_dep(tout << "m_l:", ep.m_l) << std::endl;
@@ -235,7 +238,7 @@ namespace lp {
                 return true;
             TRACE("dioph_eq", tout << "g:" << g << std::endl;);
             mpq new_c = ep.m_e.c() / g;
-            if (new_c.is_int() == false) {
+            if (!new_c.is_int()) {
                 TRACE("dioph_eq",
                       print_term_o(ep.m_e, tout << "conflict m_e:") << std::endl;
                       tout << "g:" << g << std::endl;
@@ -246,6 +249,21 @@ namespace lp {
                           print_term_o(t.m_value, tout) << std::endl;
                       }
                     );
+                /* We have ep.m_e/g = 0, or sum((coff_i/g)*x_i) + new_c = 0,
+                or sum((coeff_i/g)*x_i) = -new_c, where new_c is not an integer
+                Then sum((coeff_i/g)*x_i) <= floor(-new_c) or sum((coeff_i/g)*x_i) >= ceil(-new_c)
+                */
+                if (lra.settings().stats().m_dio_conflicts % lra.settings().dio_cut_from_proof_period() == 0) {
+                    // prepare int_solver for reporting
+                    lar_term& t = lia.get_term();
+                    for (auto& p: ep.m_e) {
+                        t.add_monomial(p.coeff()/g, p.j());
+                    }
+                    lia.offset() = floor(-new_c);
+                    lia.is_upper() = true;
+                    m_report_branch = true;
+                    TRACE("dioph_eq", tout << "prepare branch:"; print_lar_term_L(t, tout) << " <= " << lia.offset() << std::endl;);
+                }
                 return false;
             } else {
                 for (auto& p: ep.m_e.coeffs()) {
@@ -467,13 +485,20 @@ namespace lp {
         lia_move check() {
             init();
             while(m_f.size()) {
-                if (!normalize_by_gcd())
+                if (!normalize_by_gcd()) {
+                    lra.settings().stats().m_dio_conflicts++;
+                    if (m_report_branch) {
+                        m_report_branch = false;
+                        return lia_move::branch;
+                    }
                     return lia_move::conflict;
+                }
                 rewrite_eqs();
             }
             TRACE("dioph_eq", print_S(tout););
             lia_move ret = tighten_with_S();
             if (ret == lia_move::conflict) {
+                lra.settings().stats().m_dio_conflicts++;
                 return lia_move::conflict;
             }
             return lia_move::undef;
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index 968333723..3fe939170 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -171,9 +171,8 @@ namespace lp {
             if (r == lia_move::conflict) {
                 de.explain(*this->m_ex);
                 return lia_move::conflict;
-            } else if (r == lia_move::sat)  {
-                return lia_move::undef;
-                NOT_IMPLEMENTED_YET();
+            } else if (r == lia_move::branch)  {
+                return lia_move::branch;
             }
 
             return lia_move::undef;
@@ -186,11 +185,12 @@ namespace lp {
         }
 
         bool should_gomory_cut() {
-            return m_number_of_calls % settings().m_int_gomory_cut_period == 0;
+            return !settings().dio_cuts()
+                && m_number_of_calls % settings().m_int_gomory_cut_period == 0;
         }
 
         bool should_solve_dioph_eq() {
-            return lia.settings().dioph_eq() && m_number_of_calls % settings().m_dioph_eq_period == 0;
+            return lia.settings().dio_eqs() && m_number_of_calls % settings().m_dioph_eq_period == 0;
         }
 
         bool should_hnf_cut() {
diff --git a/src/math/lp/lp_settings.cpp b/src/math/lp/lp_settings.cpp
index 0048f8079..7c65b9b15 100644
--- a/src/math/lp/lp_settings.cpp
+++ b/src/math/lp/lp_settings.cpp
@@ -32,5 +32,7 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
     report_frequency = p.arith_rep_freq();
     m_simplex_strategy = static_cast<lp::simplex_strategy_enum>(p.arith_simplex_strategy());
     m_nlsat_delay = p.arith_nl_delay();
-    m_dioph_eq = p.arith_lp_dioph_eq();
+    m_dio_eqs = p.arith_lp_dio_eqs();
+    m_dio_cuts = m_dio_eqs && p.arith_lp_dio_cuts();
+    m_dio_cut_from_proof_period = p.arith_lp_dio_cut_from_proof_period();
 }
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index e22a078b4..ec19012b0 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -128,6 +128,7 @@ struct statistics {
     unsigned m_grobner_conflicts = 0;
     unsigned m_offset_eqs = 0;
     unsigned m_fixed_eqs = 0;
+    unsigned m_dio_conflicts = 0;
     ::statistics m_st = {};
 
     void reset() {
@@ -217,7 +218,7 @@ public:
     unsigned         column_number_threshold_for_using_lu_in_lar_solver = 4000;
     unsigned         m_int_gomory_cut_period = 4;
     unsigned         m_int_find_cube_period = 4;
-    unsigned         m_dioph_eq_period = 4;
+    unsigned         m_dioph_eq_period = 1;
 private:
     unsigned         m_hnf_cut_period = 4;
     bool             m_int_run_gcd_test = true;
@@ -229,7 +230,10 @@ private:
     bool             m_enable_hnf = true;
     bool             m_print_external_var_name = false;
     bool             m_propagate_eqs = false;
-    bool             m_dioph_eq = false;
+    bool             m_dio_eqs = false;
+    bool             m_dio_cuts = false;
+    unsigned         m_dio_cut_from_proof_period = 3;
+
 public:
     bool print_external_var_name() const { return m_print_external_var_name; }
     bool propagate_eqs() const { return m_propagate_eqs;}
@@ -237,9 +241,11 @@ public:
     void set_hnf_cut_period(unsigned period) { m_hnf_cut_period = period;  }
     unsigned random_next() { return m_rand(); }
     unsigned random_next(unsigned u ) { return m_rand(u); }
-    bool dioph_eq() { return m_dioph_eq; }
+    bool dio_eqs() { return m_dio_eqs; }
+    bool dio_cuts() { return m_dio_eqs && m_dio_cuts; }
+    unsigned dio_cut_from_proof_period() { return m_dio_cut_from_proof_period; }
     void set_random_seed(unsigned s) { m_rand.set_seed(s); }
-
+    
     bool bound_progation() const { 
         return m_bound_propagation;
     }
diff --git a/src/smt/params/smt_params_helper.pyg b/src/smt/params/smt_params_helper.pyg
index 47a249f25..8dc6b8574 100644
--- a/src/smt/params/smt_params_helper.pyg
+++ b/src/smt/params/smt_params_helper.pyg
@@ -57,7 +57,9 @@ def_module_params(module_name='smt',
                           ('bv.solver', UINT, 0, 'bit-vector solver engine: 0 - bit-blasting, 1 - polysat, 2 - intblast, requires sat.smt=true'),
                           ('arith.random_initial_value', BOOL, False, 'use random initial values in the simplex-based procedure for linear arithmetic'),
                           ('arith.solver', UINT, 6, 'arithmetic solver: 0 - no solver, 1 - bellman-ford based solver (diff. logic only), 2 - simplex based solver, 3 - floyd-warshall based solver (diff. logic only) and no theory combination 4 - utvpi, 5 - infinitary lra, 6 - lra solver'),
-                          ('arith.lp.dioph_eq', BOOL, True, 'use Diophantine equalities'),
+                          ('arith.lp.dio_eqs', BOOL, True, 'use Diophantine equalities'),
+                          ('arith.lp.dio_cut_from_proof_period', UINT, 3, 'Period of creating a cut from proof in dioph equations'),
+                          ('arith.lp.dio_cuts', BOOL, True, 'use cuts from Diophantine equalities conflics instead of Gomory cuts, only works when dioph_eq is true'),                          
                           ('arith.nl', BOOL, True, '(incomplete) nonlinear arithmetic support based on Groebner basis and interval propagation, relevant only if smt.arith.solver=2'),
                           ('arith.nl.nra', BOOL, True, 'call nra_solver when incremental linearization does not produce a lemma, this option is ignored when arith.nl=false, relevant only if smt.arith.solver=6'),
                           ('arith.nl.branching', BOOL, True, 'branching on integer variables in non linear clusters'),