From cd0af999a824cbcc2478639b400e2c788b9393d1 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Mon, 29 Aug 2022 14:32:13 -0700 Subject: [PATCH] fix #6302 crash due to not checking for dead rows. non-termination due to solving div and mod separately. To ensure termination one needs to at least process them simultaneously, otherwise the metric of number-of-terms x under number of mod/div does not decrease. Substituting in K*y + z under either a mod or div increases the number of terms under a mod/div when eliminating only one of the kinds. Currently handling divides constraints separately because pre-existing solution uses the model to determine z as a constant between 0 and K-1. The treatment of mod/div is supposed to be more general and use a variable while at the same time reducing the mod/div terms where the eliminated variable is used (the variable z is not added under the mod/div terms, but instead the model is used to determine cut-offs to calculate mod/div directly. --- src/math/simplex/model_based_opt.cpp | 290 ++++++++++++--------------- src/math/simplex/model_based_opt.h | 4 +- 2 files changed, 125 insertions(+), 169 deletions(-) diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp index 2e899c301..b90851582 100644 --- a/src/math/simplex/model_based_opt.cpp +++ b/src/math/simplex/model_based_opt.cpp @@ -380,40 +380,37 @@ namespace opt { m_below.reset(); for (unsigned row_id : row_ids) { SASSERT(row_id != m_objective_id); - if (visited.contains(row_id)) { + if (visited.contains(row_id)) continue; - } visited.insert(row_id); row& r = m_rows[row_id]; - if (r.m_alive) { - rational a = get_coefficient(row_id, x); - if (a.is_zero()) { - // skip - } - else if (a.is_pos() == is_pos || r.m_type == t_eq) { - rational value = x_val - (r.m_value/a); - if (bound_row_index == UINT_MAX) { - lub_val = value; - bound_row_index = row_id; - bound_coeff = a; - } - else if ((value == lub_val && r.m_type == opt::t_lt) || - (is_pos && value < lub_val) || - - (!is_pos && value > lub_val)) { - m_above.push_back(bound_row_index); - lub_val = value; - bound_row_index = row_id; - bound_coeff = a; - } - else { - m_above.push_back(row_id); - } - } - else { - m_below.push_back(row_id); - } + if (!r.m_alive) + continue; + rational a = get_coefficient(row_id, x); + if (a.is_zero()) { + // skip } + else if (a.is_pos() == is_pos || r.m_type == t_eq) { + rational value = x_val - (r.m_value/a); + if (bound_row_index == UINT_MAX) { + lub_val = value; + bound_row_index = row_id; + bound_coeff = a; + } + else if ((value == lub_val && r.m_type == opt::t_lt) || + (is_pos && value < lub_val) || + + (!is_pos && value > lub_val)) { + m_above.push_back(bound_row_index); + lub_val = value; + bound_row_index = row_id; + bound_coeff = a; + } + else + m_above.push_back(row_id); + } + else + m_below.push_back(row_id); } return bound_row_index != UINT_MAX; } @@ -1062,23 +1059,18 @@ namespace opt { unsigned eq_row = UINT_MAX; // select the lub and glb. for (unsigned row_id : row_ids) { - if (visited.contains(row_id)) { + if (visited.contains(row_id)) continue; - } visited.insert(row_id); row& r = m_rows[row_id]; - if (!r.m_alive) { + if (!r.m_alive) continue; - } rational a = get_coefficient(row_id, x); - if (a.is_zero()) { + if (a.is_zero()) continue; - } - if (r.m_type == t_eq) { + if (r.m_type == t_eq) eq_row = row_id; - continue; - } - if (r.m_type == t_mod) + else if (r.m_type == t_mod) mod_rows.push_back(row_id); else if (r.m_type == t_div) div_rows.push_back(row_id); @@ -1111,15 +1103,12 @@ namespace opt { } } - if (!mod_rows.empty()) - return solve_mod(x, mod_rows, compute_def); - - if (!div_rows.empty()) - return solve_div(x, div_rows, compute_def); - if (!divide_rows.empty()) return solve_divides(x, divide_rows, compute_def); + if (!div_rows.empty() || !mod_rows.empty()) + return solve_mod_div(x, mod_rows, div_rows, compute_def); + if (eq_row != UINT_MAX) return solve_for(eq_row, x, compute_def); @@ -1223,94 +1212,7 @@ namespace opt { // - 0 <= g*z.value + w.value < K*(g+1) // - add g*z + w - v - k*K = 0 for suitable k from 0 .. g based on model // - - model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& _mod_rows, bool compute_def) { - def result; - unsigned_vector mod_rows(_mod_rows); - rational K(1); - for (unsigned ri : mod_rows) - K = lcm(K, m_rows[ri].m_mod); - - rational x_value = m_var2value[x]; - rational y_value = div(x_value, K); - rational z_value = mod(x_value, K); - SASSERT(x_value == K * y_value + z_value); - SASSERT(0 <= z_value && z_value < K); - // add new variables - unsigned z = add_var(z_value, true); - unsigned y = add_var(y_value, true); - - uint_set visited; - unsigned j = 0; - for (unsigned ri : mod_rows) { - if (visited.contains(ri)) - continue; - m_rows[ri].m_alive = false; - visited.insert(ri); - mod_rows[j++] = ri; - } - mod_rows.shrink(j); - - // replace x by K*y + z in other rows. - for (unsigned ri : m_var2row_ids[x]) { - if (visited.contains(ri)) - continue; - replace_var(ri, x, K, y, rational::one(), z); - visited.insert(ri); - normalize(ri); - } - - // add bounds for z - add_lower_bound(z, rational::zero()); - add_upper_bound(z, K - 1); - - for (unsigned ri : mod_rows) { - rational a = get_coefficient(ri, x); - replace_var(ri, x, rational::zero()); - - // add w = b mod K - vector coeffs = m_rows[ri].m_vars; - rational coeff = m_rows[ri].m_coeff; - unsigned v = m_rows[ri].m_id; - rational v_value = m_var2value[v]; - - unsigned w = UINT_MAX; - rational offset(0); - if (coeffs.empty() || K == 1) - offset = mod(coeff, K); - else - w = add_mod(coeffs, coeff, K); - - - rational w_value = w == UINT_MAX ? offset : m_var2value[w]; - - // add v = a*z + w - V, for k = (a*z_value + w_value) div K - // claim: (= (mod x K) (- x (* K (div x K)))))) is a theorem for every x, K != 0 - rational V = v_value - a * z_value - w_value; - vector mod_coeffs; - mod_coeffs.push_back(var(v, rational::minus_one())); - mod_coeffs.push_back(var(z, a)); - if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one())); - add_constraint(mod_coeffs, V + offset, t_eq); - add_lower_bound(v, rational::zero()); - add_upper_bound(v, K - 1); - - retire_row(ri); - - project(v, false); - } - - def y_def = project(y, compute_def); - def z_def = project(z, compute_def); - - if (compute_def) { - result = (y_def * K) + z_def; - m_var2value[x] = eval(result); - } - TRACE("opt", display(tout << "solve_mod\n")); - return result; - } - + // // // Given v = a*x + b div K // Replace x |-> K*y + z @@ -1332,14 +1234,17 @@ namespace opt { // where k is between 0 and g // when gcd(a, K) = 1, then there are only two cases. // - model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& _div_rows, bool compute_def) { + model_based_opt::def model_based_opt::solve_mod_div(unsigned x, unsigned_vector const& _mod_rows, unsigned_vector const& _div_rows, bool compute_def) { def result; - unsigned_vector div_rows(_div_rows); - SASSERT(!div_rows.empty()); + unsigned_vector div_rows(_div_rows), mod_rows(_mod_rows); + SASSERT(!div_rows.empty() || !mod_rows.empty()); + TRACE("opt", display(tout << "solve_div " << x << "\n")); rational K(1); for (unsigned ri : div_rows) K = lcm(K, m_rows[ri].m_mod); + for (unsigned ri : mod_rows) + K = lcm(K, m_rows[ri].m_mod); rational x_value = m_var2value[x]; rational z_value = mod(x_value, K); @@ -1362,6 +1267,17 @@ namespace opt { div_rows[j++] = ri; } div_rows.shrink(j); + + j = 0; + for (unsigned ri : mod_rows) { + if (visited.contains(ri)) + continue; + m_rows[ri].m_alive = false; + visited.insert(ri); + mod_rows[j++] = ri; + } + mod_rows.shrink(j); + // replace x by K*y + z in other rows. for (unsigned ri : m_var2row_ids[x]) { @@ -1376,9 +1292,10 @@ namespace opt { add_lower_bound(z, rational::zero()); add_upper_bound(z, K - 1); - TRACE("opt", display(tout)); - // solve for x_value = K*y_value + z_value, 0 <= z_value < K. + // solve for x_value = K*y_value + z_value, 0 <= z_value < K. + + unsigned_vector vs; for (unsigned ri : div_rows) { @@ -1458,10 +1375,48 @@ namespace opt { add_constraint(bound_coeffs, k * K - offset, t_le); // allow to recycle row. retire_row(ri); - project(v, false); + vs.push_back(v); } - TRACE("opt", display(tout << "solve_div reduced " << y << " " << z << "\n")); + for (unsigned ri : mod_rows) { + rational a = get_coefficient(ri, x); + replace_var(ri, x, rational::zero()); + + // add w = b mod K + vector coeffs = m_rows[ri].m_vars; + rational coeff = m_rows[ri].m_coeff; + unsigned v = m_rows[ri].m_id; + rational v_value = m_var2value[v]; + + unsigned w = UINT_MAX; + rational offset(0); + if (coeffs.empty() || K == 1) + offset = mod(coeff, K); + else + w = add_mod(coeffs, coeff, K); + + + rational w_value = w == UINT_MAX ? offset : m_var2value[w]; + + // add v = a*z + w - V, for k = (a*z_value + w_value) div K + // claim: (= (mod x K) (- x (* K (div x K)))))) is a theorem for every x, K != 0 + rational V = v_value - a * z_value - w_value; + vector mod_coeffs; + mod_coeffs.push_back(var(v, rational::minus_one())); + mod_coeffs.push_back(var(z, a)); + if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one())); + add_constraint(mod_coeffs, V + offset, t_eq); + add_lower_bound(v, rational::zero()); + add_upper_bound(v, K - 1); + + retire_row(ri); + vs.push_back(v); + } + + + for (unsigned v : vs) + project(v, false); + // project internal variables. def y_def = project(y, compute_def); @@ -1525,12 +1480,12 @@ namespace opt { unsigned_vector const& row_ids = m_var2row_ids[x]; uint_set visited; for (unsigned row_id : row_ids) { - if (!visited.contains(row_id)) { - // x |-> D*y + u - replace_var(row_id, x, D, y, u); - visited.insert(row_id); - normalize(row_id); - } + if (visited.contains(row_id)) + continue; + // x |-> D*y + u + replace_var(row_id, x, D, y, u); + visited.insert(row_id); + normalize(row_id); } TRACE("opt1", display(tout << "tableau after replace x by y := v" << y << "\n");); def result = project(y, compute_def); @@ -1635,25 +1590,28 @@ namespace opt { uint_set visited; visited.insert(row_id1); for (unsigned row_id2 : row_ids) { - if (!visited.contains(row_id2)) { - visited.insert(row_id2); - b = get_coefficient(row_id2, x); - if (b.is_zero()) - continue; - row& dst = m_rows[row_id2]; - switch (dst.m_type) { - case t_eq: - case t_lt: - case t_le: - solve(row_id1, a, row_id2, x); - break; - case t_divides: - case t_mod: - case t_div: - // mod reduction already done. - UNREACHABLE(); - break; - } + if (visited.contains(row_id2)) + continue; + visited.insert(row_id2); + row& r = m_rows[row_id2]; + if (!r.m_alive) + continue; + b = get_coefficient(row_id2, x); + if (b.is_zero()) + continue; + row& dst = m_rows[row_id2]; + switch (dst.m_type) { + case t_eq: + case t_lt: + case t_le: + solve(row_id1, a, row_id2, x); + break; + case t_divides: + case t_mod: + case t_div: + // mod reduction already done. + UNREACHABLE(); + break; } } def result; diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h index 96b6508de..5ba9bd619 100644 --- a/src/math/simplex/model_based_opt.h +++ b/src/math/simplex/model_based_opt.h @@ -166,9 +166,7 @@ namespace opt { def solve_divides(unsigned x, unsigned_vector const& divide_rows, bool compute_def); - def solve_mod(unsigned x, unsigned_vector const& divide_rows, bool compute_def); - - def solve_div(unsigned x, unsigned_vector const& divide_rows, bool compute_def); + def solve_mod_div(unsigned x, unsigned_vector const& mod_rows, unsigned_vector const& divide_rows, bool compute_def); bool is_int(unsigned x) const { return m_var2is_int[x]; }