diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ
index 8e0462b7a..f6e675f18 100644
--- a/docs/paper/reductions.typ
+++ b/docs/paper/reductions.typ
@@ -12343,6 +12343,156 @@ _Solution extraction and bounds._ Copy the subset bits; aggregate recovery is id
  The canonical instance and target coefficients are shown below. Its source witness is #repr(ex.solutions.at(0).source_config); the displayed target witness decodes to it and meets the objective or feasibility certificate described above.
 ]
 
+// ComparativeContainment → DecisionQUBO
+#let cc-qubo = load-example("ComparativeContainment", "DecisionQUBO",
+  source-variant: (weight: "One"), target-variant: (weight: "i64"))
+#let cc-source = cc-qubo.source.instance
+#let cc-target = cc-qubo.target.instance
+#let cc-sol = cc-qubo.solutions.at(0)
+#let cc-n = cc-source.universe_size
+#let cc-r = cc-source.r_sets.len()
+#let cc-s = cc-source.s_sets.len()
+#let cc-m = cc-r + 1
+#let cc-outside(a) = range(cc-n).filter(v => not a.contains(v))
+#let cc-g = (cc-source.r_sets + cc-source.s_sets).map(a => calc.max(0, cc-outside(a).len()-1)).sum(default: 0)
+#let cc-selected(bits) = range(cc-n).filter(v => bits.at(v))
+#let cc-count(family, bits) = family.filter(a => cc-selected(bits).all(v => a.contains(v))).len()
+#let cc-energy(bits) = cc-target.inner.entries.map(e => if bits.at(e.at(0)) and bits.at(e.at(1)) { e.at(2) } else { 0 }).sum(default: 0)
+#let cc-extend(bits) = {
+  let aux = ()
+  for a in cc-source.r_sets + cc-source.s_sets {
+    let outside = cc-outside(a)
+    if outside.len() >= 2 {
+      for k in range(2, outside.len()+1) {
+        aux.push(outside.slice(0,k).all(v => not bits.at(v)))
+      }
+    }
+  }
+  bits + aux
+}
+#assert(cc-target.inner.num_vars == cc-n + cc-g)
+#assert(cc-extend(cc-sol.source_config) == cc-sol.target_config)
+#assert(cc-source.s_sets.all(a => range(cc-n).all(v => a.contains(v))))
+#assert(cc-source.r_sets.all(a => cc-outside(a).len() == 2))
+#assert(cc-target.bound == -(cc-m * cc-g + cc-s))
+#assert(cc-count(cc-source.r_sets,cc-sol.source_config) >= cc-count(cc-source.s_sets,cc-sol.source_config))
+#assert(cc-energy(cc-sol.target_config) <= cc-target.bound)
+#let cc-faulty = cc-sol.target_config.enumerate().map(pair => if pair.at(0) == cc-n { not pair.at(1) } else { pair.at(1) })
+#let cc-candidate = cc-sol.source_config.enumerate().map(pair => pair.at(1) or pair.at(0) == cc-source.r_sets.at(0).at(0))
+#assert(cc-energy(cc-faulty) > cc-target.bound)
+#assert(cc-count(cc-source.r_sets,cc-candidate) < cc-count(cc-source.s_sets,cc-candidate))
+#assert(cc-energy(cc-extend(cc-candidate)) > cc-target.bound)
+#let cc-no-s = cc-source.s_sets + cc-source.s_sets
+#let cc-no-bound = cc-target.bound - cc-source.s_sets.len()
+#assert(cc-no-s.len() > cc-r)
+#assert(cc-energy(cc-sol.target_config) > cc-no-bound)
+#reduction-rule("ComparativeContainment", "DecisionQUBO",
+  example: true,
+  example-source-variant: (weight: "One"),
+  example-target-variant: (weight: "i64"),
+  extra: [
+    #pred-commands(
+      "pred create --example " + rule-spec(cc-qubo) + " -o containment.json",
+      "pred path ComparativeContainment/One DecisionQUBO/i64 -o paths.json",
+      "jq '.paths[] | select(.steps == 1)' paths.json > route.json",
+      "pred reduce containment.json --via route.json -o bundle.json",
+      "pred solve bundle.json --solver brute-force",
+      "pred evaluate containment.json --config " + cli-config(cc-sol.source_config),
+    )
+  ],
+)[
+  Unit Comparative Containment reduces to integer decision QUBO by chaining
+  containment indicators. The indicator viewpoint follows the
+  #link("https://github.com/CodingThrust/problem-reductions/issues/401")[antecedent containment formulation];
+  the quadratic product penalty is Rosenberg's, as presented in @borosgruber2014quadratization.
+][
+  _Construction._ Let $U={0,dots,n-1}$ be the universe, and $R,S$ the indexed
+  families of sizes $r,s$. Put the universe selection bits $x_v$ first. For
+  each entry $A$, its containment indicator is $t_A=product_(v in U without A)(1-x_v)$.
+  Enumerate its complement increasingly. An empty complement uses constant 1;
+  one element uses its complemented bit; otherwise use a fresh chain of
+  $|U without A|-1$ product bits. Process R before S in input order and never
+  share chain bits across entries. For each binary gate $z=a b$, add
+  $P(a,b,z)=a b-2a z-2b z+3z$. With $M=r+1$, form
+  $H=M sum_"gates" P + sum_(A in S)t_A-sum_(A in R)t_A=C+Q$.
+  Expand into unique upper-triangular entries (diagonals are linear, each
+  coefficient used once), discard zero coefficients, and use bound $-C$.
+
+  _Correctness._ ($arrow.r.double$) A feasible source subset extends by its
+  true prefix products. Every penalty vanishes and $H$ is its S count minus
+  its R count, at most zero; hence $Q<=-C$. ($arrow.l.double$) On any target
+  assignment, each indicator remains binary, so the signed indicator sum is
+  at least $-r$. The gate penalty is nonnegative, vanishes exactly at $z=a b$,
+  and is at least 1 otherwise. One faulty gate would give $H>=M-r=1$.
+  A target witness therefore has every gate correct. Induction along each
+  chain recovers actual containment, and $H<=0$ gives the source inequality.
+
+  _Solution extraction._ Validate the target predicate once and return its
+  first $n$ bits. Reject wrong lengths and assignments above the bound.
+  Completed `Or` values map identically, including NO without a witness.
+  An empty family entry contains only the empty selected subset; a full-universe
+  entry always contributes 1 and needs no gates. With empty universe the
+  single empty witness works exactly when $r>=s$. Empty families and repeated
+  entries retain their count semantics; repeated elements inside a set do
+  not change membership. In particular, empty R with one full-universe S entry
+  is NO: $Q=0$ and bound $-1$.
+
+  _Size and arithmetic._ Write $g=sum_A max(0,n-|A|-1)$ over all indexed
+  entries. There are exactly $n+g$ variables. A first gate has six nonconstant
+  entries and three off-diagonal terms; a later gate has five and three.
+  Each indicator adds at most one diagonal. Thus there are at most $6g+r+s$
+  stored entries and $3g$ off-diagonal terms; coefficient collection can
+  cancel entries. Since $g<=n(r+s)$, the registered bounds are $n+n(r+s)$
+  variables and $6n(r+s)+r+s$ off-diagonal terms. The latter is conservative;
+  neither registered formula asserts equality. Complement scans cost
+  $O(n(r+s))$ plus reading the raw listings; sparse collection adds a
+  logarithmic factor in the number of gate and indicator records.
+
+  Checked size arithmetic and conversions precede expansion. Require
+  $M=r+1<=2^63-1$ and $B=8M g+r+s<=2^63-1$, computing $M g$ before multiplying
+  by 8 so zero gates remain valid. Expanded gates have nonconstant absolute
+  coefficient sums at most $8M$ and constants at most $M$. Consequently the
+  nonconstant norm is at most $B$, and separately $|C|<=M g+r+s<=B$.
+  Every coefficient, accumulation prefix and evaluation on every assignment
+  fits `i64`, even for infeasible assignments. All operations remain checked;
+  an unrepresentable target returns a typed error.
+
+  *YES example.* The loaded source has universe size #cc-n, R entries
+  #repr(cc-source.r_sets) and S entries #repr(cc-source.s_sets). Here
+  $r=#cc-r$, $s=#cc-s$, $M=#cc-m$, $g=#cc-g$: the #cc-n universe bits and
+  #cc-g chain bits give #cc-target.inner.num_vars variables. The target has
+  #cc-target.inner.entries.len() stored entries, of which
+  #(cc-target.inner.entries.filter(e => e.at(0) != e.at(1)).len()) are
+  off-diagonal. Its constant is #(-cc-target.bound), bound #cc-target.bound,
+  coefficient norm #(cc-target.inner.entries.map(e => calc.abs(e.at(2))).sum(default: 0)),
+  and arithmetic certificate #(8*cc-m*cc-g+cc-r+cc-s).
+  Its R complement pairs are #repr(cc-source.r_sets.map(cc-outside)).
+  The full polynomial is specified by triples #repr(cc-target.inner.entries).
+  The selected subset #repr(cc-selected(cc-sol.source_config)) has R count
+  #cc-count(cc-source.r_sets,cc-sol.source_config) and S count
+  #cc-count(cc-source.s_sets,cc-sol.source_config). Its extension
+  #repr(cc-sol.target_config) has $Q=#cc-energy(cc-sol.target_config)$ and
+  $H=#(cc-energy(cc-sol.target_config)-cc-target.bound)$, and extracts back
+  to #repr(cc-sol.source_config). Flipping its first chain bit gives
+  $Q=#cc-energy(cc-faulty)$ and $H=#(cc-energy(cc-faulty)-cc-target.bound)$,
+  which fails the bound. The subset #repr(cc-selected(cc-candidate)) has
+  counts #cc-count(cc-source.r_sets,cc-candidate) and
+  #cc-count(cc-source.s_sets,cc-candidate); its true extension has
+  $Q=#cc-energy(cc-extend(cc-candidate))$, also rejected.
+
+  *NO example.* Double the loaded S family, keeping the same universe and R.
+  Its #cc-no-s.len() full-universe entries always contribute that count,
+  while at most #cc-r R entries can contain any selected subset. The count
+  difference is therefore at least #(cc-no-s.len()-cc-r), strictly positive.
+  There are still #cc-g gates and #cc-target.inner.num_vars variables, and
+  the same nonconstant Q. The constant becomes #(-cc-no-bound), bound
+  #cc-no-bound, and arithmetic certificate #(8*cc-m*cc-g+cc-r+cc-no-s.len()).
+  Every consistent target misses the bound by at least
+  #(cc-no-s.len()-cc-r); every inconsistent target has $H>=M-r=1$.
+  The previous vector still gives $Q=#cc-energy(cc-sol.target_config)$,
+  missing the new bound by #(cc-energy(cc-sol.target_config)-cc-no-bound).
+]
+
 #reduction-rule("ThreeDimensionalMatching", "DecisionQUBO", example: true,
  extra: [#direct-qubo-commands("ThreeDimensionalMatching", "DecisionQUBO")],
 )[
diff --git a/docs/paper/references.bib b/docs/paper/references.bib
index 0dda81fd2..678d6a0b5 100644
--- a/docs/paper/references.bib
+++ b/docs/paper/references.bib
@@ -2230,6 +2230,7 @@ @misc{borosgruber2014quadratization
   eprint = {1404.6538},
   archivePrefix = {arXiv},
   primaryClass = {math.OC},
+  url = {https://arxiv.org/abs/1404.6538},
   note = {Paper presented at ISAIM 2012; arXiv version 2014},
 }
 
diff --git a/src/rules/comparativecontainment_qubo.rs b/src/rules/comparativecontainment_qubo.rs
new file mode 100644
index 000000000..0921f8268
--- /dev/null
+++ b/src/rules/comparativecontainment_qubo.rs
@@ -0,0 +1,169 @@
+//! Unit containment indicators encoded by penalized binary product chains.
+use crate::models::{ComparativeContainment, Decision, QUBO};
+use crate::rules::qubo_helpers::QuboBuilder;
+use crate::rules::{ReduceTo, ReductionError, ReductionResult};
+use crate::types::One;
+use std::collections::HashSet;
+
+type Builder = QuboBuilder<ComparativeContainment<One>, Decision<QUBO<i64>>>;
+
+#[derive(Debug, Clone)]
+pub struct ReductionComparativeContainmentToQUBO {
+    target: Decision<QUBO<i64>>,
+    universe_size: usize,
+}
+
+impl ReductionResult for ReductionComparativeContainmentToQUBO {
+    type Source = ComparativeContainment<One>;
+    type Target = Decision<QUBO<i64>>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    fn extract_solution(&self, bits: &Vec<bool>) -> crate::rules::ExtractionResult<Vec<bool>> {
+        crate::rules::traits::validate_target_witness(
+            &self.target,
+            bits,
+            |value| value.0,
+            "QUBO does not meet the containment threshold",
+        )?;
+        Ok(bits[..self.universe_size].to_vec())
+    }
+}
+
+#[crate::aggregate_reduction(identity)]
+impl crate::rules::AggregateReductionResult for ReductionComparativeContainmentToQUBO {}
+
+#[crate::reduction(transform = upper_bound {
+    num_vars = "universe_size + universe_size * (num_r_sets + num_s_sets)",
+    num_quadratic_terms = "6 * universe_size * (num_r_sets + num_s_sets) + num_r_sets + num_s_sets",
+})]
+impl ReduceTo<Decision<QUBO<i64>>> for ComparativeContainment<One> {
+    type Result = ReductionComparativeContainmentToQUBO;
+
+    fn reduce_to(&self) -> Result<Self::Result, ReductionError> {
+        construct(self, i64::MAX)
+    }
+}
+
+fn checked_sizes(
+    n: usize,
+    r: usize,
+    s: usize,
+    g: usize,
+    limit: i64,
+) -> Result<(usize, i64), ReductionError> {
+    let count = r.checked_add(s).ok_or_else(Builder::overflow)?;
+    let num_vars = n.checked_add(g).ok_or_else(Builder::overflow)?;
+    let r = i64::try_from(r).map_err(|_| Builder::overflow())?;
+    let count = i64::try_from(count).map_err(|_| Builder::overflow())?;
+    let g = i64::try_from(g).map_err(|_| Builder::overflow())?;
+    let penalty = r.checked_add(1).ok_or_else(Builder::overflow)?;
+    // Multiply by the gate count first: a gate-free target needs no 8*M.
+    let certificate = penalty
+        .checked_mul(g)
+        .and_then(|v| v.checked_mul(8))
+        .and_then(|v| v.checked_add(count))
+        .ok_or_else(Builder::overflow)?;
+    if penalty > limit || certificate > limit {
+        return Err(Builder::overflow());
+    }
+    Ok((num_vars, penalty))
+}
+
+fn construct(
+    source: &ComparativeContainment<One>,
+    limit: i64,
+) -> Result<ReductionComparativeContainmentToQUBO, ReductionError> {
+    let n = source.universe_size();
+    let mut g = 0usize;
+    for entry in source.r_sets().iter().chain(source.s_sets()) {
+        let distinct: HashSet<_> = entry.iter().collect();
+        let k = n - distinct.len();
+        g = g
+            .checked_add(k.saturating_sub(1))
+            .ok_or_else(Builder::overflow)?;
+    }
+    let (num_vars, penalty) = checked_sizes(n, source.num_r_sets(), source.num_s_sets(), g, limit)?;
+    let p = i128::from(penalty);
+    let mut q = Builder::new(num_vars);
+    let mut next = n;
+    for (family, sign) in [(source.r_sets(), -1i128), (source.s_sets(), 1)] {
+        for entry in family {
+            let present: HashSet<_> = entry.iter().copied().collect();
+            let mut complement = (0..n).filter(|v| !present.contains(v));
+            let Some(u) = complement.next() else {
+                q.add_constant(sign)?;
+                continue;
+            };
+            let Some(v) = complement.next() else {
+                q.add_constant(sign)?;
+                q.add(u, u, -sign)?;
+                continue;
+            };
+            let mut z = next;
+            next = next.checked_add(1).ok_or_else(Builder::overflow)?;
+            // P(1-x_u,1-x_v,z) = 1-x_u-x_v+x_u*x_v-z+2*x_u*z+2*x_v*z.
+            q.add_constant(p)?;
+            for (a, b, c) in [
+                (u, u, -p),
+                (v, v, -p),
+                (u, v, p),
+                (z, z, -p),
+                (u, z, 2 * p),
+                (v, z, 2 * p),
+            ] {
+                q.add(a, b, c)?;
+            }
+            for v in complement {
+                let previous = z;
+                z = next;
+                next = next.checked_add(1).ok_or_else(Builder::overflow)?;
+                // P(w,1-x_v,z) = w-w*x_v-2*w*z+z+2*x_v*z.
+                for (a, b, c) in [
+                    (previous, previous, p),
+                    (previous, v, -p),
+                    (previous, z, -2 * p),
+                    (z, z, p),
+                    (v, z, 2 * p),
+                ] {
+                    q.add(a, b, c)?;
+                }
+            }
+            q.add(z, z, sign)?;
+        }
+    }
+    let (target, constant) = q.finish()?;
+    let bound = constant.checked_neg().ok_or_else(Builder::overflow)?;
+    Ok(ReductionComparativeContainmentToQUBO {
+        target: Decision::new(target, bound),
+        universe_size: n,
+    })
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "comparativecontainment_qubo",
+        build: || {
+            let source = ComparativeContainment::<One>::new(
+                4,
+                vec![vec![0, 1], vec![1, 2], vec![1, 3]],
+                vec![vec![0, 1, 2, 3]; 2],
+            )
+            .expect("canonical containment instance must be valid");
+            crate::example_db::specs::rule_example_with_witness::<_, Decision<QUBO<i64>>>(
+                source,
+                crate::export::SolutionPair {
+                    source_config: serde_json::json!([false, true, false, false]),
+                    target_config: serde_json::json!([false, true, false, false, true, true, true]),
+                },
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/comparativecontainment_qubo.rs"]
+mod tests;
diff --git a/src/rules/mod.rs b/src/rules/mod.rs
index 7606f4923..ccf94474d 100644
--- a/src/rules/mod.rs
+++ b/src/rules/mod.rs
@@ -45,6 +45,7 @@ pub(crate) mod twodimensionalconsecutivesets_ilp;
 pub(crate) use graph::{RecoveryStep, recover_completed_result};
 mod binpacking_casts;
 mod comparativecontainment_casts;
+pub(crate) mod comparativecontainment_qubo;
 pub(crate) mod graph_helpers;
 pub(crate) mod graphpartitioning_maxcut;
 pub(crate) mod graphpartitioning_qubo;
@@ -498,6 +499,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
         optimallineararrangement_sequencingtominimizeweightedcompletiontime::canonical_rule_example_specs(),
     );
     specs.extend(setsplitting_betweenness::canonical_rule_example_specs());
+    specs.extend(comparativecontainment_qubo::canonical_rule_example_specs());
     specs.extend(setsplitting_naesatisfiability::canonical_rule_example_specs());
     specs.extend(satisfiability_integralflowhomologousarcs::canonical_rule_example_specs());
     specs.extend(satisfiability_naesatisfiability::canonical_rule_example_specs());
diff --git a/src/unit_tests/rules/comparativecontainment_qubo.rs b/src/unit_tests/rules/comparativecontainment_qubo.rs
new file mode 100644
index 000000000..ac07e1c62
--- /dev/null
+++ b/src/unit_tests/rules/comparativecontainment_qubo.rs
@@ -0,0 +1,359 @@
+use super::*;
+use crate::rules::AggregateReductionResult;
+use crate::rules::test_helpers::{
+    assert_parameter_predictions, assert_satisfaction_round_trip_from_satisfaction_target,
+};
+use crate::solvers::BruteForce;
+use crate::traits::Problem;
+use crate::types::{Min, Or};
+
+fn example(s: usize) -> ComparativeContainment<One> {
+    ComparativeContainment::new(
+        4,
+        vec![vec![0, 1], vec![1, 2], vec![1, 3]],
+        vec![vec![0, 1, 2, 3]; s],
+    )
+    .unwrap()
+}
+
+#[test]
+fn test_comparativecontainment_to_qubo_closed_loop() {
+    for s in [2, 4] {
+        let source = example(s);
+        let reduction = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source).unwrap();
+        if s == 2 {
+            assert_satisfaction_round_trip_from_satisfaction_target(
+                &source,
+                &reduction,
+                "containment",
+            );
+        }
+        assert_parameter_predictions(&source, &reduction);
+        assert_eq!(reduction.extract_value(Or(s == 2)), Or(s == 2));
+        assert_eq!(BruteForce::new().solve(&source).unwrap().is_some(), s == 2);
+        assert_eq!(
+            BruteForce::new()
+                .solve(ReductionResult::target_problem(&reduction))
+                .unwrap()
+                .is_some(),
+            s == 2
+        );
+    }
+}
+
+#[test]
+fn test_canonical_polynomial_and_rejected_candidates() {
+    let reduction = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&example(2)).unwrap();
+    let target = ReductionResult::target_problem(&reduction);
+    assert_eq!(target.inner().num_vars(), 7);
+    assert_eq!(*target.bound(), -14);
+    assert_eq!(target.inner().num_quadratic_terms(), 9);
+    assert_eq!(
+        target.inner().entries(),
+        &[
+            (0, 0, -8),
+            (0, 2, 4),
+            (0, 3, 4),
+            (0, 5, 8),
+            (0, 6, 8),
+            (2, 2, -8),
+            (2, 3, 4),
+            (2, 4, 8),
+            (2, 6, 8),
+            (3, 3, -8),
+            (3, 4, 8),
+            (3, 5, 8),
+            (4, 4, -5),
+            (5, 5, -5),
+            (6, 6, -5),
+        ]
+    );
+    let witness = vec![false, true, false, false, true, true, true];
+    assert_eq!(target.inner().evaluate(&witness).unwrap(), Min(Some(-15)));
+    assert_eq!(reduction.extract_solution(&witness).unwrap(), witness[..4]);
+    let mut faulty = witness.clone();
+    faulty[4] = false;
+    assert_eq!(target.inner().evaluate(&faulty).unwrap(), Min(Some(-10)));
+    assert!(reduction.extract_solution(&faulty).is_err());
+    let candidate = vec![true, true, false, false, true, false, false];
+    assert_eq!(target.inner().evaluate(&candidate).unwrap(), Min(Some(-13)));
+    assert!(reduction.extract_solution(&candidate).is_err());
+    for bad in [vec![], vec![false; 6], vec![false; 8]] {
+        assert!(reduction.extract_solution(&bad).is_err());
+    }
+    let no = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&example(4)).unwrap();
+    assert_eq!(*ReductionResult::target_problem(&no).bound(), -16);
+    assert!(no.extract_solution(&witness).is_err());
+}
+
+fn bits(mask: usize, n: usize) -> Vec<bool> {
+    (0..n).map(|i| mask & (1 << i) != 0).collect()
+}
+
+fn check_instance(source: ComparativeContainment<One>) {
+    let n = source.universe_size();
+    let reduction = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source).unwrap();
+    let target = ReductionResult::target_problem(&reduction);
+    let g: usize = source
+        .r_sets()
+        .iter()
+        .chain(source.s_sets())
+        .map(|a| {
+            let k = (0..n).filter(|v| !a.contains(v)).count();
+            k.saturating_sub(1)
+        })
+        .sum();
+    assert_eq!(target.inner().num_vars(), n + g);
+    assert!(target.inner().num_quadratic_terms() <= 3 * g);
+    assert!(target.inner().entries().len() <= 6 * g + source.num_r_sets() + source.num_s_sets());
+    let mut source_yes = false;
+    for mask in 0..(1 << n) {
+        let x = bits(mask, n);
+        let mut extension = x.clone();
+        for a in source.r_sets().iter().chain(source.s_sets()) {
+            let outside: Vec<_> = (0..n).filter(|v| !a.contains(v)).collect();
+            if outside.len() >= 2 {
+                let mut prefix = !x[outside[0]];
+                for &v in &outside[1..] {
+                    prefix &= !x[v];
+                    extension.push(prefix);
+                }
+            }
+        }
+        // Source sums are evaluated by the model, independently of QUBO gates.
+        let valid = source.evaluate(&x).unwrap().0;
+        source_yes |= valid;
+        let difference =
+            source.s_weight_sum(&x).unwrap().unwrap() - source.r_weight_sum(&x).unwrap().unwrap();
+        assert_eq!(
+            target.inner().evaluate(&extension).unwrap().unwrap() - target.bound(),
+            difference
+        );
+        assert_eq!(target.evaluate(&extension).unwrap().0, valid);
+    }
+    let mut target_yes = false;
+    for mask in 0..(1 << target.inner().num_vars()) {
+        let y = bits(mask, target.inner().num_vars());
+        if target.evaluate(&y).unwrap().0 {
+            target_yes = true;
+            let x = reduction.extract_solution(&y).unwrap();
+            assert_eq!(x, y[..n]);
+            assert!(source.evaluate(&x).unwrap().0);
+        } else {
+            assert!(reduction.extract_solution(&y).is_err());
+        }
+    }
+    assert_eq!(source_yes, target_yes);
+}
+
+#[test]
+fn test_exhaustive_two_element_family_pairs() {
+    for n in 0..=2 {
+        let sets: Vec<Vec<_>> = (0..1 << n)
+            .map(|mask| (0..n).filter(|v| mask & (1 << v) != 0).collect())
+            .collect();
+        let mut families = vec![vec![]];
+        for a in &sets {
+            families.push(vec![a.clone()]);
+            for b in &sets {
+                families.push(vec![a.clone(), b.clone()]);
+            }
+        }
+        for r in &families {
+            for s in &families {
+                check_instance(ComparativeContainment::new(n, r.clone(), s.clone()).unwrap());
+            }
+        }
+    }
+}
+
+fn check_all_two_entries(n: usize) {
+    let sets: Vec<Vec<_>> = (0..1 << n)
+        .map(|mask| (0..n).filter(|v| mask & (1 << v) != 0).collect())
+        .collect();
+    for r in 0..=2 {
+        for a in &sets {
+            for b in &sets {
+                let entries = [a.clone(), b.clone()];
+                check_instance(
+                    ComparativeContainment::new(n, entries[..r].to_vec(), entries[r..].to_vec())
+                        .unwrap(),
+                );
+            }
+        }
+    }
+}
+
+#[test]
+fn test_exhaustive_three_element_two_entry_families() {
+    check_all_two_entries(3);
+}
+#[test]
+fn test_exhaustive_four_element_two_entry_families() {
+    check_all_two_entries(4);
+}
+#[test]
+fn test_exhaustive_five_element_two_entry_families() {
+    check_all_two_entries(5);
+}
+
+#[test]
+fn test_empty_full_duplicate_and_cancelled_sets() {
+    for n in 0..=5 {
+        let full: Vec<_> = (0..n).collect();
+        for (r, s) in [
+            (vec![], vec![]),
+            (vec![], vec![full.clone()]),
+            (vec![vec![]], vec![]),
+            (vec![full.clone()], vec![]),
+            (vec![full.clone()], vec![full.clone()]),
+            (vec![vec![], vec![]], vec![full.clone()]),
+            (vec![full.clone(); 2], vec![full.clone(); 3]),
+        ] {
+            let source = ComparativeContainment::<One>::new(n, r, s).unwrap();
+            let reduction = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source).unwrap();
+            assert_parameter_predictions(&source, &reduction);
+            check_instance(source);
+        }
+    }
+    let a = ComparativeContainment::<One>::new(3, vec![vec![0, 2]], vec![vec![1]]).unwrap();
+    let b =
+        ComparativeContainment::<One>::new(3, vec![vec![2, 0, 2, 0]], vec![vec![1, 1]]).unwrap();
+    let ta = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&a).unwrap();
+    let tb = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&b).unwrap();
+    assert_eq!(
+        serde_json::to_value(ReductionResult::target_problem(&ta)).unwrap(),
+        serde_json::to_value(ReductionResult::target_problem(&tb)).unwrap()
+    );
+    check_instance(b);
+    let source = ComparativeContainment::<One>::new(3, vec![vec![0, 1]], vec![vec![0, 1]]).unwrap();
+    let reduced = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source).unwrap();
+    let target = ReductionResult::target_problem(&reduced);
+    assert_eq!(target.inner().num_vars(), 3);
+    assert!(target.inner().entries().is_empty());
+    assert_eq!(*target.bound(), 0);
+}
+
+#[test]
+fn test_numeric_certificate_and_index_boundaries() {
+    let limit = i64::MAX;
+    for (r, s) in [(0, 0), (1, 0), (2, 3)] {
+        let g = usize::try_from((limit - r - s) / (8 * (r + 1))).unwrap();
+        let (nvars, penalty) = checked_sizes(
+            0,
+            usize::try_from(r).unwrap(),
+            usize::try_from(s).unwrap(),
+            g,
+            limit,
+        )
+        .unwrap();
+        assert_eq!((nvars, penalty), (g, r + 1));
+        assert!(matches!(
+            checked_sizes(
+                0,
+                usize::try_from(r).unwrap(),
+                usize::try_from(s).unwrap(),
+                g + 1,
+                limit
+            ),
+            Err(ReductionError::IntegerOverflow { .. })
+        ));
+    }
+    let max = usize::try_from(limit).unwrap();
+    assert_eq!(checked_sizes(0, max - 1, 1, 0, limit).unwrap(), (0, limit));
+    assert_eq!(checked_sizes(0, 0, max, 0, limit).unwrap(), (0, 1));
+    for counts in [
+        (0, max, 0, 0),
+        (0, max + 1, 0, 0),
+        (0, max - 1, 2, 0),
+        (0, 0, max + 1, 0),
+        (usize::MAX, 0, 0, 1),
+        (0, usize::MAX, 1, 0),
+        (0, 0, 0, usize::MAX),
+    ] {
+        let (n, r, s, g) = counts;
+        assert!(matches!(
+            checked_sizes(n, r, s, g, limit),
+            Err(ReductionError::IntegerOverflow { .. })
+        ));
+    }
+    assert_eq!(
+        checked_sizes(usize::MAX - 1, 0, 0, 1, limit).unwrap(),
+        (usize::MAX, 1)
+    );
+    // A large penalty with zero gates must not compute 8*M first.
+    assert!(checked_sizes(0, max - 1, 0, 0, limit).is_ok());
+}
+
+#[test]
+fn test_constructed_target_at_coefficient_limit() {
+    let empty = ComparativeContainment::<One>::new(0, vec![], vec![]).unwrap();
+    let zero = construct(&empty, 1).unwrap();
+    assert_eq!(
+        ReductionResult::target_problem(&zero)
+            .evaluate(&vec![])
+            .unwrap(),
+        Or(true)
+    );
+    assert!(matches!(
+        construct(&empty, 0),
+        Err(ReductionError::IntegerOverflow { .. })
+    ));
+
+    let mut r = vec![vec![0, 1]; 9];
+    r.push(vec![]);
+    let source = ComparativeContainment::<One>::new(2, r, vec![]).unwrap();
+    // r=10, M=11, g=1: 8*M*g+r+s=98, while an actual coefficient is 22.
+    let reduction = construct(&source, 98).unwrap();
+    let target = ReductionResult::target_problem(&reduction);
+    assert_eq!(*target.bound(), -2);
+    assert_eq!(
+        target.inner().entries(),
+        &[
+            (0, 0, -11),
+            (0, 1, 11),
+            (0, 2, 22),
+            (1, 1, -11),
+            (1, 2, 22),
+            (2, 2, -12),
+        ]
+    );
+    assert!(matches!(
+        construct(&source, 97),
+        Err(ReductionError::IntegerOverflow { .. })
+    ));
+    for mask in 0..8 {
+        let y = bits(mask, 3);
+        assert_eq!(
+            target.evaluate(&y).unwrap(),
+            ReductionResult::target_problem(
+                &ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source).unwrap()
+            )
+            .evaluate(&y)
+            .unwrap()
+        );
+    }
+}
+
+#[test]
+fn test_sparse_dimension_and_gate_count_overflow() {
+    let empty = ComparativeContainment::<One>::new(usize::MAX, vec![], vec![]).unwrap();
+    let reduced = ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&empty).unwrap();
+    assert_eq!(
+        ReductionResult::target_problem(&reduced).inner().num_vars(),
+        usize::MAX
+    );
+    assert!(
+        ReductionResult::target_problem(&reduced)
+            .inner()
+            .entries()
+            .is_empty()
+    );
+    for r in [vec![vec![]], vec![vec![], vec![]]] {
+        let source = ComparativeContainment::<One>::new(usize::MAX, r, vec![]).unwrap();
+        assert!(matches!(
+            ReduceTo::<Decision<QUBO<i64>>>::reduce_to(&source),
+            Err(ReductionError::IntegerOverflow { .. })
+        ));
+    }
+}