Lcube (#9858)

Implemented the largest cube heuristic from Bromberger and Weidenbach's
paper on cubes. Also fixes an overflow bug in mzp.
Use vswhere to find the visual studio version on windows in the build's ymls.
---------

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

.github/workflows/Windows.yml

@@ -37,9 +37,10 @@ jobs:
           ${{ matrix.cmd1 }}
           ${{ matrix.cmd2 }}
           ${{ matrix.cmd3 }}
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }}
-          cmake ${{ matrix.bindings }} -G "NMake Makefiles" ../
-          nmake    
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} || exit /b 1
+          cmake ${{ matrix.bindings }} -G "NMake Makefiles" ../ || exit /b 1
+          nmake || exit /b 1
           cd ..
         shell: cmd
       - name: Run Regressions
@@ -52,7 +53,8 @@ jobs:
         if: ${{ matrix.test }}
         run: |
            pushd build
-           call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }}
+           for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+           call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} || exit /b 1
            pushd build\python
            python z3test.py z3
            python z3test.py z3num

.github/workflows/nightly.yml

@@ -552,7 +552,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1
           python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip
       
       - name: Upload artifact
@@ -578,7 +579,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1
           python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip
       
       - name: Upload artifact
@@ -604,7 +606,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1
           python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ env.MAJOR }}.${{ env.MINOR }}.${{ env.PATCH }} --zip
       
       - name: Upload artifact

.github/workflows/nuget-build.yml

@@ -30,7 +30,8 @@ jobs:
       - name: Build Windows x64
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1
           python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip
       
       - name: Upload Windows x64 artifact
@@ -54,7 +55,8 @@ jobs:
       - name: Build Windows x86
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1
           python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip
       
       - name: Upload Windows x86 artifact
@@ -78,7 +80,8 @@ jobs:
       - name: Build Windows ARM64
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1
           python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip
       
       - name: Upload Windows ARM64 artifact

.github/workflows/release.yml

@@ -562,7 +562,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1
           python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip
       
       - name: Upload artifact
@@ -588,7 +589,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1
           python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip
       
       - name: Upload artifact
@@ -614,7 +616,8 @@ jobs:
       - name: Build
         shell: cmd
         run: |
-          call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64
+          for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i"
+          call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1
           python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ env.RELEASE_VERSION }} --zip
       
       - name: Upload artifact

src/math/lp/int_cube.cpp

@@ -16,6 +16,9 @@ Author:
 Revision History:
 --*/
 
+#include <algorithm>
+#include <unordered_map>
+
 #include "math/lp/int_solver.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/int_cube.h"
@@ -81,6 +84,252 @@ namespace lp {
         SASSERT(lp_status::OPTIMAL == lra.get_status() || lp_status::FEASIBLE == lra.get_status());
     }
 
