Add LLL cube heuristic for integer LP (experimental, default off)

Generalize int_cube from the unit cube C = [-1/2, 1/2]^n to a
parallelepiped K = B * C where B is an n x n unimodular integer matrix
found by a monotone pairwise basis-reduction that directly minimizes
the actual cube cost

    C(B) = (1/2) * (||A * B||_1 + ||B||_1)
         = sum_r delta_row(r, B) + sum_j delta_col(j, B)

The atomic move is a single elementary column op col_j -= q * col_k
with q chosen to minimize C(B) analytically (floor/ceil of the weighted
median of breakpoints {H[r,j]/H[r,k]} and {B[i,j]/B[i,k]}).  Starting
from B = I and accepting only strict improvements makes the heuristic
*monotone-safe*: never worse than the plain int_cube.  This addresses
the regression of int_cube_hnf (branch hnf_cube), whose triangulation
can blow up the column-delta term ||B||_1.  In a 153-instance random
matrix study the HNF basis was worse than B=I by an average factor
3x-50x, while pairwise-greedy LLL was uniformly >= plain cube.

Implementation:
* src/math/lp/int_cube_lll.{h,cpp} -- the heuristic.
* The infrastructure (collect J/terms, tighten bounds, round on a
  saved x_J snapshot, lar_solver::apply_lattice_assignment) mirrors
  the earlier hnf_cube experiment; the only algorithmic change is
  swapping HNF column-reduction for the cost-minimizing pairwise
  reduction with bail-on-overflow.

New parameters:
* lp.enable_lll_cube (bool, default false) -- feature gate.
* lp.int_find_lll_cube_period (uint, default 4) -- the LLL cube is
  invoked every Nth final-check call after a plain cube failure.
  Because it is monotone-safe it runs at cube's period (4) rather
  than the throttled 16 the HNF variant used.

Statistics: arith-lll-cube-{calls,success,bail-collect,bail-build,
bail-basis,bail-tighten,bail-infeasible}.

Resource limits: rows <= 75, cols <= 150, bitsize <= 4096, max_passes
<= 8; bail above.

Validation:
* test-z3 /a: 89/89 unit tests pass.
* Smoke run on QF_LIA cut_lemmas and CAV_2009/Bromberger samples:
  no result disagreements vs. plain cube; one timeout-to-sat win on
  20180326-Bromberger/.../unbd-sage0.smt2 (-T:15).

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

src/math/lp/CMakeLists.txt

