diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ
index 8e0462b7a..608dcdb02 100644
--- a/docs/paper/reductions.typ
+++ b/docs/paper/reductions.typ
@@ -13540,6 +13540,152 @@ These constructions use the checked sparse penalty and extraction contracts abov
   caption: [Ising gadgets for logic gates. Ground states match truth tables.]
 ) <tab:gadgets>
 
+// SubsetProduct → CircuitSAT
+#let sp_cs = load-example("SubsetProduct", "CircuitSAT")
+#let sp_cs_sol = sp_cs.solutions.at(0)
+#let sp_sizes = sp_cs.source.instance.sizes.map(int)
+#let sp_target = int(sp_cs.source.instance.target)
+#let sp_bits(a) = {
+  let bits = ()
+  while a > 0 { bits.push(calc.rem(a, 2)); a = calc.floor(a / 2) }
+  bits
+}
+#let sp_n = sp_sizes.len()
+#let sp_d = sp_sizes.map(a => sp_bits(a).len()).sum()
+#let sp_h = sp_sizes.map(a => sp_bits(a).sum()).sum()
+#let sp_w = calc.max(1 + sp_d, sp_bits(sp_target).len())
+#let sp_vars = circuit-num-variables(sp_cs.target.instance)
+#let sp_assignments = sp_cs.target.instance.circuit.assignments
+#let sp_nodes(expr) = {
+  let (kind, value) = expr.op.pairs().first()
+  if kind == "Var" or kind == "Const" { 1 }
+  else if kind == "Not" { 1 + sp_nodes(value) }
+  else { 1 + value.map(sp_nodes).sum(default: 0) }
+}
+#let sp_word(value) = range(sp_w).map(b => str(calc.rem(calc.floor(value / calc.pow(2, b)), 2))).join(", ")
+#let sp_node_count = sp_assignments.map(a => sp_nodes(a.expr)).sum()
+#let sp_output_count = sp_assignments.map(a => a.outputs.len()).sum()
+#let sp_prefix = (1,)
+#for (i, a) in sp_sizes.enumerate() {
+  sp_prefix.push(sp_prefix.last() * if sp_cs_sol.source_config.at(i) { a } else { 1 })
+}
+#let sp_wire(name) = {
+  let coordinate = sp_cs.target.instance.variables.position(v => v == name)
+  assert(coordinate != none)
+  sp_cs_sol.target_config.at(coordinate)
+}
+#for (i, expected) in sp_prefix.enumerate() {
+  let actual = range(sp_w).map(b => if sp_wire("p_" + str(i) + "_" + str(b)) { calc.pow(2, b) } else { 0 }).sum()
+  assert(actual == expected)
+}
+#let sp_products = (1,)
+#for a in sp_sizes { sp_products += sp_products.map(p => p * a) }
+#let sp_no = range(sp_sizes.sorted().last() + 1, sp_target).find(b => not sp_products.contains(b))
+#assert(sp_vars == sp_w * (1 + sp_n + 5 * sp_h) + sp_n)
+#assert(sp_assignments.len() == sp_w * (2 + sp_n + 5 * sp_h))
+#assert(sp_output_count == sp_assignments.len())
+#assert(sp_node_count == sp_w * (2 + 8 * sp_n + 15 * sp_h))
+#assert(sp_prefix.last() == sp_target)
+#assert(sp_no > 0 and not sp_products.contains(sp_no))
+#assert(calc.max(1 + sp_d, sp_bits(sp_no).len()) == sp_w)
+#reduction-rule("SubsetProduct", "CircuitSAT",
+  example: true,
+  example-caption: [Select factors from $(#sp_sizes.map(str).join(", "))$ to obtain $#sp_target$],
+  extra: [
+    #pred-commands(
+      "pred create --example " + rule-spec(sp_cs) + " -o subset-product.json",
+      "pred path SubsetProduct CircuitSAT -o paths.json",
+      "jq '.paths[] | select(.path | length == 1)' paths.json > route.json",
+      "pred reduce subset-product.json --via route.json -o bundle.json",
+      "pred solve bundle.json",
+      "pred evaluate subset-product.json --config " + cli-config(sp_cs_sol.source_config),
+    )
+    *Worked YES example.* The fixture has sizes $(#sp_sizes.map(str).join(", "))$
+    and target $#sp_target$. Thus $D=#sp_d$, $H=#sp_h$, and $W=#sp_w$.
+    There are $#sp_n$ selectors, $#(sp_n + 1)$ product words,
+    $#(5 * sp_w * sp_h)$ adder outputs, $#(sp_n * sp_w)$ multiplexers,
+    $#sp_w$ initial constants and $#sp_w$ final pins. The exported target has
+    $#sp_vars$ variables, $#sp_assignments.len()$ assignments,
+    $#sp_output_count$ outputs and $#sp_node_count$ expression nodes.
+
+    The canonical selectors $(#fmt-values(sp_cs_sol.source_config))$ give the
+    prefix products $(#sp_prefix.map(str).join(", "))$. At each stage:
+    #table(columns: 4, [Factor], [Incoming word], [Shifted summands], [Selected word],
+      ..sp_sizes.enumerate().map(((i, a)) => {
+        let summands = sp_bits(a).enumerate().filter(((j, b)) => b == 1)
+          .map(((j, b)) => sp_prefix.at(i) * calc.pow(2, j))
+        assert(summands.sum() == sp_prefix.at(i) * a)
+        (str(a), str(sp_prefix.at(i)), summands.map(str).join(" + "), str(sp_prefix.at(i + 1)))
+      }).flatten())
+    The final little-endian word is $(#sp_word(sp_target))$.
+    Named extraction returns the canonical selectors, whose product is $#sp_target$.
+
+    *Worked NO example.* Keep the fixture's sizes and replace the target by
+    $B=#sp_no$. The width and all four counts remain unchanged. The complete
+    subset-product list is $(#sp_products.sorted().map(str).join(", "))$;
+    none equals $#sp_no$. The final pins demand $(#sp_word(sp_no))$.
+    By backward correctness there can be no satisfying target assignment.
+  ],
+)[
+  Positive arbitrary-precision subset multiplication reduces to simultaneous
+  Boolean equations by a cascade of constant multipliers and selectors. Binary
+  summands and adders are classical arithmetic circuitry @wallace1964fastmultiplier;
+  the ripple cascade and decision-equivalence argument below give this construction.
+  It uses polynomial space and time in the binary input length, without factoring.
+][
+  _Construction._ For positive sizes $a_1, dots, a_n$ and target $B$, define
+  $ell(a)=floor(log_2(a))+1$, $D=sum_i ell(a_i)$, $H=sum_i "popcount"(a_i)$,
+  and $W=max(1+D,ell(B))$. Words are $W$ bits in little-endian order.
+  Introduce free selectors $x_0, dots, x_(n-1)$ and pin $P_0=1$.
+  For each stage $i$, initialize an accumulator to zero. For every set-bit
+  position $j$ of $a_i$, in increasing order, add $P_(i-1) 2^j$ to it with
+  a fresh ripple row and incoming carry zero. Shifted columns below $j$ are
+  literal false. For each column with accumulator bit $u$, shifted bit $v$
+  and carry $c$, emit five fresh equations:
+  $ h=u xor v, quad s=h xor c, quad f=u and v, quad g=h and c, quad d=f or g. $
+  The next column uses carry $d$; retain even the last carry equation.
+  The row's sum bits become the accumulator. Call the final word $T_i$ and emit
+  $P_(i,b)=(x_(i-1) and T_(i,b)) or ((not x_(i-1)) and P_(i-1,b))$.
+  Finally append $W$ constant equations pinning $P_n=B$, reusing its names.
+  Product names are `p_i_b`; adder names include stage, shift, column and gate.
+  The disjoint name families prevent collisions. Repeated final outputs are
+  simultaneous constraints, not overwrites.
+
+  _Correctness._ Each full adder obeys $u+v+c=s+2d$ on all eight input triples,
+  with uniquely determined outputs. Hence a ripple row computes addition modulo
+  $2^W$. Let $Q_i=product_(t=1)^i a_t$. Positivity bounds any selected prefix
+  product by $Q_i$. At stage $i$, every shifted summand and partial accumulation
+  is at most $P_(i-1)a_i<=Q_i<=Q_n<2^D<2^W$ when $n>0$. Thus no shift loses
+  a high bit and every discarded carry is zero. This bound includes branches
+  discarded by selectors. With no items, $1<2^W$ and the repeated pins require
+  exactly $B=1$.
+
+  ($arrow.r.double$) Given a feasible subset, assign its selectors and evaluate
+  fresh gates in construction order. Exact addition and the multiplexers force
+  $P_i=product_(t<=i, x_(t-1)=1) a_t$, so every final pin holds.
+  ($arrow.l.double$) Any satisfying assignment must obey the initial pins and
+  every deterministic gate equation. Induction forces the same exact prefix
+  products, even for unselected multiplication branches. The final pins compare
+  all $W>=ell(B)$ bits, so the extracted subset product is exactly $B$.
+  This also proves unique extension of each feasible selector vector. Ones,
+  repeated items, single items, full-product targets and arbitrarily wide
+  magnitudes require no special cases. Targets wider than $D$ cannot wrap around.
+
+  _Solution extraction._ Validate the complete target witness once, rejecting
+  wrong lengths and unsatisfied equations. Read selector coordinates by names
+  `x_0`, …, `x_(n-1)` in source order; lexicographic variable order is not numeric
+  selector order. Missing coordinates are errors. Completed `Or` values map identically.
+
+  _Size._ The exact counts are $W(1+n+5H)+n$ variables, $W(2+n+5H)$ assignments
+  and outputs, and $W(2+8n+15H)$ expression nodes. Only `num_elements` is a
+  registered source parameter, so all four parameter transforms are unavailable:
+  bit lengths and popcounts cannot be bounded in terms of that parameter alone.
+  For $L=D+ell(B)+n+1$, the circuit has $O(L^2)$ expression nodes; its named
+  encoding takes $O(L^2 log L)$ bits and comparison sorting gives a conservative
+  $O(L^2 log^2 L)$ construction bound. Machine counts, index conversions and
+  allocation sizes are checked; failures are construction errors, never NO answers.
+]
+
 #let fact_cs = load-example("Factoring", "CircuitSAT")
 #let fact_cs_sol = fact_cs.solutions.at(0)
 #let fact-nbf = fact_cs.source.instance.m
diff --git a/docs/paper/references.bib b/docs/paper/references.bib
index 0dda81fd2..123eba49b 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{wallace1964fastmultiplier,
+  author = {Wallace, C. S.},
+  title = {A Suggestion for a Fast Multiplier},
+  journal = {IEEE Transactions on Electronic Computers},
+  volume = {EC-13},
+  number = {1},
+  pages = {14--17},
+  year = {1964},
+  doi = {10.1109/PGEC.1964.263830}
+}
diff --git a/src/models/formula/circuit.rs b/src/models/formula/circuit.rs
index 2e078c074..3006fcf1a 100644
--- a/src/models/formula/circuit.rs
+++ b/src/models/formula/circuit.rs
@@ -254,8 +254,59 @@ pub struct CircuitSAT {
 impl CircuitSAT {
     /// Create a new CircuitSAT problem.
     pub fn new(circuit: Circuit) -> Self {
-        let variables = circuit.variables();
-        Self { circuit, variables }
+        Self::try_new(circuit).unwrap_or_else(|error| panic!("{error}"))
+    }
+
+    /// Create a CircuitSAT problem, reporting failures to allocate its variable table.
+    /// Repeated outputs remain simultaneous equations, and names are sorted.
+    pub fn try_new(circuit: Circuit) -> Result<Self, crate::registry::ConstructionError> {
+        fn allocation(
+            error: std::collections::TryReserveError,
+        ) -> crate::registry::ConstructionError {
+            crate::registry::ConstructionError::Conversion(format!(
+                "cannot allocate circuit variables: {error}"
+            ))
+        }
+        fn collect<'a>(
+            expr: &'a BooleanExpr,
+            names: &mut Vec<&'a str>,
+        ) -> Result<(), crate::registry::ConstructionError> {
+            match &expr.op {
+                BooleanOp::Var(name) => {
+                    names.try_reserve(1).map_err(allocation)?;
+                    names.push(name);
+                }
+                BooleanOp::Const(_) => {}
+                BooleanOp::Not(inner) => collect(inner, names)?,
+                BooleanOp::And(args) | BooleanOp::Or(args) | BooleanOp::Xor(args) => {
+                    for arg in args {
+                        collect(arg, names)?;
+                    }
+                }
+            }
+            Ok(())
+        }
+        let mut names = Vec::new();
+        for assignment in &circuit.assignments {
+            names
+                .try_reserve(assignment.outputs.len())
+                .map_err(allocation)?;
+            names.extend(assignment.outputs.iter().map(String::as_str));
+            collect(&assignment.expr, &mut names)?;
+        }
+        names.sort_unstable();
+        names.dedup();
+        let mut variables = Vec::new();
+        variables
+            .try_reserve_exact(names.len())
+            .map_err(allocation)?;
+        for name in names {
+            let mut owned = String::new();
+            owned.try_reserve_exact(name.len()).map_err(allocation)?;
+            owned.push_str(name);
+            variables.push(owned);
+        }
+        Ok(Self { circuit, variables })
     }
 
     /// Get the circuit.
diff --git a/src/rules/mod.rs b/src/rules/mod.rs
index 7606f4923..67ee25cf8 100644
--- a/src/rules/mod.rs
+++ b/src/rules/mod.rs
@@ -172,6 +172,7 @@ mod spinglass_casts;
 pub(crate) mod spinglass_maxcut;
 pub(crate) mod spinglass_qubo;
 mod steinertree_casts;
+pub(crate) mod subsetproduct_circuitsat;
 pub(crate) mod subsetsum_closestvectorproblem;
 pub(crate) mod subsetsum_integerexpressionmembership;
 pub(crate) mod subsetsum_partition;
@@ -502,6 +503,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
     specs.extend(satisfiability_integralflowhomologousarcs::canonical_rule_example_specs());
     specs.extend(satisfiability_naesatisfiability::canonical_rule_example_specs());
     specs.extend(sat_circuitsat::canonical_rule_example_specs());
+    specs.extend(subsetproduct_circuitsat::canonical_rule_example_specs());
     specs.extend(sat_coloring::canonical_rule_example_specs());
     specs.extend(sat_ksat::canonical_rule_example_specs());
     specs.extend(sat_maximumindependentset::canonical_rule_example_specs());
diff --git a/src/rules/subsetproduct_circuitsat.rs b/src/rules/subsetproduct_circuitsat.rs
new file mode 100644
index 000000000..f8a2dc6bf
--- /dev/null
+++ b/src/rules/subsetproduct_circuitsat.rs
@@ -0,0 +1,234 @@
+//! Reduction from SubsetProduct to CircuitSAT by exact binary multiplication.
+
+use crate::models::formula::{Assignment, BooleanExpr, Circuit};
+use crate::models::{CircuitSAT, SubsetProduct};
+use crate::reduction;
+use crate::rules::{ExtractionError, ReduceTo, ReductionError, ReductionResult};
+
+/// Circuit and named selector coordinates for subset-product extraction.
+#[derive(Debug, Clone)]
+pub struct ReductionSubsetProductToCircuitSAT {
+    target: CircuitSAT,
+    source_var_indices: Vec<usize>,
+}
+
+impl ReductionResult for ReductionSubsetProductToCircuitSAT {
+    type Source = SubsetProduct;
+    type Target = CircuitSAT;
+
+    fn target_problem(&self) -> &CircuitSAT {
+        &self.target
+    }
+
+    fn extract_solution(&self, solution: &Vec<bool>) -> crate::rules::ExtractionResult<Vec<bool>> {
+        crate::rules::traits::validate_target_witness(
+            &self.target,
+            solution,
+            |value| value.0,
+            "target circuit is not satisfied",
+        )?;
+        self.source_var_indices
+            .iter()
+            .map(|&index| {
+                solution
+                    .get(index)
+                    .copied()
+                    .ok_or_else(|| ExtractionError::invalid("missing selector coordinate"))
+            })
+            .collect()
+    }
+}
+
+#[crate::aggregate_reduction(identity)]
+impl crate::rules::AggregateReductionResult for ReductionSubsetProductToCircuitSAT {}
+
+fn checked<T>(value: Option<T>) -> Result<T, ReductionError> {
+    value.ok_or_else(|| {
+        ReductionError::integer_overflow::<SubsetProduct, CircuitSAT>(
+            "sizing the multiplication circuit",
+        )
+    })
+}
+
+fn sum_sizes(values: impl IntoIterator<Item = usize>) -> Result<usize, ReductionError> {
+    values
+        .into_iter()
+        .try_fold(0usize, |sum, value| checked(sum.checked_add(value)))
+}
+
+fn index(value: u64) -> Result<usize, ReductionError> {
+    checked(usize::try_from(value).ok())
+}
+
+// Width, variables, assignments (also outputs), and expression nodes.
+fn dimensions(
+    n: usize,
+    d: usize,
+    h: usize,
+    target_bits: usize,
+) -> Result<(usize, usize, usize, usize), ReductionError> {
+    let w = checked(d.checked_add(1))?.max(target_bits);
+    let gates = checked(h.checked_mul(5))?;
+    let variables = checked(checked(w.checked_mul(sum_sizes([1, n, gates])?))?.checked_add(n))?;
+    let assignments = checked(w.checked_mul(sum_sizes([2, n, gates])?))?;
+    let nodes = checked(w.checked_mul(sum_sizes([
+        2,
+        checked(n.checked_mul(8))?,
+        checked(h.checked_mul(15))?,
+    ])?))?;
+    checked(assignments.checked_mul(std::mem::size_of::<Assignment>()))?;
+    checked(nodes.checked_mul(std::mem::size_of::<BooleanExpr>()))?;
+    checked(variables.checked_mul(std::mem::size_of::<String>()))?;
+    Ok((w, variables, assignments, nodes))
+}
+
+fn reserve_assignments(count: usize) -> Result<Vec<Assignment>, ReductionError> {
+    let mut assignments = Vec::new();
+    assignments.try_reserve_exact(count).map_err(|error| {
+        <SubsetProduct as ReduceTo<CircuitSAT>>::target_construction(
+            crate::registry::ConstructionError::Conversion(format!(
+                "cannot allocate circuit assignments: {error}"
+            )),
+        )
+    })?;
+    Ok(assignments)
+}
+
+#[reduction(transform = unavailable {
+    num_variables = "Integer bit lengths and popcounts are not bounded by num_elements",
+    num_assignments = "Integer bit lengths and popcounts are not bounded by num_elements",
+    num_assignment_outputs = "Integer bit lengths and popcounts are not bounded by num_elements",
+    num_expression_nodes = "Integer bit lengths and popcounts are not bounded by num_elements",
+})]
+impl ReduceTo<CircuitSAT> for SubsetProduct {
+    type Result = ReductionSubsetProductToCircuitSAT;
+
+    fn reduce_to(&self) -> Result<Self::Result, ReductionError> {
+        let d = self
+            .sizes()
+            .iter()
+            .try_fold(0usize, |sum, a| checked(sum.checked_add(index(a.bits())?)))?;
+        let h = self.sizes().iter().try_fold(0usize, |sum, a| {
+            checked(sum.checked_add(index(a.count_ones())?))
+        })?;
+        let (w, _, count, _) = dimensions(self.sizes().len(), d, h, index(self.target().bits())?)?;
+        let mut assignments = reserve_assignments(count)?;
+        for b in 0..w {
+            assignments.push(Assignment::new(
+                vec![format!("p_0_{b}")],
+                BooleanExpr::constant(b == 0),
+            ));
+        }
+        for (i, a) in self.sizes().iter().enumerate() {
+            // The width bounds every branch, including products discarded by a selector.
+            let stage = checked(i.checked_add(1))?;
+            let mut accumulator = vec![BooleanExpr::constant(false); w];
+            for shift in (0..a.bits()).filter(|&shift| a.bit(shift)) {
+                let j = index(shift)?;
+                let mut carry = BooleanExpr::constant(false);
+                for (b, sum) in accumulator.iter_mut().enumerate() {
+                    let u = sum.clone();
+                    let v = if b >= j {
+                        BooleanExpr::var(&format!("p_{i}_{}", b - j))
+                    } else {
+                        BooleanExpr::constant(false)
+                    };
+                    let name = |kind| format!("a_{stage}_{j}_{b}_{kind}");
+                    let h = BooleanExpr::var(&name('h'));
+                    let f = BooleanExpr::var(&name('f'));
+                    let g = BooleanExpr::var(&name('g'));
+                    for (kind, expr) in [
+                        ('h', BooleanExpr::xor(vec![u.clone(), v.clone()])),
+                        ('s', BooleanExpr::xor(vec![h.clone(), carry.clone()])),
+                        ('f', BooleanExpr::and(vec![u, v])),
+                        ('g', BooleanExpr::and(vec![h, carry])),
+                        ('d', BooleanExpr::or(vec![f, g])),
+                    ] {
+                        assignments.push(Assignment::new(vec![name(kind)], expr));
+                    }
+                    *sum = BooleanExpr::var(&name('s'));
+                    carry = BooleanExpr::var(&name('d'));
+                }
+            }
+            for (b, product) in accumulator.into_iter().enumerate() {
+                let selector = BooleanExpr::var(&format!("x_{i}"));
+                assignments.push(Assignment::new(
+                    vec![format!("p_{stage}_{b}")],
+                    BooleanExpr::or(vec![
+                        BooleanExpr::and(vec![selector.clone(), product]),
+                        BooleanExpr::and(vec![
+                            BooleanExpr::not(selector),
+                            BooleanExpr::var(&format!("p_{i}_{b}")),
+                        ]),
+                    ]),
+                ));
+            }
+        }
+        for b in 0..w {
+            let bit = checked(u64::try_from(b).ok())?;
+            assignments.push(Assignment::new(
+                vec![format!("p_{}_{b}", self.sizes().len())],
+                BooleanExpr::constant(self.target().bit(bit)),
+            ));
+        }
+        let target =
+            CircuitSAT::try_new(Circuit::new(assignments)).map_err(Self::target_construction)?;
+        let source_var_indices = (0..self.sizes().len())
+            .map(|i| {
+                target
+                    .variable_names()
+                    .binary_search(&format!("x_{i}"))
+                    .map_err(|_| {
+                        ReductionError::invalid_target::<Self, CircuitSAT>(
+                            "missing selector variable",
+                        )
+                    })
+            })
+            .collect::<Result<_, _>>()?;
+        Ok(ReductionSubsetProductToCircuitSAT {
+            target,
+            source_var_indices,
+        })
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "subsetproduct_to_circuitsat",
+        build: || {
+            let source = SubsetProduct::new(vec![2u32, 3, 5], 30u32);
+            let reduction = ReduceTo::<CircuitSAT>::reduce_to(&source).expect("example reduction");
+            let target = reduction.target_problem();
+            let mut values = std::collections::HashMap::from([
+                ("x_0".to_owned(), true),
+                ("x_1".to_owned(), true),
+                ("x_2".to_owned(), true),
+            ]);
+            for assignment in &target.circuit().assignments {
+                let value = assignment.expr.evaluate(&values);
+                for name in &assignment.outputs {
+                    values.entry(name.clone()).or_insert(value);
+                }
+            }
+            let witness: Vec<bool> = target
+                .variable_names()
+                .iter()
+                .map(|name| values[name])
+                .collect();
+            assert_eq!(reduction.extract_solution(&witness).unwrap(), vec![true; 3]);
+            crate::example_db::specs::assemble_rule_example(
+                &source,
+                target,
+                vec![crate::export::SolutionPair {
+                    source_config: serde_json::json!([true, true, true]),
+                    target_config: serde_json::json!(witness),
+                }],
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/subsetproduct_circuitsat.rs"]
+mod tests;
diff --git a/src/unit_tests/models/formula/circuit.rs b/src/unit_tests/models/formula/circuit.rs
index 6fe536bc5..5f56564c0 100644
--- a/src/unit_tests/models/formula/circuit.rs
+++ b/src/unit_tests/models/formula/circuit.rs
@@ -308,3 +308,16 @@ fn test_circuit_sat_paper_example() {
     let all = solver.find_all_witnesses(&problem).unwrap();
     assert_eq!(all.len(), 4);
 }
+
+#[test]
+fn test_fallible_circuit_construction_preserves_simultaneous_equations() {
+    let circuit = Circuit::new(vec![
+        Assignment::new(vec!["z".into()], BooleanExpr::var("a")),
+        Assignment::new(vec!["z".into()], BooleanExpr::constant(true)),
+    ]);
+    let target = CircuitSAT::try_new(circuit).unwrap();
+    assert_eq!(target.variable_names(), &["a", "z"]);
+    assert!(target.is_valid_solution(&[true, true]).unwrap());
+    assert!(!target.is_valid_solution(&[false, true]).unwrap());
+    assert!(!target.is_valid_solution(&[false, false]).unwrap());
+}
diff --git a/src/unit_tests/rules/subsetproduct_circuitsat.rs b/src/unit_tests/rules/subsetproduct_circuitsat.rs
new file mode 100644
index 000000000..e15bbed55
--- /dev/null
+++ b/src/unit_tests/rules/subsetproduct_circuitsat.rs
@@ -0,0 +1,459 @@
+use crate::models::{CircuitSAT, SubsetProduct};
+use crate::rules::test_helpers::assert_satisfaction_round_trip_from_satisfaction_target;
+use crate::rules::{ReduceTo, ReductionResult};
+use crate::types::Or;
+use crate::{BruteForce, Problem};
+
+#[test]
+fn test_subsetproduct_to_circuitsat_closed_loop() {
+    for target in [1u32, 2] {
+        let source = SubsetProduct::new(Vec::<u32>::new(), target);
+        let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+        if target == 1 {
+            assert_satisfaction_round_trip_from_satisfaction_target(
+                &source,
+                &result,
+                "empty product",
+            );
+        }
+        assert_eq!(
+            BruteForce::new()
+                .solve(result.target_problem())
+                .unwrap()
+                .is_some(),
+            target == 1
+        );
+        assert_eq!(
+            crate::rules::AggregateReductionResult::extract_value(&result, Or(target == 1)),
+            Or(target == 1)
+        );
+    }
+}
+
+#[test]
+fn test_subsetproduct_to_circuitsat_paper_example() {
+    let source = SubsetProduct::new(vec![2u32, 3, 5], 30u32);
+    let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    let target = result.target_problem();
+    assert_eq!(target.num_variables(), 235);
+    assert_eq!(target.num_assignments(), 240);
+    assert_eq!(target.num_assignment_outputs(), 240);
+    assert_eq!(target.num_expression_nodes(), 808);
+}
+
+#[test]
+fn test_subsetproduct_to_circuitsat_rejects_invalid_witnesses() {
+    let source = SubsetProduct::new(vec![2u32], 2u32);
+    let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    let size = result.target_problem().num_variables();
+    for length in [0, size - 1, size, size + 1] {
+        assert!(result.extract_solution(&vec![false; length]).is_err());
+    }
+}
+
+use super::{dimensions, reserve_assignments, sum_sizes};
+use num_bigint::BigUint;
+use num_traits::{One, Zero};
+use std::collections::HashMap;
+
+// Evaluate fresh equations only; repeated outputs remain constraints in evaluate().
+fn extension(target: &CircuitSAT, selectors: &[bool]) -> (Vec<bool>, HashMap<String, bool>) {
+    let mut values: HashMap<_, _> = selectors
+        .iter()
+        .enumerate()
+        .map(|(i, &v)| (format!("x_{i}"), v))
+        .collect();
+    for assignment in &target.circuit().assignments {
+        for name in assignment.expr.variables() {
+            assert!(values.contains_key(&name), "undefined input {name}");
+        }
+        let value = assignment.expr.evaluate(&values);
+        for name in &assignment.outputs {
+            values.entry(name.clone()).or_insert(value);
+        }
+    }
+    let witness = target
+        .variable_names()
+        .iter()
+        .map(|name| values[name])
+        .collect();
+    (witness, values)
+}
+
+fn word(values: &HashMap<String, bool>, prefix: &str, suffix: &str, width: usize) -> BigUint {
+    (0..width)
+        .filter(|b| values[&format!("{prefix}{b}{suffix}")])
+        .fold(BigUint::zero(), |a, b| a + (BigUint::one() << b))
+}
+
+fn audit(source: &SubsetProduct) {
+    let result = ReduceTo::<CircuitSAT>::reduce_to(source).unwrap();
+    let target = result.target_problem();
+    let d: u64 = source.sizes().iter().map(BigUint::bits).sum();
+    let h: u64 = source.sizes().iter().map(BigUint::count_ones).sum();
+    let w = usize::try_from((d + 1).max(source.target().bits())).unwrap();
+    let n = source.sizes().len();
+    let h = usize::try_from(h).unwrap();
+    assert_eq!(target.num_variables(), w * (1 + n + 5 * h) + n);
+    assert_eq!(target.num_assignments(), w * (2 + n + 5 * h));
+    assert_eq!(target.num_assignment_outputs(), target.num_assignments());
+    assert_eq!(target.num_expression_nodes(), w * (2 + 8 * n + 15 * h));
+    for mask in 0..(1usize << n) {
+        let selectors: Vec<_> = (0..n).map(|i| mask & (1 << i) != 0).collect();
+        let (witness, values) = extension(target, &selectors);
+        let mut prefix = BigUint::one();
+        assert_eq!(word(&values, "p_0_", "", w), prefix);
+        for (i, a) in source.sizes().iter().enumerate() {
+            let mut sum = BigUint::zero();
+            for j in (0..a.bits()).filter(|&j| a.bit(j)) {
+                sum += &prefix << usize::try_from(j).unwrap();
+                assert!(sum.bits() <= u64::try_from(w).unwrap());
+                assert_eq!(word(&values, &format!("a_{}_{j}_", i + 1), "_s", w), sum);
+                assert!(!values[&format!("a_{}_{j}_{}_d", i + 1, w - 1)]);
+            }
+            assert_eq!(sum, &prefix * a);
+            if selectors[i] {
+                prefix *= a;
+            }
+            assert_eq!(word(&values, &format!("p_{}_", i + 1), "", w), prefix);
+        }
+        let feasible = prefix == *source.target();
+        assert_eq!(source.evaluate(&selectors).unwrap(), Or(feasible));
+        assert_eq!(target.evaluate(&witness).unwrap(), Or(feasible));
+        if feasible {
+            assert_eq!(result.extract_solution(&witness).unwrap(), selectors);
+        } else {
+            assert!(result.extract_solution(&witness).is_err());
+        }
+    }
+}
+
+#[test]
+fn test_carries_repeated_factors_and_unreachable_products() {
+    for (sizes, target) in [
+        (vec![3, 3], 9),
+        (vec![3, 3], 7),
+        (vec![3, 3], 15),
+        (vec![7, 7], 49),
+        (vec![2, 2, 2], 8),
+        (vec![2, 3, 5], 7),
+        (vec![1, 1, 1], 1),
+        (vec![1], 2),
+        (vec![5], 5),
+        (vec![5], 1),
+        (vec![5], 3),
+        (vec![2, 3], 4096),
+    ] {
+        audit(&SubsetProduct::new(sizes, target));
+    }
+}
+
+fn exhaustive(n: usize, masks: std::ops::Range<usize>) {
+    for mask in masks {
+        let sizes: Vec<u32> = (0..n)
+            .map(|i| if mask & (1 << i) == 0 { 1 } else { 2 })
+            .collect();
+        for target in 1..=(1u32 << n) + 1 {
+            audit(&SubsetProduct::new(sizes.clone(), target));
+        }
+    }
+}
+
+#[test]
+fn test_exhaustive_zero_to_two_elements() {
+    for n in 0..=2 {
+        exhaustive(n, 0..(1 << n));
+    }
+}
+#[test]
+fn test_exhaustive_three_elements() {
+    exhaustive(3, 0..8);
+}
+#[test]
+fn test_exhaustive_four_elements() {
+    exhaustive(4, 0..16);
+}
+// Partition the exhaustive five-element family to keep each test below five seconds.
+macro_rules! five_element_tests {
+    ($($name:ident: $start:expr),* $(,)?) => { $(
+        #[test] fn $name() { exhaustive(5, $start..$start + 4); }
+    )* };
+}
+five_element_tests! {
+    test_five_elements_ending_in_one_one_one: 0,
+    test_five_elements_ending_in_two_one_one: 4,
+    test_five_elements_ending_in_one_two_one: 8,
+    test_five_elements_ending_in_two_two_one: 12,
+    test_five_elements_ending_in_one_one_two: 16,
+    test_five_elements_ending_in_two_one_two: 20,
+    test_five_elements_ending_in_one_two_two: 24,
+    test_five_elements_ending_in_two_two_two: 28,
+}
+
+#[test]
+fn test_wide_selected_and_unselected_branches() {
+    let power = BigUint::one() << 128usize;
+    audit(&SubsetProduct::new(vec![power.clone()], power.clone()));
+    for target in [BigUint::from(3u32), (&power + 1u32) * 3u32] {
+        audit(&SubsetProduct::new(
+            vec![&power + 1u32, BigUint::from(3u32)],
+            target,
+        ));
+    }
+    audit(&SubsetProduct::new(
+        vec![BigUint::one()],
+        (BigUint::one() << 257usize) + 1u32,
+    ));
+}
+
+#[test]
+fn test_selector_names_preserve_numeric_order() {
+    let source = SubsetProduct::new(vec![1u32; 12], 1u32);
+    let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    let selectors: Vec<_> = (0..12).map(|i| [2, 9, 11].contains(&i)).collect();
+    let (witness, _) = extension(result.target_problem(), &selectors);
+    assert_eq!(result.extract_solution(&witness).unwrap(), selectors);
+}
+
+#[test]
+fn test_every_canonical_witness_bit_is_necessary() {
+    let source = SubsetProduct::new(vec![2u32, 3, 5], 30u32);
+    let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    let (mut witness, _) = extension(result.target_problem(), &[true; 3]);
+    assert_eq!(result.extract_solution(&witness).unwrap(), [true; 3]);
+    for i in 0..witness.len() {
+        witness[i] = !witness[i];
+        assert!(
+            result.extract_solution(&witness).is_err(),
+            "corruption at {i}"
+        );
+        witness[i] = !witness[i];
+    }
+}
+
+#[test]
+fn test_size_sums_and_all_dimension_boundaries() {
+    assert_eq!(sum_sizes([usize::MAX - 1, 1]).unwrap(), usize::MAX);
+    assert!(matches!(
+        sum_sizes([usize::MAX, 1]),
+        Err(crate::rules::ReductionError::IntegerOverflow { .. })
+    ));
+    let max = BigUint::from(usize::MAX);
+    let mut values = vec![0, 1, 2, usize::MAX];
+    for divisor in [
+        2,
+        5,
+        8,
+        15,
+        std::mem::size_of::<crate::models::formula::Assignment>(),
+        std::mem::size_of::<crate::models::formula::BooleanExpr>(),
+        std::mem::size_of::<String>(),
+    ] {
+        let boundary = usize::MAX / divisor;
+        values.extend([boundary - 1, boundary, boundary + 1]);
+    }
+    for &value in &values {
+        for (n, d, h, b) in [
+            (value, 0, 0, 1),
+            (0, value, 0, 1),
+            (0, 0, value, 1),
+            (0, 0, 0, value),
+            (1, 1, value, 2),
+        ] {
+            let width = (BigUint::from(d) + 1u32).max(BigUint::from(b));
+            let vars = &width * (BigUint::from(n) + BigUint::from(h) * 5u32 + 1u32) + n;
+            let assignments = &width * (BigUint::from(n) + BigUint::from(h) * 5u32 + 2u32);
+            let nodes = &width * (BigUint::from(n) * 8u32 + BigUint::from(h) * 15u32 + 2u32);
+            let fits = [&width, &vars, &assignments, &nodes]
+                .into_iter()
+                .all(|v| v <= &max)
+                && &assignments * std::mem::size_of::<crate::models::formula::Assignment>() <= max
+                && &nodes * std::mem::size_of::<crate::models::formula::BooleanExpr>() <= max
+                && &vars * std::mem::size_of::<String>() <= max;
+            match dimensions(n, d, h, b) {
+                Ok(actual) => {
+                    assert!(fits, "accepted overflowing dimensions {:?}", (n, d, h, b));
+                    assert_eq!(
+                        (
+                            BigUint::from(actual.0),
+                            BigUint::from(actual.1),
+                            BigUint::from(actual.2),
+                            BigUint::from(actual.3)
+                        ),
+                        (width, vars, assignments, nodes)
+                    );
+                }
+                Err(error) => {
+                    assert!(
+                        !fits,
+                        "rejected representable dimensions {:?}: {error}",
+                        (n, d, h, b)
+                    );
+                    assert!(matches!(
+                        error,
+                        crate::rules::ReductionError::IntegerOverflow { .. }
+                    ));
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn test_assignment_allocation_failure_is_typed() {
+    assert!(reserve_assignments(1).is_ok());
+    assert!(matches!(
+        reserve_assignments(usize::MAX),
+        Err(crate::rules::ReductionError::Construction { .. })
+    ));
+}
+
+#[test]
+fn test_missing_selector_coordinate_is_rejected() {
+    let source = SubsetProduct::new(vec![1u32], 1u32);
+    let mut result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    let (witness, _) = extension(result.target_problem(), &[true]);
+    result.source_var_indices[0] = usize::MAX;
+    assert!(result.extract_solution(&witness).is_err());
+}
+
+#[test]
+fn test_registered_decision_mapping_and_unavailable_magnitudes() {
+    let entries = crate::rules::registry::reduction_entries();
+    let entry = entries
+        .iter()
+        .find(|e| e.source_name == "SubsetProduct" && e.target_name == "CircuitSAT")
+        .unwrap();
+    let contract = entry.parameter_contract().unwrap();
+    assert!(contract.transform().is_none());
+    assert_eq!(contract.unavailable().len(), 4);
+    for field in [
+        "num_variables",
+        "num_assignments",
+        "num_assignment_outputs",
+        "num_expression_nodes",
+    ] {
+        assert!(
+            contract
+                .unavailable()
+                .iter()
+                .any(|entry| entry.field == field && !entry.reason.is_empty())
+        );
+    }
+    let source = SubsetProduct::new(vec![1u32], 1u32);
+    let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+    for value in [false, true] {
+        assert_eq!(
+            crate::rules::AggregateReductionResult::extract_value(&result, Or(value)),
+            Or(value)
+        );
+    }
+}
+
+#[test]
+fn test_all_target_witnesses_for_a_free_unit_selector() {
+    for target in [1u32, 2] {
+        let source = SubsetProduct::new(vec![1u32], target);
+        let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+        let witnesses = BruteForce::new()
+            .find_all_witnesses(result.target_problem())
+            .unwrap();
+        let mut extracted: Vec<_> = witnesses
+            .iter()
+            .map(|w| result.extract_solution(w).unwrap())
+            .collect();
+        extracted.sort();
+        let expected = if target == 1 {
+            vec![vec![false], vec![true]]
+        } else {
+            vec![]
+        };
+        assert_eq!(extracted, expected);
+        assert_eq!(
+            crate::rules::AggregateReductionResult::extract_value(
+                &result,
+                Or(!witnesses.is_empty())
+            ),
+            Or(target == 1)
+        );
+    }
+}
+
+#[cfg(feature = "example-db")]
+#[test]
+fn test_canonical_export_has_a_complete_valid_witness() {
+    let examples = super::canonical_rule_example_specs();
+    let example = (examples[0].build)();
+    let source: SubsetProduct = serde_json::from_value(example.source.instance).unwrap();
+    let target: CircuitSAT = serde_json::from_value(example.target.instance).unwrap();
+    let source_witness: Vec<bool> =
+        serde_json::from_value(example.solutions[0].source_config.clone()).unwrap();
+    let target_witness: Vec<bool> =
+        serde_json::from_value(example.solutions[0].target_config.clone()).unwrap();
+    assert_eq!(source.sizes(), &[2u32.into(), 3u32.into(), 5u32.into()]);
+    assert_eq!(source.target(), &30u32.into());
+    assert_eq!(source.evaluate(&source_witness).unwrap(), Or(true));
+    assert_eq!(target.evaluate(&target_witness).unwrap(), Or(true));
+    assert_eq!(
+        ReduceTo::<CircuitSAT>::reduce_to(&source)
+            .unwrap()
+            .extract_solution(&target_witness)
+            .unwrap(),
+        source_witness
+    );
+}
+
+fn random_extensions(n: usize) {
+    use rand::SeedableRng;
+    use rand::rngs::SmallRng;
+    use rand::seq::IndexedRandom;
+    let mut rng = SmallRng::seed_from_u64(42 + u64::try_from(n).unwrap());
+    for _ in 0..300 {
+        let sizes: Vec<u32> = (0..n)
+            .map(|_| *[1, 2, 3, 5, 7, 9, 15, 31].choose(&mut rng).unwrap())
+            .collect();
+        let selectors: Vec<bool> = (0..n)
+            .map(|_| *[false, true].choose(&mut rng).unwrap())
+            .collect();
+        let product: u64 = sizes
+            .iter()
+            .zip(&selectors)
+            .filter(|(_, selected)| **selected)
+            .map(|(&a, _)| u64::from(a))
+            .product();
+        let target = product + u64::from(*[false, true].choose(&mut rng).unwrap());
+        let source = SubsetProduct::new(sizes, target);
+        let result = ReduceTo::<CircuitSAT>::reduce_to(&source).unwrap();
+        let (witness, _) = extension(result.target_problem(), &selectors);
+        assert_eq!(
+            result.target_problem().evaluate(&witness).unwrap(),
+            Or(product == target)
+        );
+        if product == target {
+            assert_eq!(result.extract_solution(&witness).unwrap(), selectors);
+        } else {
+            assert!(result.extract_solution(&witness).is_err());
+        }
+    }
+}
+
+#[test]
+fn test_random_single_element_extensions() {
+    random_extensions(1);
+}
+#[test]
+fn test_random_two_element_extensions() {
+    random_extensions(2);
+}
+#[test]
+fn test_random_three_element_extensions() {
+    random_extensions(3);
+}
+#[test]
+fn test_random_four_element_extensions() {
+    random_extensions(4);
+}
+#[test]
+fn test_random_five_element_extensions() {
+    random_extensions(5);
+}