+    // The largest cube test of Bromberger and Weidenbach:
+    //     maximize x_e subject to Ax + a'(x_e/2) <= b, x_e >= 0, where a'_i = ||a_i||_1,
+    // with the 1-norm taken over the integer variables of the row.
+    // The solution is the center z of a largest cube contained in the polyhedron.
+    // If the maximal edge length is at least 1, then the rounding of z is
+    // an integer solution; otherwise the rounding is checked, and possibly repaired,
+    // against the original constraints.
+    lia_move int_cube::find_largest_cube() {
+        lia.settings().stats().m_lcube_calls++;
+        TRACE(cube,
+              for (unsigned j = 0; j < lra.number_of_vars(); ++j)
+                  lia.display_column(tout, j);
+              tout << lra.constraints();
+              );
+
+        lra.push();
+        // The edge rows are ephemeral: suppress the add-term callback,
+        // dioph_eq's reaction to it is not undone by pop().
+        auto add_term_cb = lra.m_add_term_callback;
+        lra.m_add_term_callback = nullptr;
+        unsigned x_e = lra.add_var(UINT_MAX, false); // the edge length of the cube
+        lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(0));
+        bool ok = add_cube_edge_rows(x_e);
+        lra.m_add_term_callback = add_term_cb;
+        if (!ok) {
+            lra.pop();
+            lra.set_status(lp_status::OPTIMAL);
+            return lia_move::undef;
+        }
+
+        lp_status st = lra.find_feasible_solution();
+        if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
+            TRACE(cube, tout << "cannot find a feasible solution";);
+            lra.pop();
+            lra.move_non_basic_columns_to_bounds();
+            // it can happen that we found an integer solution here
+            return !lra.r_basis_has_inf_int()? lia_move::sat: lia_move::undef;
+        }
+
+        impq e; // the maximal edge length
+        st = lra.maximize_term(x_e, e, /*fix_int_cols*/ false);
+        if (lia.settings().get_cancel_flag()) {
+            lra.pop();
+            return lia_move::undef;
+        }
+        if (st == lp_status::UNBOUNDED) {
+            // infinite lattice width: the polyhedron contains cubes of arbitrary edge length
+            lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(1));
+            st = lra.find_feasible_solution();
+            if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
+                lra.pop();
+                return lia_move::undef;
+            }
+            lra.pop();
+            return sat_after_rounding();
+        }
+        TRACE(cube, tout << "max edge length = " << e << "\n";);
+        if (e >= impq(mpq(1))) {
+            lra.pop();
+            return sat_after_rounding();
+        }
+        // the largest cube is smaller than the unit cube:
+        // the rounded center is only a candidate
+        lra.pop();
+        return round_and_repair();
+    }
+
+    bool int_cube::add_cube_edge_rows(unsigned x_e) {
+        // snapshot the term columns: add_edge_rows_for_term appends to lra.terms()
+        svector<unsigned> term_columns;
+        for (const lar_term* t : lra.terms())
+            term_columns.push_back(t->j());
+        for (unsigned j : term_columns)
+            if (!add_edge_rows_for_term(j, x_e)) {
+                TRACE(cube, tout << "cannot add the edge rows";);
+                return false;
+            }
+        return true;
+    }
+
+    // i is the column index having the term
+    bool int_cube::add_edge_rows_for_term(unsigned i, unsigned x_e) {
+        if (!lra.column_associated_with_row(i))
+            return true;
+        const lar_term& t = lra.get_term(i);
+        impq delta = get_cube_delta_for_term(t);
+        TRACE(cube, lra.print_term_as_indices(t, tout); tout << ", delta = " << delta << "\n";);
+        if (is_zero(delta))
+            return true;
+        if (!is_zero(delta.y))
+            // the infinitesimal delta does not scale with x_e: tighten statically,
+            // it is sound for any edge length
+            return lra.tighten_term_bounds_by_delta(i, delta);
+        if (lra.column_has_upper_bound(i)) {
+            impq u = lra.get_upper_bound(i); // copy: add_term invalidates bound references
+            vector<std::pair<mpq, unsigned>> coeffs;
+            coeffs.push_back(std::make_pair(mpq(1), i));
+            coeffs.push_back(std::make_pair(delta.x, x_e));
+            unsigned s = lra.add_term(coeffs, UINT_MAX);
+            lra.add_var_bound(s, is_zero(u.y) ? lconstraint_kind::LE : lconstraint_kind::LT, u.x);
+        }
+        if (lra.column_has_lower_bound(i)) {
+            impq l = lra.get_lower_bound(i); // copy: add_term invalidates bound references
+            vector<std::pair<mpq, unsigned>> coeffs;
+            coeffs.push_back(std::make_pair(mpq(1), i));
+            coeffs.push_back(std::make_pair(-delta.x, x_e));
+            unsigned s = lra.add_term(coeffs, UINT_MAX);
+            lra.add_var_bound(s, is_zero(l.y) ? lconstraint_kind::GE : lconstraint_kind::GT, l.x);
+        }
+        return true;
+    }
+
+    lia_move int_cube::sat_after_rounding() {
+        lra.round_to_integer_solution();
+        lra.set_status(lp_status::FEASIBLE);
+        SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible());
+        TRACE(cube, tout << "largest cube success";);
+        lia.settings().stats().m_lcube_success++;
+        return lia_move::sat;
+    }
+
+    lia_move int_cube::round_and_repair() {
+        lra.backup_x(); // remember the cube center
+        vector<flip_candidate> flips;
+        for (unsigned j = 0; j < lra.column_count(); ++j) {
+            if (!lra.column_is_int(j) || lra.column_has_term(j))
+                continue;
+            const impq& v = lra.get_column_value(j);
+            if (v.is_int())
+                continue;
+            flip_candidate f;
+            f.m_j = j;
+            f.m_lo = floor(v);
+            flips.push_back(f);
+        }
+        lra.round_to_integer_solution();
+        for (auto& f : flips)
+            f.m_at_hi = lra.get_column_value(f.m_j).x > f.m_lo;
+        if (repair_rounded_candidate(flips)) {
+            lra.set_status(lp_status::FEASIBLE);
+            SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible());
+            TRACE(cube, tout << "largest cube success";);
+            lia.settings().stats().m_lcube_success++;
+            return lia_move::sat;
+        }
+        // return to the cube center: an interior point of the polyhedron
+        lra.restore_x();
+        lra.set_status(lp_status::FEASIBLE);
+        return lia_move::undef;
+    }
+
+    // Checks the rounded center against the original constraints. On failure
+    // searches the vertices of the lattice cell around the center greedily:
+    // flip a coordinate between floor and ceiling to maximally decrease the
+    // total bound violation, within a budget.
+    bool int_cube::repair_rounded_candidate(vector<flip_candidate>& flips) {
+        vector<bounded_row> rows;
+        for (const lar_term* t : lra.terms()) {
+            unsigned j = t->j();
+            if (!lra.column_associated_with_row(j))
+                continue;
+            if (!lra.column_has_upper_bound(j) && !lra.column_has_lower_bound(j))
+                continue;
+            bounded_row r;
+            r.m_j = j;
+            r.m_val = t->apply(lra.r_x());
+            rows.push_back(r);
+        }
+        auto row_violation = [&](unsigned ri, const impq& v) {
+            impq w;
+            unsigned j = rows[ri].m_j;
+            if (lra.column_has_upper_bound(j) && v > lra.get_upper_bound(j))
+                w += v - lra.get_upper_bound(j);
+            if (lra.column_has_lower_bound(j) && v < lra.get_lower_bound(j))
+                w += lra.get_lower_bound(j) - v;
+            return w;
+        };
+        impq violation;
+        for (unsigned ri = 0; ri < rows.size(); ++ri)
+            violation += row_violation(ri, rows[ri].m_val);
+        if (is_zero(violation))
+            return true; // the rounded center fits as it is
+        if (flips.empty())
+            return false;
+
+        std::unordered_map<unsigned, unsigned> flip_of_var;
+        for (unsigned fi = 0; fi < flips.size(); ++fi)
+            flip_of_var[flips[fi].m_j] = fi;
+        // occurrences of the flip candidates in the bounded rows
+        vector<vector<std::pair<unsigned, mpq>>> occs;
+        occs.resize(flips.size());
+        for (unsigned ri = 0; ri < rows.size(); ++ri) {
+            const lar_term& t = lra.get_term(rows[ri].m_j);
+            for (lar_term::ival p : t) {
+                auto it = flip_of_var.find(p.j());
+                if (it != flip_of_var.end())
+                    occs[it->second].push_back(std::make_pair(ri, p.coeff()));
+            }
+        }
+
+        unsigned budget = std::min(2 * flips.size(), lia.settings().lcube_flips());
+        bool flipped = false;
+        while (!is_zero(violation) && budget-- > 0) {
+            unsigned best_fi = UINT_MAX;
+            impq best_gain;
+            for (unsigned fi = 0; fi < flips.size(); ++fi) {
+                if (occs[fi].empty())
+                    continue;
+                mpq step = flips[fi].m_at_hi ? mpq(-1) : mpq(1);
+                impq gain;
+                for (const auto& o : occs[fi]) {
+                    const impq& v = rows[o.first].m_val;
+                    gain += row_violation(o.first, v + impq(step * o.second)) - row_violation(o.first, v);
+                }
+                if (gain < best_gain) {
+                    best_gain = gain;
+                    best_fi = fi;
+                }
+            }
+            if (best_fi == UINT_MAX)
+                return false; // no flip decreases the violation
+            mpq step = flips[best_fi].m_at_hi ? mpq(-1) : mpq(1);
+            for (const auto& o : occs[best_fi])
+                rows[o.first].m_val += impq(step * o.second);
+            flips[best_fi].m_at_hi = !flips[best_fi].m_at_hi;
+            violation += best_gain;
+            flipped = true;
+            TRACE(cube, tout << "flipped column " << flips[best_fi].m_j << ", violation = " << violation << "\n";);
+        }
+        if (!is_zero(violation))
+            return false;
+
+        // apply the repaired candidate
+        for (const auto& f : flips)
+            lra.set_column_value(f.m_j, impq(f.m_at_hi ? f.m_lo + 1 : f.m_lo));
+        for (const lar_term* t : lra.terms()) {
+            unsigned j = t->j();
+            if (!lra.column_associated_with_row(j))
+                continue;
+            lra.set_column_value(j, t->apply(lra.r_x()));
+        }
+        if (flipped)
+            lia.settings().stats().m_lcube_flip_success++;
+        return true;
+    }
+
     impq int_cube::get_cube_delta_for_term(const lar_term& t) const {
         if (t.size() == 2) {
             bool seen_minus = false;

src/math/lp/int_cube.h

@@ -10,9 +10,15 @@ Abstract:
     Cube finder
 
     This routine attempts to find a feasible integer solution
-    by tightnening bounds and running an LRA solver on the 
+    by tightnening bounds and running an LRA solver on the
     tighter system.
 
+    find_largest_cube() implements the largest cube test of
+    Bromberger and Weidenbach (Fast Cube Tests for LIA Constraint
+    Solving, IJCAR 2016): a fresh variable x_e for the cube edge
+    length is introduced and maximized; the center of the largest
+    cube is rounded to a candidate integer solution.
+
 Author:
     Nikolaj Bjorner (nbjorner)
     Lev Nachmanson (levnach)
@@ -21,7 +27,10 @@ Revision History:
 --*/
 #pragma once
 
+#include "util/vector.h"
 #include "math/lp/lia_move.h"
+#include "math/lp/numeric_pair.h"
+#include "math/lp/lar_term.h"
 
 namespace lp {
     class int_solver;
@@ -29,12 +38,30 @@ namespace lp {
     class int_cube {
         class int_solver& lia;
         class lar_solver& lra;
+        // a fractional integer coordinate of the cube center:
+        // the candidate value is m_lo or m_lo + 1
+        struct flip_candidate {
+            unsigned m_j = 0;
+            mpq      m_lo;
+            bool     m_at_hi = false;
+        };
+        // a term column with at least one bound, tracked during the repair
+        struct bounded_row {
+            unsigned m_j = 0;
+            impq     m_val;
+        };
         bool tighten_term_for_cube(unsigned i);
         bool tighten_terms_for_cube();
         void find_feasible_solution();
         impq get_cube_delta_for_term(const lar_term& t) const;
+        bool add_edge_rows_for_term(unsigned i, unsigned x_e);
+        bool add_cube_edge_rows(unsigned x_e);
+        lia_move sat_after_rounding();
+        lia_move round_and_repair();
+        bool repair_rounded_candidate(vector<flip_candidate>& flips);
     public:
         int_cube(int_solver& lia);
         lia_move operator()();
+        lia_move find_largest_cube();
     };
 }

src/math/lp/int_solver.cpp

@@ -44,7 +44,8 @@ namespace lp {
         dioph_eq            m_dio;  
         int_gcd_test        m_gcd;
         unsigned            m_initial_dio_calls_period;
-        
+        unsigned            m_lcube_period;
+
         bool column_is_int_inf(unsigned j) const {
             return lra.column_is_int(j) && (!lia.value_is_int(j));
         }
@@ -52,7 +53,8 @@ namespace lp {
         imp(int_solver& lia): lia(lia), lra(lia.lra), lrac(lia.lrac), m_hnf_cutter(lia), m_dio(lia), m_gcd(lia) {
             m_hnf_cut_period = settings().hnf_cut_period();
             m_initial_dio_calls_period = settings().dio_calls_period();
-        } 
+            m_lcube_period = settings().m_int_find_cube_period;
+        }
 
         bool has_lower(unsigned j) const {
             switch (lrac.m_column_types()[j]) {
@@ -196,6 +198,25 @@ namespace lp {
             return m_number_of_calls % settings().m_int_find_cube_period == 0;
         }
 
+        // The largest cube test is throttled exponentially: when the polyhedron
+        // does not contain a large enough cube it is unlikely to contain one
+        // later, after more constraints are added, so each failure doubles the
+        // period and a success resets it.
+        bool should_find_lcube() {
+            return settings().lcube() && m_number_of_calls % m_lcube_period == 0;
+        }
+
+        lia_move find_lcube() {
+            lia_move r = int_cube(lia).find_largest_cube();
+            if (r == lia_move::undef) {
+                if (m_lcube_period < (1u << 30))
+                    m_lcube_period *= 2;
+            }
+            else
+                m_lcube_period = settings().m_int_find_cube_period;
+            return r;
+        }
+
         bool should_gomory_cut() {
             bool dio_allows_gomory = !settings().dio() || settings().dio_enable_gomory_cuts() ||
                                       m_dio.some_terms_are_ignored();
@@ -246,6 +267,7 @@ namespace lp {
             ++m_number_of_calls;
             if (r == lia_move::undef) r = patch_basic_columns();
             if (r == lia_move::undef && should_find_cube()) r = int_cube(lia)();
+            if (r == lia_move::undef && should_find_lcube()) r = find_lcube();
             if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds();
             if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut();
             if (r == lia_move::undef && should_gomory_cut()) r = gomory(lia).get_gomory_cuts(2);

src/math/lp/lar_solver.cpp

@@ -864,7 +864,7 @@ namespace lp {
     }
 
     lp_status lar_solver::maximize_term(unsigned j,
-        impq& term_max) {
+        impq& term_max, bool fix_int_cols) {
         TRACE(lar_solver, print_values(tout););
         SASSERT(get_core_solver().m_r_solver.calc_current_x_is_feasible_include_non_basis());
         lar_term term = get_term_to_maximize(j);
@@ -879,6 +879,11 @@ namespace lp {
             return lp_status::UNBOUNDED;
         }
 
+        if (!fix_int_cols) {
+            set_status(lp_status::OPTIMAL);
+            return lp_status::OPTIMAL;
+        }
+
         impq opt_val = term_max;
 
         bool change = false;

src/math/lp/lar_solver.h

@@ -205,7 +205,9 @@ public:
         set_column_value(j, v);
     }
 
-    lp_status maximize_term(unsigned j_or_term, impq& term_max);
+    // fix_int_cols: after maximizing try to move the integer columns to integer values;
+    // pass false to keep the optimal (possibly fractional) vertex intact, e.g., for the largest cube test
+    lp_status maximize_term(unsigned j_or_term, impq& term_max, bool fix_int_cols);
 
     core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const;
     

src/math/lp/lp_params_helper.pyg

@@ -8,6 +8,8 @@ def_module_params(module_name='lp',
                           ('dio_cuts_enable_hnf', BOOL, True, 'enable hnf cuts together with Diophantine cuts, only relevant when dioph_eq is true'),
                           ('dio_ignore_big_nums', BOOL, True, 'Ignore the terms with big numbers in the Diophantine handler, only relevant when dioph_eq is true'),
                           ('dio_calls_period', UINT, 1, 'Period of calling the Diophantine handler in the final_check()'),
-                          ('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'),                          
+                          ('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'),
+                          ('lcube', BOOL, True, 'use the largest cube test for integer feasibility'),
+                          ('lcube_flips', UINT, 16, 'maximal number of coordinate flips when repairing the rounded largest cube center, only relevant when lcube is true'),
                          ))
                          

src/math/lp/lp_settings.cpp

@@ -43,5 +43,7 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
     m_dio_ignore_big_nums = lp_p.dio_ignore_big_nums();
     m_dio_calls_period = lp_p.dio_calls_period();
     m_dio_run_gcd = lp_p.dio_run_gcd();
+    m_lcube = lp_p.lcube();
+    m_lcube_flips = lp_p.lcube_flips();
     m_max_conflicts = p.max_conflicts();
 }

src/math/lp/lp_settings.h

@@ -112,6 +112,9 @@ struct statistics {
     unsigned m_gcd_conflicts = 0;
     unsigned m_cube_calls = 0;
     unsigned m_cube_success = 0;
+    unsigned m_lcube_calls = 0;
+    unsigned m_lcube_success = 0;
+    unsigned m_lcube_flip_success = 0;
     unsigned m_patches = 0;
     unsigned m_patches_success = 0;
     unsigned m_hnf_cutter_calls = 0;
@@ -152,6 +155,9 @@ struct statistics {
         st.update("arith-gcd-conflict", m_gcd_conflicts);
         st.update("arith-cube-calls", m_cube_calls);
         st.update("arith-cube-success", m_cube_success);
+        st.update("arith-lcube-calls", m_lcube_calls);
+        st.update("arith-lcube-success", m_lcube_success);
+        st.update("arith-lcube-flip-success", m_lcube_flip_success);
         st.update("arith-patches", m_patches);
         st.update("arith-patches-success", m_patches_success);
         st.update("arith-hnf-calls", m_hnf_cutter_calls);
@@ -258,7 +264,11 @@ private:
     bool             m_dio_ignore_big_nums = false;
     unsigned         m_dio_calls_period = 4;
     bool             m_dio_run_gcd = true;
+    bool             m_lcube = true;
+    unsigned         m_lcube_flips = 16;
 public:
+    bool lcube() const { return m_lcube; }
+    unsigned lcube_flips() const { return m_lcube_flips; }
     unsigned dio_calls_period() const { return m_dio_calls_period; }
     unsigned & dio_calls_period() { return m_dio_calls_period; }
     bool print_external_var_name() const { return m_print_external_var_name; }

src/smt/theory_lra.cpp

@@ -4029,7 +4029,7 @@ public:
         if (!lp().is_feasible() || lp().has_changed_columns())
             make_feasible();
         vi = get_lpvar(v);
-        auto st = lp().maximize_term(vi, term_max);
+        auto st = lp().maximize_term(vi, term_max, /*fix_int_cols*/ true);
         if (has_int() && lp().has_inf_int()) {
             st = lp::lp_status::FEASIBLE;
             lp().restore_x();

src/test/CMakeLists.txt

@@ -79,6 +79,7 @@ add_executable(test-z3
   "${CMAKE_CURRENT_BINARY_DIR}/install_tactic.cpp"
   interval.cpp
   karr.cpp
+  lcube.cpp
   list.cpp
   main.cpp
   map.cpp

src/test/lcube.cpp

@@ -0,0 +1,261 @@
+/*++
+  Copyright (c) 2020 Microsoft Corporation
+
+  Module Name:
+
+  lcube.cpp
+
+  Abstract:
+
+  Tests for the largest cube test of Bromberger and Weidenbach
+  (Fast Cube Tests for LIA Constraint Solving, IJCAR 2016),
+  implemented in int_cube::find_largest_cube().
+
+  This file lives directly under src/test (not src/test/lp) so that the
+  scripts/mk_make.py build, which only compiles the top-level test
+  directory, links tst_lcube().
+
+  Author:
+
+  Lev Nachmanson (levnach)
+
+  --*/
+
+#include <initializer_list>
+#include <iostream>
+#include <utility>
+
+#include "util/debug.h"
+#include "util/params.h"
+#include "math/lp/int_cube.h"
+#include "math/lp/int_solver.h"
+#include "math/lp/lar_solver.h"
+#include "math/lp/numeric_pair.h"
+
+namespace lp {
+
+// tests for the largest cube test of Bromberger and Weidenbach
+namespace lcube_test {
+
+    struct ineq {
+        vector<std::pair<mpq, unsigned>> m_coeffs;
+        lconstraint_kind m_kind;
+        mpq m_rs;
+    };
+
+    // builds for every inequality a term with the bound and solves
+    static void setup(lar_solver& solver, const vector<ineq>& ineqs, svector<unsigned>* term_columns = nullptr) {
+        unsigned term_ext = 1000;
+        for (const auto& in : ineqs) {
+            unsigned t = solver.add_term(in.m_coeffs, term_ext++);
+            solver.add_var_bound(t, in.m_kind, in.m_rs);
+            if (term_columns)
+                term_columns->push_back(t);
+        }
+        auto st = solver.solve();
+        VERIFY(st == lp_status::OPTIMAL || st == lp_status::FEASIBLE);
+    }
+
+    static void verify_model(const lar_solver& solver, const vector<ineq>& ineqs) {
+        for (const auto& in : ineqs) {
+            impq v;
+            for (const auto& p : in.m_coeffs)
+                v += solver.get_column_value(p.second) * p.first;
+            switch (in.m_kind) {
+            case lconstraint_kind::LE: VERIFY(v <= impq(in.m_rs)); break;
+            case lconstraint_kind::GE: VERIFY(v >= impq(in.m_rs)); break;
+            default: VERIFY(false);
+            }
+        }
+    }
+
+    static void verify_int_values(const lar_solver& solver, std::initializer_list<unsigned> vars) {
+        for (unsigned j : vars)
+            VERIFY(solver.get_column_value(j).is_int());
+    }
+
+    // The example of Bromberger and Weidenbach: 3x1 - x2 <= 0, -2x1 - x2 <= -2, -2x1 + x2 <= 1.
+    // The largest cube is smaller than the unit cube, the rounded center is not a solution,
+    // no coordinate flip repairs it (the only integer solution lies outside the lattice cell
+    // of the center): expect undef and an intact solver state.
+    static void test_paper_example_undef() {
+        std::cout << "lcube: paper example, expecting undef\n";
+        lar_solver solver;
+        unsigned x1 = solver.add_named_var(0, true, "x1");
+        unsigned x2 = solver.add_named_var(1, true, "x2");
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(3), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(0)});
+        ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(-2)});
+        ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(1), x2}}, lconstraint_kind::LE, mpq(1)});
+        setup(solver, ineqs);
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        lia_move m = int_cube(i_s).find_largest_cube();
+        std::cout << "lcube returned " << lia_move_to_string(m) << "\n";
+        VERIFY(m == lia_move::undef);
+        VERIFY(solver.ax_is_correct());
+    }
+
+    // 3x1 - x2 <= 0, -2x1 - x2 <= -1, -2x1 + x2 <= 1: the largest cube has
+    // edge 4/17 with the center (3/17, 1) that rounds to the solution (0, 1),
+    // while the unit cube test fails: the largest cube test is stronger here.
+    static void test_beats_unit_cube() {
+        std::cout << "lcube: beating the unit cube test\n";
+        lar_solver solver;
+        unsigned x1 = solver.add_named_var(0, true, "x1");
+        unsigned x2 = solver.add_named_var(1, true, "x2");
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(3), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(0)});
+        ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(-1)});
+        ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(1), x2}}, lconstraint_kind::LE, mpq(1)});
+        svector<unsigned> tcols;
+        setup(solver, ineqs, &tcols);
+        // move the solution to a fractional feasible point: the cube tests only
+        // run when the current solution has fractional integer variables
+        solver.set_column_value_test(x1, impq(mpq(1, 4)));
+        solver.set_column_value_test(x2, impq(mpq(5, 4)));
+        solver.set_column_value_test(tcols[0], impq(mpq(-1, 2)));
+        solver.set_column_value_test(tcols[1], impq(mpq(-7, 4)));
+        solver.set_column_value_test(tcols[2], impq(mpq(3, 4)));
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        lia_move m = int_cube(i_s)();
+        std::cout << "unit cube returned " << lia_move_to_string(m) << "\n";
+        VERIFY(m == lia_move::undef);
+        m = int_cube(i_s).find_largest_cube();
+        std::cout << "lcube returned " << lia_move_to_string(m) << "\n";
+        VERIFY(m == lia_move::sat);
+        verify_int_values(solver, {x1, x2});
+        verify_model(solver, ineqs);
+    }
+
+    // 9/10 <= x + y + r <= 11/10, -11/10 <= x - y + r <= 1/10, 0 <= r <= 1/10,
+    // x, y integer, r real. The real variable keeps the terms and their bounds
+    // non-integer. A fractional center, e.g. (1/2, 1/2), rounds to an infeasible
+    // point that is repaired by flipping one coordinate: expect sat.
+    static void test_flip_repair() {
+        std::cout << "lcube: rounding repair\n";
+        lar_solver solver;
+        unsigned x = solver.add_named_var(0, true, "x");
+        unsigned y = solver.add_named_var(1, true, "y");
+        unsigned r = solver.add_named_var(2, false, "r");
+        solver.add_var_bound(r, lconstraint_kind::GE, mpq(0));
+        solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10));
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(9, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(11, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(-11, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(1, 10)});
+        setup(solver, ineqs);
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        lia_move m = int_cube(i_s).find_largest_cube();
+        std::cout << "lcube returned " << lia_move_to_string(m)
+                  << ", flip successes: " << solver.settings().stats().m_lcube_flip_success << "\n";
+        VERIFY(m == lia_move::sat);
+        verify_int_values(solver, {x, y});
+        verify_model(solver, ineqs);
+    }
+
+    // 3x + 5y >= 7 alone has infinite lattice width: the edge length is
+    // unbounded and any cube center of edge >= 1 rounds to a solution.
+    static void test_infinite_lattice_width() {
+        std::cout << "lcube: infinite lattice width\n";
+        lar_solver solver;
+        unsigned x = solver.add_named_var(0, true, "x");
+        unsigned y = solver.add_named_var(1, true, "y");
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(3), x}, {mpq(5), y}}, lconstraint_kind::GE, mpq(7)});
+        setup(solver, ineqs);
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        lia_move m = int_cube(i_s).find_largest_cube();
+        std::cout << "lcube returned " << lia_move_to_string(m) << "\n";
+        VERIFY(m == lia_move::sat);
+        verify_int_values(solver, {x, y});
+        verify_model(solver, ineqs);
+    }
+
+    // 0 <= x + 2y + r <= 8, 0 <= x - 2y + r <= 8: the maximal edge length is
+    // 8/3 >= 1, so the rounded center is guaranteed to be a solution.
+    static void test_edge_at_least_one() {
+        std::cout << "lcube: edge length at least 1\n";
+        lar_solver solver;
+        unsigned x = solver.add_named_var(0, true, "x");
+        unsigned y = solver.add_named_var(1, true, "y");
+        unsigned r = solver.add_named_var(2, false, "r");
+        solver.add_var_bound(r, lconstraint_kind::GE, mpq(0));
+        solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10));
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(2), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(0)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(2), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(8)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-2), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(0)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-2), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(8)});
+        setup(solver, ineqs);
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        lia_move m = int_cube(i_s).find_largest_cube();
+        std::cout << "lcube returned " << lia_move_to_string(m) << "\n";
+        VERIFY(m == lia_move::sat);
+        verify_int_values(solver, {x, y});
+        verify_model(solver, ineqs);
+    }
+
+    // runs the flip-repair instance through int_solver::check() with the
+    // lp.lcube parameter set and the cube period lowered to 1: verifies the
+    // dispatch and the parameter plumbing
+    static void test_dispatch() {
+        std::cout << "lcube: dispatch through int_solver::check\n";
+        lar_solver solver;
+        params_ref p;
+        p.set_bool("lcube", true);
+        solver.settings().updt_params(p);
+        VERIFY(solver.settings().lcube());
+        solver.settings().m_int_find_cube_period = 1;
+        unsigned x = solver.add_named_var(0, true, "x");
+        unsigned y = solver.add_named_var(1, true, "y");
+        unsigned r = solver.add_named_var(2, false, "r");
+        solver.add_var_bound(r, lconstraint_kind::GE, mpq(0));
+        solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10));
+        vector<ineq> ineqs;
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(9, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(11, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(-11, 10)});
+        ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(1, 10)});
+        svector<unsigned> tcols;
+        setup(solver, ineqs, &tcols);
+        // a fractional feasible point, so that check() does not return sat right away
+        solver.set_column_value_test(x, impq(mpq(1, 2)));
+        solver.set_column_value_test(y, impq(mpq(1, 2)));
+        solver.set_column_value_test(r, impq(mpq(1, 10)));
+        solver.set_column_value_test(tcols[0], impq(mpq(11, 10)));
+        solver.set_column_value_test(tcols[1], impq(mpq(11, 10)));
+        solver.set_column_value_test(tcols[2], impq(mpq(1, 10)));
+        solver.set_column_value_test(tcols[3], impq(mpq(1, 10)));
+        int_solver i_s(solver);
+        solver.set_int_solver(&i_s);
+        explanation ex;
+        lia_move m = i_s.check(&ex);
+        std::cout << "check returned " << lia_move_to_string(m)
+                  << ", lcube calls: " << solver.settings().stats().m_lcube_calls << "\n";
+        VERIFY(m == lia_move::sat);
+        VERIFY(solver.settings().stats().m_lcube_calls >= 1);
+        verify_int_values(solver, {x, y});
+        verify_model(solver, ineqs);
+    }
+
+    static void run() {
+        test_paper_example_undef();
+        test_beats_unit_cube();
+        test_flip_repair();
+        test_infinite_lattice_width();
+        test_edge_at_least_one();
+        test_dispatch();
+        std::cout << "lcube tests passed\n";
+    }
+} // namespace lcube_test
+}  // namespace lp
+
+void tst_lcube() {
+    lp::lcube_test::run();
+}