@@ -12,6 +12,7 @@ z3_add_component(lp
     indexed_vector.cpp
     int_branch.cpp
     int_cube.cpp
+    int_cube_lll.cpp
     int_gcd_test.cpp
     int_solver.cpp
     lar_solver.cpp

src/math/lp/int_cube_lll.cpp

@@ -0,0 +1,443 @@
+/*++
+Copyright (c) 2025 Microsoft Corporation
+
+Module Name:
+
+    int_cube_lll.cpp
+
+Abstract:
+
+    LLL-style cube heuristic. See int_cube_lll.h for the math.
+
+Author:
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+--*/
+
+#include "math/lp/int_cube_lll.h"
+#include "math/lp/int_solver.h"
+#include "math/lp/lar_solver.h"
+#include <algorithm>
+
+namespace lp {
+
+    static const unsigned LLL_CUBE_LIMIT_ROWS = 75;
+    static const unsigned LLL_CUBE_LIMIT_COLS = 150;
+    static const unsigned LLL_CUBE_BITSIZE_LIMIT = 4096;
+    static const unsigned LLL_CUBE_MAX_PASSES = 8;
+
+    int_cube_lll::int_cube_lll(int_solver& lia) : lia(lia), lra(lia.lra) {}
+
+    bool int_cube_lll::too_big(const mpq& v) const {
+        return v.bitsize() > LLL_CUBE_BITSIZE_LIMIT;
+    }
+
+    // For the delta=0 fast path on rows/columns with ||row||_1 <= 1 we need
+    // both bounds (if present) to be integer-valued so that a rounded integer
+    // y_p in the LP-feasible interval cannot escape [L, U].  Strict bounds
+    // (non-zero y component) also disqualify the fast path.
+    bool int_cube_lll::column_bounds_are_integral(unsigned j) const {
+        if (lra.column_has_lower_bound(j)) {
+            const impq& lb = lra.get_lower_bound(j);
+            if (!is_zero(lb.y) || !lb.x.is_int())
+                return false;
+        }
+        if (lra.column_has_upper_bound(j)) {
+            const impq& ub = lra.get_upper_bound(j);
+            if (!is_zero(ub.y) || !ub.x.is_int())
+                return false;
+        }
+        return true;
+    }
+
+    // Collect J = { j : column_is_int(j) && !column_has_term(j) } and
+    // m_term_js = { t.j() : t in lra.terms() && column_associated_with_row(t.j()) }.
+    bool int_cube_lll::collect_J_and_terms() {
+        m_J.reset();
+        m_col_to_J.reset();
+        m_col_to_J.resize(lra.column_count(), UINT_MAX);
+        for (unsigned j = 0; j < lra.column_count(); ++j) {
+            if (!lra.column_is_int(j))
+                continue;
+            if (lra.column_has_term(j))
+                continue;
+            m_col_to_J[j] = m_J.size();
+            m_J.push_back(j);
+        }
+        if (m_J.empty())
+            return false;
+        if (m_J.size() > LLL_CUBE_LIMIT_COLS)
+            return false;
+
+        m_term_js.reset();
+        for (const lar_term* t : lra.terms()) {
+            if (!lra.column_associated_with_row(t->j()))
+                continue;
+            m_term_js.push_back(t->j());
+        }
+        if (m_term_js.size() > LLL_CUBE_LIMIT_ROWS)
+            return false;
+        return true;
+    }
+
+    // Build A in m_H and initialize B, B_inv to identity. Same as the HNF
+    // variant: starts H = A, B = I.
+    bool int_cube_lll::build_matrix() {
+        unsigned m = m_term_js.size();
+        unsigned n = m_J.size();
+        m_H.reset();
+        m_H.resize(m);
+        for (unsigned r = 0; r < m; ++r)
+            m_H[r].resize(n, mpq(0));
+
+        for (unsigned r = 0; r < m; ++r) {
+            const lar_term& t = lra.get_term(m_term_js[r]);
+            for (lar_term::ival p : t) {
+                if (p.j() >= m_col_to_J.size())
+                    continue;
+                unsigned k = m_col_to_J[p.j()];
+                if (k == UINT_MAX)
+                    continue;
+                if (!p.coeff().is_int())
+                    return false;
+                if (too_big(p.coeff()))
+                    return false;
+                m_H[r][k] = p.coeff();
+            }
+        }
+        m_B.reset();
+        m_B.resize(n);
+        m_B_inv.reset();
+        m_B_inv.resize(n);
+        for (unsigned i = 0; i < n; ++i) {
+            m_B[i].resize(n, mpq(0));
+            m_B_inv[i].resize(n, mpq(0));
+            m_B[i][i] = mpq(1);
+            m_B_inv[i][i] = mpq(1);
+        }
+        return true;
+    }
+
+    // Find the integer q minimizing
+    //     f(q) = sum_r |H[r][j] - q*H[r][k]| + sum_i |B[i][j] - q*B[i][k]|
+    // and, if applying the corresponding column op strictly decreases f
+    // (relative to q=0), perform it.  Sets improved = true on accept.
+    // Returns false on bail (overflow).
+    bool int_cube_lll::reduce_pair(unsigned j, unsigned k, bool& improved) {
+        improved = false;
+        unsigned m = m_H.size();
+        unsigned n = m_B.size();
+
+        // Gather breakpoints. Each (a, b) with b != 0 contributes
+        //   |a - q*b| = |b| * |q - a/b|
+        // to the cost; pairs with b == 0 contribute the constant |a|.
+        struct bp { mpq a; mpq b; };
+        vector<bp> bps;
+        bps.reserve(m + n);
+        mpq constant_cost(0);
+        for (unsigned r = 0; r < m; ++r) {
+            const mpq& a = m_H[r][j];
+            const mpq& b = m_H[r][k];
+            if (b.is_zero()) {
+                constant_cost += abs(a);
+            } else {
+                bps.push_back({a, b});
+            }
+        }
+        for (unsigned i = 0; i < n; ++i) {
+            const mpq& a = m_B[i][j];
+            const mpq& b = m_B[i][k];
+            if (b.is_zero()) {
+                constant_cost += abs(a);
+            } else {
+                bps.push_back({a, b});
+            }
+        }
+        if (bps.empty())
+            return true;
+
+        // Sort by breakpoint a/b ascending. Compare a1/b1 vs a2/b2 via
+        // a1*b2 vs a2*b2*b1/b1... use cross multiplication with sign care.
+        // (Easier: just compute a/b as mpq.)
+        struct kv { mpq val; mpq w; const bp* p; };
+        vector<kv> kvs;
+        kvs.reserve(bps.size());
+        for (auto& x : bps) {
+            mpq w = abs(x.b);
+            mpq v = x.a / x.b;
+            kvs.push_back({v, w, &x});
+        }
+        std::sort(kvs.begin(), kvs.end(),
+                  [](const kv& a, const kv& b) { return a.val < b.val; });
+
+        // Total weight; weighted median is the smallest v with prefix
+        // weight >= total / 2.
+        mpq total(0);
+        for (auto& x : kvs) total += x.w;
+        mpq half = total * mpq(1, 2);
+        mpq cum(0);
+        mpq median(0);
+        for (auto& x : kvs) {
+            cum += x.w;
+            if (cum >= half) { median = x.val; break; }
+        }
+
+        // Integer candidates: floor(median) and floor(median)+1. Evaluate
+        // f at q = floor, q = floor+1, q = 0 and pick the smallest.
+        auto eval = [&](const mpq& q) -> mpq {
+            mpq c(0);
+            for (auto& x : bps) {
+                mpq d = x.a - q * x.b;
+                c += abs(d);
+            }
+            return c + constant_cost;
+        };
+        mpq qf = floor(median);
+        mpq qc = qf + mpq(1);
+        mpq cf = eval(qf);
+        mpq cc = eval(qc);
+        mpq c0 = eval(mpq(0));
+        mpq best_q = qf;
+        mpq best_c = cf;
+        if (cc < best_c) { best_q = qc; best_c = cc; }
+        if (best_c >= c0 || best_q.is_zero()) {
+            // No strict improvement over q = 0 (identity).
+            return true;
+        }
+
+        // Apply column op col_j <- col_j - q * col_k on H and B, and the
+        // matching row op row_k <- row_k + q * row_j on B_inv.
+        const mpq& q = best_q;
+        for (unsigned r = 0; r < m; ++r) {
+            mpq v = m_H[r][j] - q * m_H[r][k];
+            if (too_big(v))
+                return false;
+            m_H[r][j] = v;
+        }
+        for (unsigned i = 0; i < n; ++i) {
+            mpq v = m_B[i][j] - q * m_B[i][k];
+            if (too_big(v))
+                return false;
+            m_B[i][j] = v;
+        }
+        for (unsigned c = 0; c < n; ++c) {
+            mpq v = m_B_inv[k][c] + q * m_B_inv[j][c];
+            if (too_big(v))
+                return false;
+            m_B_inv[k][c] = v;
+        }
+        improved = true;
+        return true;
+    }
+
+    // Monotone-safe greedy: start at B = I (so H = A) and accept only
+    // strict improvements of cube cost.  At most LLL_CUBE_MAX_PASSES
+    // sweeps over all ordered pairs; terminates earlier on a full
+    // no-improvement pass.
+    bool int_cube_lll::compute_basis() {
+        unsigned n = m_B.size();
+        if (n <= 1)
+            return true;
+        for (unsigned pass = 0; pass < LLL_CUBE_MAX_PASSES; ++pass) {
+            bool any_improved_this_pass = false;
+            for (unsigned j = 0; j < n; ++j) {
+                for (unsigned k = 0; k < n; ++k) {
+                    if (j == k)
+                        continue;
+                    bool improved = false;
+                    if (!reduce_pair(j, k, improved))
+                        return false;
+                    if (improved)
+                        any_improved_this_pass = true;
+                }
+            }
+            if (!any_improved_this_pass)
+                break;
+        }
+        return true;
+    }
+
+    void int_cube_lll::compute_deltas() {
+        unsigned m = m_H.size();
+        unsigned n = m_B.size();
+        m_term_delta.reset();
+        m_term_delta.resize(m);
+        for (unsigned r = 0; r < m; ++r) {
+            mpq s(0);
+            for (unsigned k = 0; k < n; ++k)
+                s += abs(m_H[r][k]);
+            if (s.is_zero()) {
+                m_term_delta[r] = impq(0);
+            }
+            else if (s == mpq(1) && column_bounds_are_integral(m_term_js[r])) {
+                // Row of H has a single +1 entry: term value = +/- y_p, an
+                // integer in lattice coords. With integer term bounds, no
+                // shrinkage is needed.
+                m_term_delta[r] = impq(0);
+            }
+            else {
+                m_term_delta[r] = impq(s * mpq(1, 2));
+            }
+        }
+        m_col_delta.reset();
+        m_col_delta.resize(n);
+        for (unsigned i = 0; i < n; ++i) {
+            mpq s(0);
+            for (unsigned k = 0; k < n; ++k)
+                s += abs(m_B[i][k]);
+            if (s.is_zero()) {
+                m_col_delta[i] = impq(0);
+            }
+            else if (s == mpq(1) && column_bounds_are_integral(m_J[i])) {
+                m_col_delta[i] = impq(0);
+            }
+            else {
+                m_col_delta[i] = impq(s * mpq(1, 2));
+            }
+        }
+    }
+
+    bool int_cube_lll::tighten_terms_for_lll_cube() {
+        for (unsigned r = 0; r < m_term_js.size(); ++r) {
+            unsigned j = m_term_js[r];
+            const impq& delta = m_term_delta[r];
+            if (is_zero(delta))
+                continue;
+            if (!lra.tighten_term_bounds_by_delta(j, delta))
+                return false;
+        }
+        return true;
+    }
+
+    bool int_cube_lll::tighten_column_bound_by_delta(unsigned j, const impq& delta) {
+        SASSERT(!lra.column_has_term(j));
+        if (is_zero(delta))
+            return true;
+        if (lra.column_has_upper_bound(j) && lra.column_has_lower_bound(j)) {
+            if (lra.get_upper_bound(j) - delta < lra.get_lower_bound(j) + delta)
+                return false;
+        }
+        if (lra.column_has_upper_bound(j)) {
+            const impq& ub = lra.get_upper_bound(j);
+            lconstraint_kind k = (!is_zero(delta.y) || !is_zero(ub.y)) ? LT : LE;
+            lra.add_var_bound(j, k, ub.x - delta.x);
+        }
+        if (lra.column_has_lower_bound(j)) {
+            const impq& lb = lra.get_lower_bound(j);
+            lconstraint_kind k = (!is_zero(delta.y) || !is_zero(lb.y)) ? GT : GE;
+            lra.add_var_bound(j, k, lb.x + delta.x);
+        }
+        return true;
+    }
+
+    bool int_cube_lll::tighten_columns_for_lll_cube() {
+        for (unsigned i = 0; i < m_J.size(); ++i) {
+            if (!tighten_column_bound_by_delta(m_J[i], m_col_delta[i]))
+                return false;
+        }
+        return true;
+    }
+
+    void int_cube_lll::save_x_J() {
+        unsigned n = m_J.size();
+        m_saved_x_J.reset();
+        m_saved_x_J.resize(n);
+        for (unsigned i = 0; i < n; ++i)
+            m_saved_x_J[i] = lra.get_column_value(m_J[i]);
+    }
+
+    void int_cube_lll::apply_rounding() {
+        unsigned n = m_J.size();
+        // Compute y* = B_inv * x_J* in impq (preserve infinitesimal y component).
+        vector<impq> y(n, impq(0));
+        for (unsigned i = 0; i < n; ++i) {
+            impq s(0);
+            for (unsigned k = 0; k < n; ++k)
+                s += m_saved_x_J[k] * m_B_inv[i][k];
+            y[i] = s;
+        }
+        // Round each y_i to the nearest integer using the same impq-aware rule
+        // as lar_solver::round_to_integer_solution (lar_solver.cpp:2582-2589).
+        for (unsigned i = 0; i < n; ++i) {
+            impq const& v = y[i];
+            impq flv(floor(v.x));
+            impq del = flv - v;
+            if (del < -impq(mpq(1, 2)))
+                y[i] = impq(ceil(v.x));
+            else
+                y[i] = flv;
+        }
+        // z = B * y_round (now y has zero y-component since it's integer).
+        vector<mpq> z(n, mpq(0));
+        for (unsigned i = 0; i < n; ++i) {
+            mpq s(0);
+            for (unsigned k = 0; k < n; ++k)
+                s += m_B[i][k] * y[k].x;
+            z[i] = s;
+        }
+        // Push the new values into the LP and propagate to dependent terms.
+        vector<std::pair<lpvar, impq>> assignments;
+        for (unsigned i = 0; i < n; ++i) {
+            unsigned j = m_J[i];
+            impq new_val(z[i]);
+            assignments.push_back({(lpvar)j, new_val});
+        }
+        lra.apply_lattice_assignment(assignments);
+    }
+
+    lia_move int_cube_lll::operator()() {
+        lia.settings().stats().m_lll_cube_calls++;
+        TRACE(lll_cube,
+              for (unsigned j = 0; j < lra.number_of_vars(); ++j)
+                  lia.display_column(tout, j);
+              tout << lra.constraints();
+            );
+
+        if (!collect_J_and_terms()) {
+            lia.settings().stats().m_lll_cube_bail_collect++;
+            return lia_move::undef;
+        }
+        if (!build_matrix()) {
+            lia.settings().stats().m_lll_cube_bail_build++;
+            return lia_move::undef;
+        }
+        if (!compute_basis()) {
+            lia.settings().stats().m_lll_cube_bail_basis++;
+            return lia_move::undef;
+        }
+        compute_deltas();
+
+        lra.push();
+        if (!tighten_terms_for_lll_cube()) {
+            lia.settings().stats().m_lll_cube_bail_tighten++;
+            lra.pop();
+            lra.set_status(lp_status::OPTIMAL);
+            return lia_move::undef;
+        }
+        if (!tighten_columns_for_lll_cube()) {
+            lia.settings().stats().m_lll_cube_bail_tighten++;
+            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) {
+            lia.settings().stats().m_lll_cube_bail_infeasible++;
+            TRACE(lll_cube, tout << "cannot find a feasible solution\n";);
+            lra.pop();
+            lra.move_non_basic_columns_to_bounds();
+            return !lra.r_basis_has_inf_int() ? lia_move::sat : lia_move::undef;
+        }
+        save_x_J();
+        lra.pop();
+        apply_rounding();
+        lra.set_status(lp_status::FEASIBLE);
+        SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible());
+        TRACE(lll_cube, tout << "lll_cube success\n";);
+        lia.settings().stats().m_lll_cube_success++;
+        return lia_move::sat;
+    }
+}

