diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ
index 8e0462b7a..ffa946957 100644
--- a/docs/paper/reductions.typ
+++ b/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.
 ][
diff --git a/docs/paper/references.bib b/docs/paper/references.bib
index 0dda81fd2..d490e5dbc 100644
--- a/docs/paper/references.bib
+++ b/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}
+}
diff --git a/src/rules/mod.rs b/src/rules/mod.rs
index 7606f4923..9b9fe454e 100644
--- a/src/rules/mod.rs
+++ b/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());
diff --git a/src/rules/simultaneousincongruences_ilp.rs b/src/rules/simultaneousincongruences_ilp.rs
new file mode 100644
index 000000000..63a38df42
--- /dev/null
+++ b/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;
diff --git a/src/solvers/pipelines.rs b/src/solvers/pipelines.rs
index b9615dead..4768c5558 100644
--- a/src/solvers/pipelines.rs
+++ b/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")]),
+}
diff --git a/src/unit_tests/rules/simultaneousincongruences_ilp.rs b/src/unit_tests/rules/simultaneousincongruences_ilp.rs
new file mode 100644
index 000000000..ffe3757aa
--- /dev/null
+++ b/src/unit_tests/rules/simultaneousincongruences_ilp.rs
@@ -0,0 +1,387 @@
+use super::*;
+use crate::models::algebraic::Comparison;
+use crate::solvers::{BruteForce, ILPSolveError, ILPSolver};
+use crate::traits::Problem;
+
+type Row = (Vec<(usize, i64)>, Comparison, i64);
+
+fn reduce(
+    pairs: &[(i64, i64)],
+) -> (
+    SimultaneousIncongruences,
+    ReductionSimultaneousIncongruencesToILP,
+) {
+    let source = SimultaneousIncongruences::new(pairs.to_vec()).unwrap();
+    let reduced = ReduceTo::<ILP<i64, i64, Bounded>>::reduce_to(&source).unwrap();
+    (source, reduced)
+}
+
+fn is_feasible(target: &ILP<i64, i64, Bounded>, values: &[i64]) -> bool {
+    target.is_feasible(values).unwrap()
+}
+
+/// Whether some point of the target is feasible, deciding it by scanning `x` and, for
+/// each pair, the local `(q, r, z)` triples (the rows meet only through `x`).
+fn target_has_point(target: &ILP<i64, i64, Bounded>, m: usize) -> bool {
+    let vars = target.variables();
+    let x_hi = vars[0].upper_bound().unwrap();
+    (0..=x_hi).any(|x| {
+        (0..m).all(|i| {
+            let (q, r, z) = (1 + i, 1 + m + i, 1 + 2 * m + i);
+            let mut v = vec![0_i64; 1 + 3 * m];
+            v[0] = x;
+            (0..=vars[q].upper_bound().unwrap()).any(|qv| {
+                (0..=vars[r].upper_bound().unwrap()).any(|rv| {
+                    (0..=1).any(|zv| {
+                        v[q] = qv;
+                        v[r] = rv;
+                        v[z] = zv;
+                        target
+                            .constraints()
+                            .iter()
+                            .filter(|row| row.terms().iter().any(|&(k, _)| [q, r, z].contains(&k)))
+                            .all(|row| row.is_satisfied(&v).unwrap())
+                    })
+                })
+            })
+        })
+    })
+}
+
+fn all_pairs(max_modulus: i64) -> Vec<(i64, i64)> {
+    (1..=max_modulus)
+        .flat_map(|b| (1..=b).map(move |a| (a, b)))
+        .collect()
+}
+
+fn source_has_witness(source: &SimultaneousIncongruences) -> bool {
+    (0..source.lcm_moduli()).any(|x| source.evaluate(&x).unwrap().0)
+}
+
+#[test]
+fn issue_example_has_the_documented_target_and_witness() {
+    let (source, reduced) = reduce(&[(2, 2), (1, 3), (2, 3)]);
+    let target = reduced.target_problem();
+    assert_eq!(source.lcm_moduli(), 6);
+    assert_eq!(target.num_vars(), 10);
+    assert_eq!(target.num_constraints(), 9);
+    assert_eq!(target.num_nonzeros(), 21);
+    let bounds: Vec<(i64, i64)> = target
+        .variables()
+        .iter()
+        .map(|v| (v.lower_bound().unwrap(), v.upper_bound().unwrap()))
+        .collect();
+    assert_eq!(
+        bounds,
+        vec![
+            (0, 5),
+            (0, 2),
+            (0, 1),
+            (0, 1),
+            (0, 1),
+            (0, 2),
+            (0, 2),
+            (0, 1),
+            (0, 1),
+            (0, 1)
+        ]
+    );
+    let rows: Vec<Row> = target
+        .constraints()
+        .iter()
+        .map(|row| (row.terms().to_vec(), row.comparison(), row.rhs()))
+        .collect();
+    let expected: Vec<Row> = vec![
+        (vec![(0, 1), (1, -2), (4, -1)], Comparison::Eq, 0),
+        (vec![(4, 1), (7, -2)], Comparison::Le, -1),
+        (vec![(4, 1), (7, -2)], Comparison::Ge, -1),
+        (vec![(0, 1), (2, -3), (5, -1)], Comparison::Eq, 0),
+        (vec![(5, 1), (8, -3)], Comparison::Le, 0),
+        (vec![(5, 1), (8, -3)], Comparison::Ge, -1),
+        (vec![(0, 1), (3, -3), (6, -1)], Comparison::Eq, 0),
+        (vec![(6, 1), (9, -3)], Comparison::Le, 1),
+        (vec![(6, 1), (9, -3)], Comparison::Ge, 0),
+    ];
+    assert_eq!(rows, expected);
+    let witness = vec![3, 1, 1, 1, 1, 0, 0, 1, 0, 0];
+    assert!(is_feasible(target, &witness));
+    assert_eq!(reduced.extract_solution(&witness).unwrap(), 3);
+    assert_eq!(
+        (0..6)
+            .filter(|x| source.evaluate(x).unwrap().0)
+            .collect::<Vec<_>>(),
+        vec![3]
+    );
+    crate::rules::test_helpers::assert_parameter_predictions(&source, &reduced);
+}
+
+#[test]
+fn exhaustive_small_lists_agree_and_every_witness_extracts() {
+    let base = all_pairs(3);
+    let mut lists = vec![vec![]];
+    let mut count = 0;
+    let mut frontier = vec![vec![]];
+    for _ in 0..5 {
+        let mut next = Vec::new();
+        for list in &frontier {
+            for &p in &base {
+                let mut l: Vec<(i64, i64)> = list.clone();
+                l.push(p);
+                next.push(l);
+            }
+        }
+        lists.extend(next.iter().cloned());
+        frontier = next;
+    }
+    for pairs in &lists {
+        let (source, reduced) = reduce(pairs);
+        let target = reduced.target_problem();
+        let m = pairs.len();
+        assert_eq!(target.num_vars(), 1 + 3 * m);
+        assert_eq!(target.num_constraints(), 3 * m);
+        assert_eq!(target.num_nonzeros(), 7 * m);
+        let expected = source_has_witness(&source);
+        assert_eq!(target_has_point(target, m), expected, "{pairs:?}");
+        for x in (0..source.lcm_moduli()).filter(|x| source.evaluate(x).unwrap().0) {
+            // Lift every witness in one period and extract it back.
+            let mut v = vec![x];
+            v.extend(pairs.iter().map(|&(_, b)| x / b));
+            v.extend(pairs.iter().map(|&(_, b)| x % b));
+            v.extend(pairs.iter().map(|&(a, b)| i64::from(x % b > a % b)));
+            assert!(is_feasible(target, &v), "{pairs:?} {v:?}");
+            assert_eq!(reduced.extract_solution(&v).unwrap(), x);
+        }
+        count += 1;
+    }
+    assert_eq!(count, 1 + 6 + 36 + 216 + 1296 + 7776);
+}
+
+#[test]
+fn test_simultaneousincongruences_to_ilp_closed_loop() {
+    for pairs in [
+        vec![],
+        vec![(1, 1)],
+        vec![(2, 2)],
+        vec![(2, 2), (1, 2)],
+        vec![(2, 2), (1, 2), (1, 3)],
+        vec![(2, 2), (1, 3), (2, 3)],
+        vec![(2, 2), (1, 3), (2, 5), (3, 7)],
+        vec![(4, 4), (1, 4), (2, 4), (3, 4)],
+        vec![(1, 4), (1, 6), (5, 6)],
+    ] {
+        let (source, reduced) = reduce(&pairs);
+        crate::rules::test_helpers::assert_parameter_predictions(&source, &reduced);
+        let expected = BruteForce::new().solve(&source).unwrap();
+        match (expected, ILPSolver::new().solve(&source)) {
+            (Some(_), Ok(witness)) => assert!(source.evaluate(&witness).unwrap().0, "{pairs:?}"),
+            (None, Err(ILPSolveError::Infeasible)) => {}
+            other => panic!("{pairs:?} disagrees: {other:?}"),
+        }
+    }
+}
+
+#[test]
+fn raw_target_domains_agree_with_the_source() {
+    for pairs in [
+        vec![(1, 1)],
+        vec![(2, 2)],
+        vec![(1, 2)],
+        vec![(3, 3)],
+        vec![(1, 2), (1, 2)],
+        vec![(2, 2), (1, 2)],
+        vec![(2, 2), (1, 3)],
+        vec![(1, 1), (1, 2)],
+        vec![(2, 3), (1, 2)],
+    ] {
+        let (source, reduced) = reduce(&pairs);
+        let target = reduced.target_problem();
+        let sizes: Vec<i64> = target
+            .variables()
+            .iter()
+            .map(|v| v.upper_bound().unwrap() + 1)
+            .collect();
+        let total: i64 = sizes.iter().product();
+        let mut feasible = false;
+        for mut index in 0..total {
+            let mut v = Vec::new();
+            for &s in &sizes {
+                v.push(index % s);
+                index /= s;
+            }
+            if is_feasible(target, &v) {
+                feasible = true;
+                assert!(
+                    source
+                        .evaluate(&reduced.extract_solution(&v).unwrap())
+                        .unwrap()
+                        .0
+                );
+            }
+        }
+        assert_eq!(feasible, source_has_witness(&source), "{pairs:?}");
+    }
+}
+
+#[test]
+fn selector_lemma_for_every_modulus_up_to_two_hundred() {
+    for b in 1..=200_i64 {
+        for rho in 0..b {
+            let a = if rho == 0 { b } else { rho };
+            let (_, reduced) = reduce(&[(a, b)]);
+            let target = reduced.target_problem();
+            for r in 0..b {
+                let some_z = (0..=1).any(|z| {
+                    let v = [0, 0, r, z];
+                    target.constraints()[1..]
+                        .iter()
+                        .all(|row| row.is_satisfied(&v).unwrap())
+                });
+                assert_eq!(some_z, r != rho, "b={b} rho={rho} r={r}");
+            }
+        }
+    }
+}
+
+#[test]
+fn special_cases_empty_modulus_one_and_residue_equal_to_modulus() {
+    let (source, reduced) = reduce(&[]);
+    let target = reduced.target_problem();
+    assert_eq!(
+        (
+            target.num_vars(),
+            target.num_constraints(),
+            target.num_nonzeros()
+        ),
+        (1, 0, 0)
+    );
+    assert_eq!(target.variables()[0].upper_bound(), Some(0));
+    assert_eq!(reduced.extract_solution(&vec![0]).unwrap(), 0);
+    assert!(source.evaluate(&0).unwrap().0);
+
+    let (source, reduced) = reduce(&[(1, 1)]);
+    let target = reduced.target_problem();
+    assert_eq!(
+        (
+            target.num_vars(),
+            target.num_constraints(),
+            target.num_nonzeros()
+        ),
+        (4, 3, 7)
+    );
+    assert!(!target_has_point(target, 1));
+    assert!(!source_has_witness(&source));
+    assert_eq!(
+        ILPSolver::new().solve(&source),
+        Err(ILPSolveError::Infeasible)
+    );
+
+    // a = b gives rho = 0: the upper row forces z = 1 and the lower row r >= 1.
+    let (_, reduced) = reduce(&[(5, 5)]);
+    let target = reduced.target_problem();
+    assert!(is_feasible(target, &[1, 0, 1, 1]));
+    assert!(!is_feasible(target, &[0, 0, 0, 0]));
+    assert!(!is_feasible(target, &[0, 0, 0, 1]));
+}
+
+#[test]
+fn infeasible_instance_has_no_target_point_and_rejects_every_assignment() {
+    let (source, reduced) = reduce(&[(2, 2), (1, 2), (1, 3)]);
+    let target = reduced.target_problem();
+    assert!(!source_has_witness(&source));
+    assert!(!target_has_point(target, 3));
+    assert_eq!(
+        ILPSolver::new().solve(&source),
+        Err(ILPSolveError::Infeasible)
+    );
+    for x in 0..6_i64 {
+        let v = vec![x, x / 2, x / 2, x / 3, x % 2, x % 2, x % 3, 0, 0, 0];
+        assert!(!is_feasible(target, &v));
+        assert!(reduced.extract_solution(&v).is_err());
+    }
+}
+
+#[test]
+fn extraction_rejects_non_solutions_and_wrong_lengths() {
+    let (_, reduced) = reduce(&[(2, 2), (1, 3), (2, 3)]);
+    let good = vec![3, 1, 1, 1, 1, 0, 0, 1, 0, 0];
+    assert_eq!(reduced.extract_solution(&good).unwrap(), 3);
+    let mut wrong_selector = good.clone();
+    wrong_selector[7] = 0;
+    let mut wrong_division = good.clone();
+    wrong_division[1] = 0;
+    let mut out_of_bounds = good.clone();
+    out_of_bounds[0] = 9;
+    let mut negative = good.clone();
+    negative[4] = -1;
+    for invalid in [
+        wrong_selector,
+        wrong_division,
+        out_of_bounds,
+        negative,
+        vec![],
+        good[..9].to_vec(),
+        [good.clone(), vec![0]].concat(),
+        vec![0; 10],
+    ] {
+        assert!(reduced.extract_solution(&invalid).is_err(), "{invalid:?}");
+    }
+}
+
+#[test]
+fn largest_accepted_period_builds_exact_bounds_without_overflow() {
+    // 2^63 - 1 = 7^2 * 73 * 127 * 337 * 92737 * 649657.
+    let b1: i64 = 49 * 73 * 127 * 337;
+    let b2: i64 = 92737 * 649657;
+    assert_eq!(b1.checked_mul(b2), Some(i64::MAX));
+    let (source, reduced) = reduce(&[(3, b1), (b2, b2)]);
+    let target = reduced.target_problem();
+    assert_eq!(source.lcm_moduli(), i64::MAX);
+    assert_eq!(target.variables()[0].upper_bound(), Some(i64::MAX - 1));
+    assert_eq!(target.variables()[1].upper_bound(), Some(b2 - 1));
+    assert_eq!(target.variables()[2].upper_bound(), Some(b1 - 1));
+    assert_eq!(target.variables()[3].upper_bound(), Some(b1 - 1));
+    assert_eq!(target.variables()[4].upper_bound(), Some(b2 - 1));
+    assert_eq!(target.max_constraint_magnitude_bits(), 63);
+    let rows = target.constraints();
+    assert_eq!(rows[0].terms(), &[(0, 1), (1, -b1), (3, -1)]);
+    assert_eq!(rows[1].rhs(), 2);
+    assert_eq!(rows[2].rhs(), 3 - b1 + 1);
+    assert_eq!(rows[4].rhs(), -1);
+    assert_eq!(rows[5].rhs(), -b2 + 1);
+    // The extreme in-domain vector keeps every division-row value inside i64.
+    let top = vec![i64::MAX - 1, 0, 0, b1 - 1, b2 - 1, 0, 0];
+    assert!(rows[0].is_satisfied(&top).is_ok());
+    let low = vec![0, b2 - 1, 0, b1 - 1, 0, 0, 0];
+    assert!(rows[0].is_satisfied(&low).is_ok());
+    crate::rules::test_helpers::assert_parameter_predictions(&source, &reduced);
+
+    // One step past the boundary is not a source instance at all.
+    assert!(SimultaneousIncongruences::new(vec![(1, b1), (1, b2 + 1)]).is_err());
+
+    // Modulus 2^63 - 1 alone: the widest row coefficients.
+    let (source, reduced) = reduce(&[(i64::MAX, i64::MAX)]);
+    let target = reduced.target_problem();
+    assert_eq!(target.variables()[1].upper_bound(), Some(0));
+    assert_eq!(target.variables()[2].upper_bound(), Some(i64::MAX - 1));
+    assert_eq!(target.constraints()[2].rhs(), -(i64::MAX - 1));
+    assert_eq!(target.max_constraint_magnitude_bits(), 63);
+    assert_eq!(source.lcm_moduli(), i64::MAX);
+}
+
+#[test]
+fn magnitude_bits_stay_within_the_period_bit_length() {
+    for pairs in [
+        vec![(1, 1)],
+        vec![(2, 2)],
+        vec![(1, 4), (1, 6)],
+        vec![(1, 1 << 40), (3, 1 << 20)],
+        vec![(5, (1 << 62) - 1)],
+        vec![(1, 1 << 62)],
+    ] {
+        let (source, reduced) = reduce(&pairs);
+        assert!(
+            reduced.target_problem().max_constraint_magnitude_bits() <= source.lcm_bits(),
+            "{pairs:?}"
+        );
+    }
+}