src/test/lp/lp.cpp

@@ -1645,7 +1645,7 @@ void test_maximize_term() {
     lia_move lm = i_solver.check(&ex);
     VERIFY(lm == lia_move::sat);
     impq term_max;
-    lp_status st = solver.maximize_term(term_2x_pl_2y, term_max);
+    lp_status st = solver.maximize_term(term_2x_pl_2y, term_max, /*fix_int_cols*/ true);
 
     std::cout << "status = " << lp_status_to_string(st) << std::endl;
     std::cout << "term_max = " << term_max << std::endl;

src/test/main.cpp

@@ -193,7 +193,8 @@
     X(ho_matcher) \
     X(finite_set) \
     X(finite_set_rewriter) \
-    X(fpa)
+    X(fpa) \
+    X(lcube)
 
 #define FOR_EACH_TEST(X, X_ARGV) \
     FOR_EACH_ALL_TEST(X, X_ARGV) \

src/util/mpz.cpp

@@ -1904,8 +1904,11 @@ std::string mpz_manager<SYNCH>::to_string(mpz const & a) const {
 
 template<bool SYNCH>
 unsigned mpz_manager<SYNCH>::hash(mpz const & a) {
-    if (is_small(a))
-        return ::abs(a.m_val);
+    if (is_small(a)) {
+        // compute abs in unsigned arithmetic: ::abs(INT_MIN) is undefined
+        unsigned u = static_cast<unsigned>(a.m_val);
+        return a.m_val < 0 ? 0u - u : u;
+    }
 #ifndef _MP_GMP
     unsigned sz = size(a);
     if (sz == 1)