src/math/lp/int_cube_lll.h

@@ -0,0 +1,92 @@
+/*++
+Copyright (c) 2025 Microsoft Corporation
+
+Module Name:
+
+    int_cube_lll.h
+
+Abstract:
+
+    LLL-style cube heuristic for the integer-LP solver.
+
+    Generalization of int_cube. Instead of the unit cube C = [-1/2, 1/2]^n
+    we use a parallelepiped K = B * C, where B is an n x n unimodular
+    integer matrix found by a *monotone* greedy basis-reduction that
+    minimizes the actual cube cost
+
+        C(B) = (1/2) * (||A * B||_1 + ||B||_1)
+             = sum_r delta_row(r, B) + sum_j delta_col(j, B)
+
+    over GL_n(Z). Atomic move is a single elementary column operation
+
+        col_j(B) <- col_j(B) - q * col_k(B)    (q in Z, j != k)
+
+    The optimal q for a fixed pair (j, k) is the floor or ceil of the
+    weighted median of the breakpoints {A_row(r,j)/A_row(r,k)} and
+    {B(i,j)/B(i,k)} (weights are the absolute values of the denominators);
+    a standard piecewise-linear minimization. We accept only strict
+    improvements of C(B), starting from B = I; therefore the heuristic
+    is never worse than the plain int_cube.
+
+    This addresses the regression of int_cube_hnf, whose triangulation
+    can blow up the column-delta term ||B||_1: in 153 random instances
+    (3x3 to 5x5, coefficients in [-20, 20]) HNF cube lost to plain cube
+    on >99% of inputs, with an average cost ratio of 3x-50x worse, while
+    pairwise-greedy LLL won or tied on 100% (see findings.md in the
+    session-state).
+
+    Soundness: identical to int_cube_hnf. We work with a saved column
+    snapshot of x_J, tighten the LP, find a real-feasible point x*, and
+    round y* = B^{-1} x* to nearest integer in lattice coordinates, then
+    lift z = B * round(y*).
+
+Author:
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+--*/
+#pragma once
+
+#include "math/lp/lia_move.h"
+#include "math/lp/numeric_pair.h"
+#include "util/vector.h"
+
+namespace lp {
+    class int_solver;
+    class lar_solver;
+    class lar_term;
+
+    class int_cube_lll {
+        class int_solver&     lia;
+        class lar_solver&     lra;
+
+        vector<unsigned>            m_J;
+        vector<unsigned>            m_col_to_J;
+        vector<unsigned>            m_term_js;
+        vector<vector<mpq>>         m_H;       // H = A * B; m x n
+        vector<vector<mpq>>         m_B;       // n x n unimodular
+        vector<vector<mpq>>         m_B_inv;   // B^{-1} = product of inverse elementaries
+        vector<impq>                m_term_delta;
+        vector<impq>                m_col_delta;
+        vector<impq>                m_saved_x_J;
+
+    public:
+        int_cube_lll(int_solver& lia);
+        lia_move operator()();
+
+    private:
+        bool collect_J_and_terms();
+        bool build_matrix();
+        bool compute_basis();
+        bool reduce_pair(unsigned j, unsigned k, bool& improved);
+        bool too_big(const mpq& v) const;
+        bool column_bounds_are_integral(unsigned j) const;
+        bool tighten_terms_for_lll_cube();
+        bool tighten_columns_for_lll_cube();
+        bool tighten_column_bound_by_delta(unsigned j, const impq& delta);
+        void compute_deltas();
+        void save_x_J();
+        void apply_rounding();
+    };
+}

