Fix #1197: [Rule] SimultaneousIncongruences to bounded integer ILP

docs/paper/reductions.typ

@@ -14867,6 +14867,87 @@ The following reductions to Integer Linear Programming are straightforward formu
   _Solution extraction._ Assign each element its position in the sorted list, rather than returning the list of elements.
 ]
 
+#let si_ilp = load-example("SimultaneousIncongruences", "ILP")
+#let si_ilp_sol = si_ilp.solutions.at(0)
+#let si_ilp_pairs = si_ilp.source.instance.pairs
+#let si_ilp_m = si_ilp_pairs.len()
+#let si_ilp_rho = si_ilp_pairs.map(p => calc.rem(p.at(0), p.at(1)))
+#let si_ilp_L = si_ilp_pairs.fold(1, (l, p) => calc.lcm(l, p.at(1)))
+#let si_ilp_x = si_ilp_sol.target_config.at(0)
+#let si_ilp_no = ((2, 2), (1, 2), (1, 3))
+#reduction-rule("SimultaneousIncongruences", "ILP",
+  example: true,
+  example-caption: [$m = #si_ilp_m$ pairs, period $L = #si_ilp_L$],
+  extra: [
+    #pred-commands(
+      "pred create --example " + rule-spec(si_ilp) + " -o si.json",
+      "pred reduce si.json --via route.json -o bundle.json",
+      "pred solve bundle.json",
+      "pred evaluate si.json --config " + cli-config(si_ilp_sol.source_config),
+    )
+
+    #{
+      assert(si_ilp.target.instance.variables.len() == 1 + 3 * si_ilp_m)
+      assert(si_ilp.target.instance.constraints.len() == 3 * si_ilp_m)
+      assert(si_ilp_sol.target_config.len() == 1 + 3 * si_ilp_m)
+      assert(si_ilp_sol.source_config == si_ilp_x)
+      for (i, p) in si_ilp_pairs.enumerate() {
+        assert(calc.rem(si_ilp_x, p.at(1)) != si_ilp_rho.at(i))
+      }
+      // The stored target vector satisfies every stored row and bound.
+      let tv = si_ilp_sol.target_config
+      for (v, b) in tv.zip(si_ilp.target.instance.variables) {
+        assert(b.lower_bound <= v and v <= b.upper_bound)
+      }
+      for c in si_ilp.target.instance.constraints {
+        let lhs = c.terms.fold(0, (acc, t) => acc + t.at(1) * tv.at(t.at(0)))
+        assert(if c.comparison == "Eq" { lhs == c.rhs } else if c.comparison == "Le" { lhs <= c.rhs } else { lhs >= c.rhs })
+      }
+      // The witness is the only one in one period.
+      assert(range(si_ilp_L).filter(x => si_ilp_pairs.all(p => calc.rem(x, p.at(1)) != calc.rem(p.at(0), p.at(1)))) == (si_ilp_x,))
+      let no_L = si_ilp_no.fold(1, (l, p) => calc.lcm(l, p.at(1)))
+      let no_rho = si_ilp_no.map(p => calc.rem(p.at(0), p.at(1)))
+      assert(no_L == 6)
+      assert(no_rho == (0, 1, 1))
+      let vb = si_ilp.target.instance.variables.map(b => $[#b.lower_bound, #b.upper_bound]$)
+      let rows = si_ilp.target.instance.constraints
+      let sym = c => if c.comparison == "Eq" { $=$ } else if c.comparison == "Le" { $<=$ } else { $>=$ }
+      let lin = c => c.terms.enumerate().map(((k, t)) => {
+        let name = if t.at(0) == 0 { $x$ } else if t.at(0) <= si_ilp_m { $q_#t.at(0)$ } else if t.at(0) <= 2 * si_ilp_m { $r_#(t.at(0) - si_ilp_m)$ } else { $z_#(t.at(0) - 2 * si_ilp_m)$ }
+        let coef = t.at(1)
+        if coef == 1 { if k == 0 { name } else { $+ name$ } } else if coef == -1 { $- name$ } else if coef < 0 { $- #(-coef) name$ } else { $#coef name$ }
+      }).join()
+      assert(range(no_L).all(x => si_ilp_no.any(p => calc.rem(x, p.at(1)) == calc.rem(p.at(0), p.at(1)))))
+      [
+        *Step 1 -- Source instance.* The pairs are #si_ilp_pairs.map(p => $(#p.at(0), #p.at(1))$).join(", "), so the forbidden residues are $rho = (#fmt-values(si_ilp_rho))$ and the period is $L = #si_ilp_L$.
+
+        *Step 2 -- Bounds and rows.* The variable intervals, in order $x, q_1, dots, q_#si_ilp_m, r_1, dots, r_#si_ilp_m, z_1, dots, z_#si_ilp_m$, are #vb.join(", "). The #rows.len() rows are:
+        #for c in rows [
+          - $#lin(c) #sym(c) #c.rhs$
+        ]
+
+        *Step 3 -- YES witness.* The stored target vector is $(#fmt-values(si_ilp_sol.target_config))$. It has $x = #si_ilp_x$; the remainders are $(#range(si_ilp_m).map(i => str(calc.rem(si_ilp_x, si_ilp_pairs.at(i).at(1)))).join(", "))$, none equal to $rho$, and extraction returns $x = #si_ilp_x$.
+
+        *Step 4 -- NO instance.* For the pairs #si_ilp_no.map(p => $(#p.at(0), #p.at(1))$).join(", ") (period #no_L, $rho = (#fmt-values(no_rho))$), pair 1 forces $r_1 = 1$ and so odd $x$, while pair 2 forces $r_2 = 0$ and so even $x$. Every $x in [0,5]$ is excluded, so the target is infeasible.
+      ]
+    }
+
+    *Multiplicity:* The fixture stores one canonical witness; here the witness in one period is unique.
+  ],
+)[
+  Quotient and remainder variables with one binary selector per pair turn each incongruence into linear rows (bounded conditional linearization @glover1975). The target has $1 + 3m$ variables, $3m$ rows and $7m$ nonzeros for $m$ pairs.
+][
+  _Construction._ Let $rho_i = a_i mod b_i$ and $L = op("lcm")(b_1, dots, b_m)$. Variables are $x in [0, L-1]$ and, for each pair, $q_i in [0, L slash b_i - 1]$, $r_i in [0, b_i - 1]$, $z_i in {0, 1}$. Rows for each pair:
+  $ x - b_i q_i - r_i = 0, quad r_i - b_i z_i <= rho_i - 1, quad r_i - b_i z_i >= (rho_i - b_i) + 1. $
+  The objective is zero.
+
+  _Correctness._ For $0 <= r <= b-1$ and $0 <= rho <= b-1$, some $z in {0,1}$ satisfies both selector rows iff $r != rho$: $z = 0$ works iff $r < rho$ (the lower row holds since $r >= 0$) and $z = 1$ works iff $r > rho$ (the upper row holds since $r <= b-1$). ($arrow.r.double$) A witness reduces modulo $L$ without changing any residue, since each $b_i$ divides $L$. Take $q_i = floor(x / b_i)$, $r_i = x mod b_i$ and $z_i = 0$ if $r_i < rho_i$, else $1$; all bounds and rows hold. ($arrow.l.double$) The division row with $0 <= r_i < b_i$ makes $r_i$ the remainder of $x$ modulo $b_i$, and the selector rows then give $r_i != rho_i$, so $x$ avoids every forbidden residue. For $b_i = 1$ neither selector value is feasible, matching the source; the empty input gives one variable fixed to $0$ and no rows.
+
+  _Arithmetic._ Every coefficient, bound and right-hand side has magnitude at most $L <= 2^63 - 1$. In the division row the running sums lie in $[0, L-1]$, then $[-(L - b_i), L-1]$, then $[-(L-1), L-1]$, so no row overflows 64-bit integers.
+
+  _Solution extraction._ Validate that the target vector is feasible, then return its first coordinate $x$.
+]
+
 #reduction-rule("TwoDimensionalConsecutiveSets", "ILP", example: true)[
   Let $P=sum_(S) binom(|S|,2)$ over the input subsets. A bounded-label formulation has $n+P$ variables, $4P$ rows and $10P$ nonzeros. Empty and singleton subsets need no rows.
 ][

docs/paper/references.bib

@@ -2352,3 +2352,14 @@ @article{een2006
   pages = {1--26},
   doi = {10.3233/SAT190014},
 }
+
+@article{glover1975,
+  author = {Fred Glover},
+  title = {Improved Linear Integer Programming Formulations of Nonlinear Integer Problems},
+  journal = {Management Science},
+  volume = {22},
+  number = {4},
+  pages = {455--460},
+  year = {1975},
+  doi = {10.1287/mnsc.22.4.455}
+}

src/rules/mod.rs

@@ -41,6 +41,7 @@ pub(crate) mod factoring_circuit;
 mod graph;
 pub(crate) mod minimumweightandorgraph_ilp;
 pub(crate) mod partitionintoperfectmatchings_ilp;
+pub(crate) mod simultaneousincongruences_ilp;
 pub(crate) mod twodimensionalconsecutivesets_ilp;
 pub(crate) use graph::{RecoveryStep, recover_completed_result};
 mod binpacking_casts;
@@ -340,6 +341,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
     specs.extend(minimumweightandorgraph_ilp::canonical_rule_example_specs());
     specs.extend(boundeddiameterspanningtree_ilp::canonical_rule_example_specs());
     specs.extend(betweenness_ilp::canonical_rule_example_specs());
+    specs.extend(simultaneousincongruences_ilp::canonical_rule_example_specs());
     specs.extend(minimumcodegenerationunlimitedregisters_ilp::canonical_rule_example_specs());
     specs.extend(cyclicordering_ilp::canonical_rule_example_specs());
     specs.extend(twodimensionalconsecutivesets_ilp::canonical_rule_example_specs());

src/rules/simultaneousincongruences_ilp.rs

@@ -0,0 +1,113 @@
+//! Reduction from SimultaneousIncongruences to a bounded-integer feasibility ILP.
+//!
+//! For every pair `(a, b)` with `rho = a mod b` the target holds a quotient `q`, a
+//! remainder `r` and a binary selector `z`. The division row `x = b q + r` makes `r`
+//! the Euclidean remainder of `x`, and the two selector rows are jointly satisfiable
+//! exactly when `r != rho`.
+
+use crate::models::algebraic::{
+    Bounded, ILP, IntegerVariable, LinearConstraint, ObjectiveSense, SimultaneousIncongruences,
+};
+use crate::rules::traits::{ReduceTo, ReductionResult};
+
+/// Result of reducing SimultaneousIncongruences to ILP.
+#[derive(Debug, Clone)]
+pub struct ReductionSimultaneousIncongruencesToILP {
+    target: ILP<i64, i64, Bounded>,
+}
+
+impl ReductionResult for ReductionSimultaneousIncongruencesToILP {
+    type Source = SimultaneousIncongruences;
+    type Target = ILP<i64, i64, Bounded>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    fn extract_solution(&self, solution: &Vec<i64>) -> crate::rules::ExtractionResult<i64> {
+        crate::rules::traits::validate_target_witness(
+            &self.target,
+            solution,
+            |value| value.value.is_some(),
+            "target ILP assignment is infeasible",
+        )?;
+        Ok(solution[0])
+    }
+}
+
+#[crate::aggregate_reduction(ilp_feasibility)]
+impl crate::rules::AggregateReductionResult for ReductionSimultaneousIncongruencesToILP {}
+
+#[crate::reduction(transform = {
+    exact {
+        num_vars = "1 + 3 * num_pairs",
+        num_constraints = "3 * num_pairs",
+        num_nonzeros = "7 * num_pairs",
+    },
+    upper_bound { max_constraint_magnitude_bits = "lcm_bits", },
+})]
+impl ReduceTo<ILP<i64, i64, Bounded>> for SimultaneousIncongruences {
+    type Result = ReductionSimultaneousIncongruencesToILP;
+
+    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
+        type Target = ILP<i64, i64, Bounded>;
+        let overflow = |operation: &str| {
+            crate::rules::ReductionError::integer_overflow::<Self, Target>(operation)
+        };
+        let pairs = self.pairs();
+        let m = pairs.len();
+        let count = m
+            .checked_mul(3)
+            .and_then(|n| n.checked_add(1))
+            .ok_or_else(|| overflow("counting target variables"))?;
+        let lcm = self.lcm_moduli();
+        let mut variables = Vec::with_capacity(count);
+        variables
+            .push(IntegerVariable::new(Some(0), Some(lcm - 1)).map_err(Self::target_construction)?);
+        for &(_, b) in pairs {
+            variables.push(
+                IntegerVariable::new(Some(0), Some(lcm / b - 1))
+                    .map_err(Self::target_construction)?,
+            );
+        }
+        for &(_, b) in pairs {
+            variables.push(
+                IntegerVariable::new(Some(0), Some(b - 1)).map_err(Self::target_construction)?,
+            );
+        }
+        variables.resize(count, IntegerVariable::binary());
+
+        let mut constraints = Vec::with_capacity(3 * m);
+        for (i, &(a, b)) in pairs.iter().enumerate() {
+            let rho = a % b;
+            let (q, r, z) = (1 + i, 1 + m + i, 1 + 2 * m + i);
+            let lower_rhs = (rho - b)
+                .checked_add(1)
+                .ok_or_else(|| overflow("computing a selector bound"))?;
+            constraints.push(LinearConstraint::eq(vec![(0, 1), (q, -b), (r, -1)], 0));
+            constraints.push(LinearConstraint::le(vec![(r, 1), (z, -b)], rho - 1));
+            constraints.push(LinearConstraint::ge(vec![(r, 1), (z, -b)], lower_rhs));
+        }
+        let target = ILP::with_variables(variables, constraints, vec![], ObjectiveSense::Minimize)
+            .map_err(Self::target_construction)?;
+        crate::rules::ilp_helpers::validate_bounded_constraint_arithmetic::<Self>(&target)?;
+        Ok(Self::Result { target })
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "simultaneousincongruences_to_ilp",
+        build: || {
+            crate::example_db::specs::rule_example_via_bounded_ilp(
+                SimultaneousIncongruences::new(vec![(2, 2), (1, 3), (2, 3)])
+                    .expect("canonical pairs are valid"),
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/simultaneousincongruences_ilp.rs"]
+mod tests;

src/solvers/pipelines.rs

@@ -996,3 +996,8 @@ register_ilp_pipeline! {
     ("DecisionILP", [("variable", "bool"), ("coefficient", "i64"), ("bounds", "general")]),
     ("ILP", [("variable", "bool"), ("coefficient", "i64"), ("bounds", "general")]),
 }
+
+register_ilp_pipeline! {
+    ("SimultaneousIncongruences", []),
+    ("ILP", [("variable", "i64"), ("coefficient", "i64"), ("bounds", "bounded")]),
+}