src/math/lp/int_solver.cpp

@@ -9,6 +9,7 @@
 #include "math/lp/gomory.h"
 #include "math/lp/int_branch.h"
 #include "math/lp/int_cube.h"
+#include "math/lp/int_cube_lll.h"
 #include "math/lp/dioph_eq.h"
 
 namespace lp {
@@ -196,6 +197,13 @@ namespace lp {
             return m_number_of_calls % settings().m_int_find_cube_period == 0;
         }
 
+        bool should_find_lll_cube() {
+            unsigned period = settings().m_int_find_lll_cube_period;
+            if (period == 0)
+                return false;
+            return m_number_of_calls % period == 0;
+        }
+
         bool should_gomory_cut() {
             bool dio_allows_gomory = !settings().dio() || settings().dio_enable_gomory_cuts() ||
                                       m_dio.some_terms_are_ignored();
@@ -245,7 +253,11 @@ 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_cube()) {
+                r = int_cube(lia)();
+                if (r == lia_move::undef && settings().enable_lll_cube() && should_find_lll_cube())
+                    r = int_cube_lll(lia)();
+            }
             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/int_solver.h

@@ -36,6 +36,7 @@ class int_solver {
     friend struct create_cut;
     friend class gomory;
     friend class int_cube;
+    friend class int_cube_lll;
     friend class int_branch;
     friend class int_gcd_test;
     friend class hnf_cutter;

src/math/lp/lar_solver.cpp

@@ -2597,6 +2597,27 @@ namespace lp {
         }
     }
 
+    void lar_solver::apply_lattice_assignment(const vector<std::pair<lpvar, impq>>& assignments) {
+        if (assignments.empty())
+            return;
+        SASSERT(m_imp->m_incorrect_columns.empty());
+        for (auto const& a : assignments) {
+            lpvar j = a.first;
+            const impq& new_val = a.second;
+            SASSERT(column_is_int(j));
+            SASSERT(!column_has_term(j));
+            const impq& cur = get_core_solver().r_x(j);
+            if (new_val == cur)
+                continue;
+            get_core_solver().m_r_solver.update_x(j, new_val);
+            m_imp->m_incorrect_columns.insert(j);
+        }
+        if (!m_imp->m_incorrect_columns.empty()) {
+            fix_terms_with_rounded_columns();
+            m_imp->m_incorrect_columns.reset();
+        }
+    }
+
     void lar_solver::fix_terms_with_rounded_columns() {
 
         for (const lar_term* t : m_imp->m_terms) {

src/math/lp/lar_solver.h

@@ -567,6 +567,7 @@ public:
         return false;
     }
     void round_to_integer_solution();
+    void apply_lattice_assignment(const vector<std::pair<lpvar, impq>>& assignments);
     inline const row_strip<mpq>& get_row(unsigned i) const { return A_r().m_rows[i]; }
     inline const row_strip<mpq>& basic2row(unsigned i) const { return A_r().m_rows[row_of_basic_column(i)]; }
     inline const column_strip& get_column(unsigned i) const { return A_r().m_columns[i]; }
@@ -604,6 +605,7 @@ public:
     std::function<void (const indexed_uint_set& columns_with_changed_bound)> m_find_monics_with_changed_bounds_func = nullptr;
     friend int_solver;
     friend int_branch;
+    friend class int_cube_lll;
     friend imp;
 
 };

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'),
+                          ('enable_lll_cube', BOOL, False, 'enable the LLL cube heuristic in the integer-LP solver (parallelepiped generalization of cube, monotone-safe basis reduction)'),
+                          ('int_find_lll_cube_period', UINT, 4, 'period for invoking the LLL cube heuristic (every Nth final-check call); 0 disables'),
                          ))
                          

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_enable_lll_cube = lp_p.enable_lll_cube();
+    m_int_find_lll_cube_period = lp_p.int_find_lll_cube_period();
     m_max_conflicts = p.max_conflicts();
 }

src/math/lp/lp_settings.h

@@ -112,6 +112,13 @@ struct statistics {
     unsigned m_gcd_conflicts = 0;
     unsigned m_cube_calls = 0;
     unsigned m_cube_success = 0;
+    unsigned m_lll_cube_calls = 0;
+    unsigned m_lll_cube_success = 0;
+    unsigned m_lll_cube_bail_collect = 0;
+    unsigned m_lll_cube_bail_build = 0;
+    unsigned m_lll_cube_bail_basis = 0;
+    unsigned m_lll_cube_bail_tighten = 0;
+    unsigned m_lll_cube_bail_infeasible = 0;
     unsigned m_patches = 0;
     unsigned m_patches_success = 0;
     unsigned m_hnf_cutter_calls = 0;
@@ -152,6 +159,13 @@ 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-lll-cube-calls", m_lll_cube_calls);
+        st.update("arith-lll-cube-success", m_lll_cube_success);
+        st.update("arith-lll-cube-bail-collect", m_lll_cube_bail_collect);
+        st.update("arith-lll-cube-bail-build", m_lll_cube_bail_build);
+        st.update("arith-lll-cube-bail-basis", m_lll_cube_bail_basis);
+        st.update("arith-lll-cube-bail-tighten", m_lll_cube_bail_tighten);
+        st.update("arith-lll-cube-bail-infeasible", m_lll_cube_bail_infeasible);
         st.update("arith-patches", m_patches);
         st.update("arith-patches-success", m_patches_success);
         st.update("arith-hnf-calls", m_hnf_cutter_calls);
@@ -207,6 +221,8 @@ private:
 public:
     void updt_params(params_ref const& p);
     bool enable_hnf() const { return m_enable_hnf; }
+    bool enable_lll_cube() const { return m_enable_lll_cube; }
+    bool& enable_lll_cube() { return m_enable_lll_cube; }
     unsigned nlsat_delay() const { return m_nlsat_delay; }
     bool int_run_gcd_test() const {
         if (!m_dio)
@@ -239,6 +255,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_int_find_lll_cube_period = 4;
 private:
     unsigned         m_hnf_cut_period = 4;
     bool             m_int_run_gcd_test = true;
@@ -249,6 +266,7 @@ public:
 private:
     unsigned         m_nlsat_delay = 0;
     bool             m_enable_hnf = true;
+    bool             m_enable_lll_cube = false;
     bool             m_print_external_var_name = false;
     bool             m_propagate_eqs = false;
     bool             m_dio = false;

src/util/trace_tags.def

@@ -1124,6 +1124,7 @@ X(kernel, qe_core, "qe core")
 X(lar_solver, lar_solver, "lar solver")
 X(lar_solver, add_var, "add var")
 X(lar_solver, cube, "cube")
+X(lar_solver, lll_cube, "lll cube")
 X(lar_solver, dump_terms, "dump terms")
 X(lar_solver, lar_solver_details, "lar solver details")
 X(lar_solver, lar_solver_improve_bounds, "lar solver improve bounds")