diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ
index d83a5fa0..fd3748ce 100644
--- a/docs/paper/reductions.typ
+++ b/docs/paper/reductions.typ
@@ -501,9 +501,9 @@ where $P$ is a penalty weight large enough that any constraint violation costs m
 ]
 
 #reduction-rule("KSatisfiability", "QUBO")[
-  Given a Max-2-SAT instance with $m$ clauses over $n$ variables, construct upper-triangular $Q in RR^(n times n)$ where each clause $(ell_i or ell_j)$ contributes a penalty gadget encoding its unique falsifying assignment.
+  Given a Max-$k$-SAT instance with $m$ clauses over $n$ variables, construct a QUBO that counts unsatisfied clauses. For $k = 2$, $Q in RR^(n times n)$ directly encodes quadratic penalties. For $k = 3$, Rosenberg quadratization introduces $m$ auxiliary variables, giving $Q in RR^((n+m) times (n+m))$.
 ][
-  _Construction._ Applying the penalty method (@sec:penalty-method), each 2-literal clause has exactly one falsifying assignment (both literals false). The penalty for that assignment is a quadratic function of $x_i, x_j$:
+  *Case $k = 2$.* Applying the penalty method (@sec:penalty-method), each 2-literal clause has exactly one falsifying assignment (both literals false). The penalty for that assignment is a quadratic function of $x_i, x_j$:
 
   #table(
     columns: (auto, auto, auto, auto),
@@ -517,6 +517,12 @@ where $P$ is a penalty weight large enough that any constraint violation costs m
   )
 
   Summing over all clauses, $f(bold(x)) = sum_j "penalty"_j (bold(x))$ counts falsified clauses. Minimizers of $f$ maximize satisfied clauses.
+
+  *Case $k = 3$ (Rosenberg quadratization).* For each clause $(ell_1 or ell_2 or ell_3)$, define complement variables $y_i = overline(ell_i)$ (so $y_i = x_i$ if the literal is negated, $y_i = 1 - x_i$ if positive). The clause is violated when $y_1 y_2 y_3 = 1$. This cubic penalty is reduced to quadratic form by introducing an auxiliary variable $a$ and the substitution $a = y_1 y_2$, enforced via a Rosenberg penalty with weight $M$:
+  $ H = a dot y_3 + M (y_1 y_2 - 2 y_1 a - 2 y_2 a + 3a) $
+  where $M = 2$ suffices. For any binary assignment, if $a = y_1 y_2$ the penalty term vanishes and $H = y_1 y_2 y_3$ counts the clause violation. If $a != y_1 y_2$, the penalty $M(dots.c) >= 1$ makes this suboptimal.
+
+  Each clause adds one auxiliary variable (indices $n, n+1, ..., n+m-1$), so the total QUBO has $n + m$ variables. Solution extraction discards auxiliary variables: return $bold(x)[0..n]$.
 ]
 
 #reduction-rule("ILP", "QUBO")[
diff --git a/docs/plans/2026-02-10-improve-example-instances-impl.md b/docs/plans/2026-02-10-improve-example-instances-impl.md
index 5ae27e55..96ec6291 100644
--- a/docs/plans/2026-02-10-improve-example-instances-impl.md
+++ b/docs/plans/2026-02-10-improve-example-instances-impl.md
@@ -1,4 +1,4 @@
-# Improve Example Instances — Implementation Plan
+# Improve Example Instances — Implementation Plan (v2)
 
 > **For Claude:** REQUIRED SUB-SKILL: Use superpowers:executing-plans to implement this plan task-by-task.
 
diff --git a/examples/reduction_circuitsat_to_spinglass.rs b/examples/reduction_circuitsat_to_spinglass.rs
index e9b0a019..e9f1770e 100644
--- a/examples/reduction_circuitsat_to_spinglass.rs
+++ b/examples/reduction_circuitsat_to_spinglass.rs
@@ -6,8 +6,11 @@
 //! of gadget Hamiltonians; ground states correspond to satisfying assignments.
 //!
 //! ## This Example
-//! - Instance: Simple AND gate circuit (c = x AND y)
-//! - Source: CircuitSAT with 2 inputs
+//! - Instance: 1-bit full adder circuit (a, b, cin -> sum, cout)
+//!   - sum = a XOR b XOR cin (via intermediate t = a XOR b)
+//!   - cout = (a AND b) OR (cin AND t)
+//!   - 5 gates (2 XOR, 2 AND, 1 OR), ~8 variables
+//! - Source: CircuitSAT with 3 inputs
 //! - Target: SpinGlass
 //!
 //! ## Output
@@ -19,18 +22,40 @@ use problemreductions::prelude::*;
 use problemreductions::topology::{Graph, SimpleGraph};
 
 fn main() {
-    // 1. Create CircuitSAT instance: c = x AND y
-    //    This is a simple circuit with one AND gate.
+    // 1. Create CircuitSAT instance: 1-bit full adder
+    //    sum = a XOR b XOR cin, cout = (a AND b) OR (cin AND (a XOR b))
+    //    Decomposed into 5 gates with intermediate variables t, ab, cin_t.
     let circuit = Circuit::new(vec![
+        // Intermediate: t = a XOR b
         Assignment::new(
-            vec!["c".to_string()],
-            BooleanExpr::and(vec![BooleanExpr::var("x"), BooleanExpr::var("y")]),
+            vec!["t".to_string()],
+            BooleanExpr::xor(vec![BooleanExpr::var("a"), BooleanExpr::var("b")]),
+        ),
+        // sum = t XOR cin
+        Assignment::new(
+            vec!["sum".to_string()],
+            BooleanExpr::xor(vec![BooleanExpr::var("t"), BooleanExpr::var("cin")]),
+        ),
+        // ab = a AND b
+        Assignment::new(
+            vec!["ab".to_string()],
+            BooleanExpr::and(vec![BooleanExpr::var("a"), BooleanExpr::var("b")]),
+        ),
+        // cin_t = cin AND t
+        Assignment::new(
+            vec!["cin_t".to_string()],
+            BooleanExpr::and(vec![BooleanExpr::var("cin"), BooleanExpr::var("t")]),
+        ),
+        // cout = ab OR cin_t
+        Assignment::new(
+            vec!["cout".to_string()],
+            BooleanExpr::or(vec![BooleanExpr::var("ab"), BooleanExpr::var("cin_t")]),
         ),
     ]);
     let circuit_sat = CircuitSAT::<i32>::new(circuit);
 
     println!("=== Circuit-SAT to Spin Glass Reduction ===\n");
-    println!("Source circuit: c = x AND y");
+    println!("Source circuit: 1-bit full adder (a, b, cin -> sum, cout)");
     println!(
         "  {} variables: {:?}",
         circuit_sat.num_variables(),
@@ -51,9 +76,9 @@ fn main() {
         sg.num_spins(),
         sg.graph().num_edges()
     );
-    println!("  Each logic gate becomes a spin glass gadget.");
-    println!("  AND gadget: 3 spins with J=[1,-2,-2], h=[-1,-1,2]");
-    println!("  Ground states encode valid truth table entries.");
+    println!("  Each logic gate (AND, OR, XOR) becomes a spin glass gadget.");
+    println!("  Gadget ground states encode valid truth table entries.");
+    println!("  Full adder uses 5 gadgets for its 5 gate decomposition.");
 
     // 3. Solve the target SpinGlass problem
     let solver = BruteForce::new();
diff --git a/examples/reduction_factoring_to_circuitsat.rs b/examples/reduction_factoring_to_circuitsat.rs
index a29fdcac..455ff92b 100644
--- a/examples/reduction_factoring_to_circuitsat.rs
+++ b/examples/reduction_factoring_to_circuitsat.rs
@@ -5,12 +5,12 @@
 //! iff N can be factored within the given bit bounds.
 //!
 //! ## This Example
-//! - Instance: Factor 6 = 2 * 3 (m=2 bits, n=2 bits)
+//! - Instance: Factor 35 = 5 × 7 (m=3 bits, n=3 bits)
 //! - Reference: Based on ProblemReductions.jl factoring example
-//! - Source: Factoring(2, 2, 6)
+//! - Source: Factoring(3, 3, 35)
 //! - Target: CircuitSAT
 //!
-//! We solve the source Factoring problem directly with BruteForce (only 4 binary
+//! We solve the source Factoring problem directly with BruteForce (only 6 binary
 //! variables), then verify the reduction produces a valid CircuitSAT encoding by
 //! simulating the circuit forward from a known factorization to build a complete
 //! satisfying assignment.
@@ -40,9 +40,9 @@ fn simulate_circuit(
 }
 
 fn main() {
-    // 1. Create Factoring instance: factor 6 with 2-bit factors
-    //    Possible: 2*3=6 or 3*2=6
-    let factoring = Factoring::new(2, 2, 6);
+    // 1. Create Factoring instance: factor 35 with 3-bit factors
+    //    Possible: 5*7=35 or 7*5=35
+    let factoring = Factoring::new(3, 3, 35);
 
     println!("=== Factoring to Circuit-SAT Reduction ===\n");
     println!(
@@ -58,7 +58,7 @@ fn main() {
         factoring.n()
     );
 
-    // 2. Solve the source Factoring problem directly (only 4 binary variables)
+    // 2. Solve the source Factoring problem directly (only 6 binary variables)
     let solver = BruteForce::new();
     let factoring_solutions = solver.find_best(&factoring);
     println!("\nFactoring solutions found: {}", factoring_solutions.len());
@@ -93,7 +93,7 @@ fn main() {
         a, b, a * b, factoring_sol
     );
 
-    // Set input variables: p1, p2 for first factor, q1, q2 for second factor
+    // Set input variables: p1..p3 for first factor, q1..q3 for second factor
     let mut input_values: HashMap<String, bool> = HashMap::new();
     for (i, &bit) in factoring_sol.iter().enumerate().take(factoring.m()) {
         input_values.insert(format!("p{}", i + 1), bit == 1);
@@ -161,7 +161,7 @@ fn main() {
         });
     }
 
-    println!("\nReduction verified successfully: {} = {} * {}", factoring.target(), a, b);
+    println!("\nReduction verified successfully: 35 = 5 * 7");
 
     // 6. Export JSON
     let overhead = lookup_overhead("Factoring", "CircuitSAT")
diff --git a/examples/reduction_factoring_to_ilp.rs b/examples/reduction_factoring_to_ilp.rs
index 152eb8e9..62ca1567 100644
--- a/examples/reduction_factoring_to_ilp.rs
+++ b/examples/reduction_factoring_to_ilp.rs
@@ -10,9 +10,9 @@
 //! Objective: feasibility (minimize 0).
 //!
 //! ## This Example
-//! - Instance: Factor 15 = 3 * 5 (m=4 bits, n=4 bits)
+//! - Instance: Factor 35 = 5 × 7 (m=3 bits, n=3 bits)
 //! - NOTE: Uses ILPSolver (not BruteForce) since the ILP has many variables
-//! - Target ILP: 4+4+16+8 = 32 variables
+//! - Target ILP: ~21 variables (factor bits + product bits + carries)
 //!
 //! ## Output
 //! Exports `docs/paper/examples/factoring_to_ilp.json` for use in paper code blocks.
@@ -22,8 +22,8 @@ use problemreductions::prelude::*;
 use problemreductions::solvers::ILPSolver;
 
 fn main() {
-    // 1. Create Factoring instance: find p (4-bit) x q (4-bit) = 15
-    let problem = Factoring::new(4, 4, 15);
+    // 1. Create Factoring instance: find p (3-bit) x q (3-bit) = 35
+    let problem = Factoring::new(3, 3, 35);
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&problem);
@@ -36,18 +36,18 @@ fn main() {
 
     // 4. Solve ILP using ILPSolver (too many variables for BruteForce)
     let solver = ILPSolver::new();
-    let ilp_solution = solver.solve(ilp).expect("ILP should be feasible for 15 = 3 * 5");
+    let ilp_solution = solver.solve(ilp).expect("ILP should be feasible for 35 = 5 * 7");
     println!("\n=== Solution ===");
-    println!("ILP solution found (first 8 vars): {:?}", &ilp_solution[..8]);
+    println!("ILP solution found (first 6 vars): {:?}", &ilp_solution[..6]);
 
     // 5. Extract factoring solution
     let extracted = reduction.extract_solution(&ilp_solution);
     println!("Source Factoring solution: {:?}", extracted);
 
-    // 6. Verify: read factors and confirm p * q = 15
+    // 6. Verify: read factors and confirm p * q = 35
     let (p, q) = problem.read_factors(&extracted);
     println!("Factors: {} x {} = {}", p, q, p * q);
-    assert_eq!(p * q, 15);
+    assert_eq!(p * q, 35);
     println!("\nReduction verified successfully");
 
     // 7. Collect solutions and export JSON
diff --git a/examples/reduction_ilp_to_qubo.rs b/examples/reduction_ilp_to_qubo.rs
index c4ec98ae..d3acb820 100644
--- a/examples/reduction_ilp_to_qubo.rs
+++ b/examples/reduction_ilp_to_qubo.rs
@@ -17,12 +17,14 @@
 //! solutions with better objective values.
 //!
 //! ## This Example
-//! - Instance: 3 binary projects with values [1, 2, 3]
-//!   - Constraint 1: x0 + x1 <= 1 (projects 0 and 1 share a team)
-//!   - Constraint 2: x1 + x2 <= 1 (projects 1 and 2 share equipment)
-//!   - Objective: maximize 1*x0 + 2*x1 + 3*x2
-//! - Expected: Select projects {0, 2} for total value 4 (x0 and x2 don't
-//!   share resources)
+//! - Instance: 6-variable binary knapsack problem
+//!   - Items with weights [3, 2, 5, 4, 2, 3] and values [10, 7, 12, 8, 6, 9]
+//!   - Constraint 1: 3x0 + 2x1 + 5x2 + 4x3 + 2x4 + 3x5 <= 10 (weight capacity)
+//!   - Constraint 2: x0 + x1 + x2 <= 2 (category A limit)
+//!   - Constraint 3: x3 + x4 + x5 <= 2 (category B limit)
+//!   - Objective: maximize 10x0 + 7x1 + 12x2 + 8x3 + 6x4 + 9x5
+//! - Expected: Select items that maximize total value while satisfying all
+//!   weight and category constraints
 //!
 //! ## Outputs
 //! - `docs/paper/examples/ilp_to_qubo.json` — reduction structure
@@ -39,35 +41,51 @@ use problemreductions::prelude::*;
 fn main() {
     println!("=== ILP (Binary) -> QUBO Reduction ===\n");
 
-    // 3 projects with values 1, 2, 3
-    // Constraint 1: x0 + x1 <= 1 (shared team)
-    // Constraint 2: x1 + x2 <= 1 (shared equipment)
+    // 6-variable binary knapsack problem
+    // Items with weights [3, 2, 5, 4, 2, 3] and values [10, 7, 12, 8, 6, 9]
+    // Constraint 1: knapsack weight capacity <= 10
+    // Constraint 2: category A items (x0, x1, x2) limited to 2
+    // Constraint 3: category B items (x3, x4, x5) limited to 2
     let ilp = ILP::binary(
-        3,
+        6,
         vec![
-            LinearConstraint::le(vec![(0, 1.0), (1, 1.0)], 1.0),
-            LinearConstraint::le(vec![(1, 1.0), (2, 1.0)], 1.0),
+            // Knapsack weight constraint: 3x0 + 2x1 + 5x2 + 4x3 + 2x4 + 3x5 <= 10
+            LinearConstraint::le(
+                vec![(0, 3.0), (1, 2.0), (2, 5.0), (3, 4.0), (4, 2.0), (5, 3.0)],
+                10.0,
+            ),
+            // Category A limit: x0 + x1 + x2 <= 2
+            LinearConstraint::le(vec![(0, 1.0), (1, 1.0), (2, 1.0)], 2.0),
+            // Category B limit: x3 + x4 + x5 <= 2
+            LinearConstraint::le(vec![(3, 1.0), (4, 1.0), (5, 1.0)], 2.0),
         ],
-        vec![(0, 1.0), (1, 2.0), (2, 3.0)],
+        vec![(0, 10.0), (1, 7.0), (2, 12.0), (3, 8.0), (4, 6.0), (5, 9.0)],
         ObjectiveSense::Maximize,
     );
 
-    let project_names = ["Alpha", "Beta", "Gamma"];
+    let item_names = ["Item0", "Item1", "Item2", "Item3", "Item4", "Item5"];
+    let weights = [3, 2, 5, 4, 2, 3];
+    let values = [10, 7, 12, 8, 6, 9];
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&ilp);
     let qubo = reduction.target_problem();
 
-    println!("Source: ILP (binary) with 3 variables, 2 constraints");
-    println!("  Objective: maximize 1*x0 + 2*x1 + 3*x2");
-    println!("  Constraint 1: x0 + x1 <= 1 (shared team)");
-    println!("  Constraint 2: x1 + x2 <= 1 (shared equipment)");
+    println!("Source: ILP (binary) with 6 variables, 3 constraints");
+    println!("  Objective: maximize 10x0 + 7x1 + 12x2 + 8x3 + 6x4 + 9x5");
+    println!("  Constraint 1: 3x0 + 2x1 + 5x2 + 4x3 + 2x4 + 3x5 <= 10 (weight capacity)");
+    println!("  Constraint 2: x0 + x1 + x2 <= 2 (category A limit)");
+    println!("  Constraint 3: x3 + x4 + x5 <= 2 (category B limit)");
     println!("Target: QUBO with {} variables", qubo.num_variables());
-    println!("Q matrix ({}x{}):", qubo.matrix().len(), qubo.matrix().len());
-    for row in qubo.matrix() {
-        let formatted: Vec<String> = row.iter().map(|v| format!("{:7.1}", v)).collect();
-        println!("  [{}]", formatted.join(", "));
-    }
+    println!(
+        "  (6 original + {} slack variables for inequality constraints)",
+        qubo.num_variables() - 6
+    );
+    println!(
+        "Q matrix size: {}x{}",
+        qubo.matrix().len(),
+        qubo.matrix().len()
+    );
 
     // Solve QUBO with brute force
     let solver = BruteForce::new();
@@ -82,17 +100,31 @@ fn main() {
             .iter()
             .enumerate()
             .filter(|(_, &x)| x == 1)
-            .map(|(i, _)| project_names[i].to_string())
+            .map(|(i, _)| item_names[i].to_string())
             .collect();
-        let value = ilp.solution_size(&extracted).size;
+        let total_weight: i32 = extracted
+            .iter()
+            .enumerate()
+            .filter(|(_, &x)| x == 1)
+            .map(|(i, _)| weights[i])
+            .sum();
+        let total_value: i32 = extracted
+            .iter()
+            .enumerate()
+            .filter(|(_, &x)| x == 1)
+            .map(|(i, _)| values[i])
+            .sum();
         println!(
-            "  Selected projects: {:?} (total value: {:.0})",
-            selected, value
+            "  Selected items: {:?} (total weight: {}, total value: {})",
+            selected, total_weight, total_value
         );
 
         // Closed-loop verification: check solution is valid in original problem
         let sol_size = ilp.solution_size(&extracted);
-        assert!(sol_size.is_valid, "Solution must be valid in source problem");
+        assert!(
+            sol_size.is_valid,
+            "Solution must be valid in source problem"
+        );
 
         solutions.push(SolutionPair {
             source_config: extracted,
@@ -103,8 +135,7 @@ fn main() {
     println!("\nVerification passed: all solutions are feasible and optimal");
 
     // Export JSON
-    let overhead = lookup_overhead("ILP", "QUBO")
-        .expect("ILP -> QUBO overhead not found");
+    let overhead = lookup_overhead("ILP", "QUBO").expect("ILP -> QUBO overhead not found");
 
     let data = ReductionData {
         source: ProblemSide {
diff --git a/examples/reduction_kcoloring_to_ilp.rs b/examples/reduction_kcoloring_to_ilp.rs
index 9c0318b8..f36ed520 100644
--- a/examples/reduction_kcoloring_to_ilp.rs
+++ b/examples/reduction_kcoloring_to_ilp.rs
@@ -8,21 +8,23 @@
 //! Objective: feasibility (minimize 0).
 //!
 //! ## This Example
-//! - Instance: Triangle K3 (3 vertices, 3 edges) with 3 colors
-//! - Source KColoring: feasible, each vertex gets a distinct color
-//! - Target ILP: 9 binary variables (3 vertices * 3 colors), 12 constraints
+//! - Instance: Petersen graph (10 vertices, 15 edges) with 3 colors, χ=3
+//! - Source KColoring: feasible, each vertex gets a color such that no adjacent vertices share a color
+//! - Target ILP: 30 binary variables (10 vertices * 3 colors), many constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/kcoloring_to_ilp.json` and `kcoloring_to_ilp.result.json`.
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
-use problemreductions::solvers::BruteForceFloat;
+use problemreductions::solvers::ILPSolver;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create KColoring instance: triangle K3 with 3 colors
-    let coloring = KColoring::<3, SimpleGraph, i32>::new(3, vec![(0, 1), (1, 2), (0, 2)]);
+    // 1. Create KColoring instance: Petersen graph (10 vertices, 15 edges) with 3 colors, χ=3
+    let (num_vertices, edges) = petersen();
+    let coloring = KColoring::<3, SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&coloring);
@@ -33,17 +35,14 @@ fn main() {
     println!("Source: KColoring<3> with {} variables", coloring.num_variables());
     println!("Target: ILP with {} variables, {} constraints", ilp.num_vars, ilp.constraints.len());
 
-    // 4. Solve target ILP
-    let solver = BruteForce::new();
-    let ilp_solutions = solver.find_best_float(ilp);
+    // 4. Solve target ILP using HiGHS solver (BruteForce on 30 vars is too slow)
+    let solver = ILPSolver::new();
+    let ilp_solution = solver.solve(ilp).expect("ILP should be feasible");
     println!("\n=== Solution ===");
-    println!("ILP solutions found: {}", ilp_solutions.len());
-
-    let ilp_solution = &ilp_solutions[0].0;
     println!("ILP solution: {:?}", ilp_solution);
 
     // 5. Extract source solution
-    let coloring_solution = reduction.extract_solution(ilp_solution);
+    let coloring_solution = reduction.extract_solution(&ilp_solution);
     println!("Source Coloring solution: {:?}", coloring_solution);
 
     // 6. Verify
@@ -54,15 +53,13 @@ fn main() {
 
     // 7. Collect solutions and export JSON
     let mut solutions = Vec::new();
-    for (target_config, _score) in &ilp_solutions {
-        let source_sol = reduction.extract_solution(target_config);
-        let s = coloring.solution_size(&source_sol);
-        assert!(s.is_valid);
-        solutions.push(SolutionPair {
-            source_config: source_sol,
-            target_config: target_config.clone(),
-        });
-    }
+    let source_sol = reduction.extract_solution(&ilp_solution);
+    let s = coloring.solution_size(&source_sol);
+    assert!(s.is_valid);
+    solutions.push(SolutionPair {
+        source_config: source_sol,
+        target_config: ilp_solution,
+    });
 
     let overhead = lookup_overhead("KColoring", "ILP")
         .expect("KColoring -> ILP overhead not found");
diff --git a/examples/reduction_kcoloring_to_qubo.rs b/examples/reduction_kcoloring_to_qubo.rs
index 5de9177a..4600fdbf 100644
--- a/examples/reduction_kcoloring_to_qubo.rs
+++ b/examples/reduction_kcoloring_to_qubo.rs
@@ -14,10 +14,10 @@
 //! The QUBO has n*K variables (n vertices, K colors).
 //!
 //! ## This Example
-//! - Instance: Complete graph K3 (triangle) with 3 colors
-//! - Source: KColoring<3> on 3 vertices, 3 edges
-//! - QUBO variables: 9 (3 vertices x 3 colors, one-hot encoding)
-//! - Expected: 6 valid 3-colorings (3! = 6 permutations of 3 colors on 3 vertices)
+//! - Instance: House graph (5 vertices, 6 edges) with 3 colors, χ=3
+//! - Source: KColoring<3> on 5 vertices, 6 edges
+//! - QUBO variables: 15 (5 vertices x 3 colors, one-hot encoding)
+//! - BruteForce on 15 variables (2^15 = 32768) completes quickly
 //!
 //! ## Outputs
 //! - `docs/paper/examples/coloring_to_qubo.json` — reduction structure
@@ -30,23 +30,24 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::house;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("=== K-Coloring -> QUBO Reduction ===\n");
 
-    // Triangle K3: all 3 vertices are adjacent
-    let edges = vec![(0, 1), (1, 2), (0, 2)];
-    let kc = KColoring::<3, SimpleGraph, i32>::new(3, edges.clone());
+    // House graph: 5 vertices, 6 edges (square base + triangle roof), χ=3
+    let (num_vertices, edges) = house();
+    let kc = KColoring::<3, SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&kc);
     let qubo = reduction.target_problem();
 
     let colors = ["Red", "Green", "Blue"];
-    println!("Source: KColoring<3> on triangle K3 (3 vertices, 3 edges)");
+    println!("Source: KColoring<3> on house graph (5 vertices, 6 edges)");
     println!(
-        "Target: QUBO with {} variables (one-hot: 3 vertices x 3 colors)",
+        "Target: QUBO with {} variables (one-hot: 5 vertices x 3 colors)",
         qubo.num_variables()
     );
     println!("Q matrix:");
@@ -81,13 +82,10 @@ fn main() {
         });
     }
 
-    // K3 with 3 colors has exactly 3! = 6 valid colorings
-    assert_eq!(
-        qubo_solutions.len(),
-        6,
-        "Triangle K3 with 3 colors should have exactly 6 valid colorings"
+    println!(
+        "\nVerification passed: {} valid colorings found",
+        qubo_solutions.len()
     );
-    println!("\nVerification passed: 6 valid colorings found");
 
     // Export JSON
     let overhead = lookup_overhead("KColoring", "QUBO")
diff --git a/examples/reduction_ksatisfiability_to_qubo.rs b/examples/reduction_ksatisfiability_to_qubo.rs
index 714715c9..2c883df9 100644
--- a/examples/reduction_ksatisfiability_to_qubo.rs
+++ b/examples/reduction_ksatisfiability_to_qubo.rs
@@ -1,4 +1,4 @@
-//! # K-Satisfiability (2-SAT) to QUBO Reduction (Penalty Method)
+//! # K-Satisfiability (3-SAT) to QUBO Reduction (Penalty Method)
 //!
 //! ## Mathematical Relationship
 //! The Maximum K-Satisfiability problem maps a CNF formula with k-literal clauses
@@ -11,14 +11,20 @@
 //! The total QUBO Hamiltonian H = -sum_j H_j is minimized when the maximum
 //! number of clauses is satisfied.
 //!
+//! For 3-SAT clauses, the cubic penalty terms are quadratized using
+//! Rosenberg's substitution, introducing one auxiliary variable per clause.
+//!
 //! ## This Example
-//! - Instance: 2-SAT with 3 variables and 4 clauses
-//!   - C1: x1 OR x2
-//!   - C2: NOT x1 OR x3
-//!   - C3: x2 OR NOT x3
-//!   - C4: NOT x2 OR NOT x3
-//! - QUBO variables: 3 (one per Boolean variable)
-//! - Expected: Assignments satisfying all 4 clauses (if possible) or
+//! - Instance: 3-SAT with 5 variables and 7 clauses
+//!   - C1: x1 OR x2 OR NOT x3
+//!   - C2: NOT x1 OR x3 OR x4
+//!   - C3: x2 OR NOT x4 OR x5
+//!   - C4: NOT x2 OR x3 OR NOT x5
+//!   - C5: x1 OR NOT x3 OR x5
+//!   - C6: NOT x1 OR NOT x2 OR x4
+//!   - C7: x3 OR NOT x4 OR NOT x5
+//! - QUBO variables: 5 original + 7 auxiliary = 12 total
+//! - Expected: Assignments satisfying all 7 clauses (if possible) or
 //!   maximizing satisfied clauses
 //!
 //! ## Outputs
@@ -34,29 +40,35 @@ use problemreductions::export::*;
 use problemreductions::prelude::*;
 
 fn main() {
-    println!("=== K-Satisfiability (2-SAT) -> QUBO Reduction ===\n");
+    println!("=== K-Satisfiability (3-SAT) -> QUBO Reduction ===\n");
 
-    // 4 clauses over 3 variables
+    // 7 clauses over 5 variables
     let clauses = vec![
-        CNFClause::new(vec![1, 2]),   // x1 OR x2
-        CNFClause::new(vec![-1, 3]),  // NOT x1 OR x3
-        CNFClause::new(vec![2, -3]),  // x2 OR NOT x3
-        CNFClause::new(vec![-2, -3]), // NOT x2 OR NOT x3
+        CNFClause::new(vec![1, 2, -3]),   // x1 OR x2 OR NOT x3
+        CNFClause::new(vec![-1, 3, 4]),   // NOT x1 OR x3 OR x4
+        CNFClause::new(vec![2, -4, 5]),   // x2 OR NOT x4 OR x5
+        CNFClause::new(vec![-2, 3, -5]),  // NOT x2 OR x3 OR NOT x5
+        CNFClause::new(vec![1, -3, 5]),   // x1 OR NOT x3 OR x5
+        CNFClause::new(vec![-1, -2, 4]),  // NOT x1 OR NOT x2 OR x4
+        CNFClause::new(vec![3, -4, -5]),  // x3 OR NOT x4 OR NOT x5
     ];
     let clause_strings = [
-        "x1 OR x2".to_string(),
-        "NOT x1 OR x3".to_string(),
-        "x2 OR NOT x3".to_string(),
-        "NOT x2 OR NOT x3".to_string(),
+        "x1 OR x2 OR NOT x3".to_string(),
+        "NOT x1 OR x3 OR x4".to_string(),
+        "x2 OR NOT x4 OR x5".to_string(),
+        "NOT x2 OR x3 OR NOT x5".to_string(),
+        "x1 OR NOT x3 OR x5".to_string(),
+        "NOT x1 OR NOT x2 OR x4".to_string(),
+        "x3 OR NOT x4 OR NOT x5".to_string(),
     ];
 
-    let ksat = KSatisfiability::<2, i32>::new(3, clauses);
+    let ksat = KSatisfiability::<3, i32>::new(5, clauses);
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&ksat);
     let qubo = reduction.target_problem();
 
-    println!("Source: KSatisfiability<2> with 3 variables, 4 clauses");
+    println!("Source: KSatisfiability<3> with 5 variables, 7 clauses");
     for (i, c) in clause_strings.iter().enumerate() {
         println!("  C{}: {}", i + 1, c);
     }
@@ -102,12 +114,12 @@ fn main() {
 
     let data = ReductionData {
         source: ProblemSide {
-            problem: KSatisfiability::<2, i32>::NAME.to_string(),
-            variant: variant_to_map(KSatisfiability::<2, i32>::variant()),
+            problem: KSatisfiability::<3, i32>::NAME.to_string(),
+            variant: variant_to_map(KSatisfiability::<3, i32>::variant()),
             instance: serde_json::json!({
                 "num_vars": ksat.num_vars(),
                 "num_clauses": ksat.clauses().len(),
-                "k": 2,
+                "k": 3,
             }),
         },
         target: ProblemSide {
diff --git a/examples/reduction_maxcut_to_spinglass.rs b/examples/reduction_maxcut_to_spinglass.rs
index 44af6a6a..b3c4433b 100644
--- a/examples/reduction_maxcut_to_spinglass.rs
+++ b/examples/reduction_maxcut_to_spinglass.rs
@@ -6,9 +6,9 @@
 //! energy -sum J_{ij} s_i s_j since s_i s_j = -1 when vertices are on opposite sides.
 //!
 //! ## This Example
-//! - Instance: Triangle K3 with unit edge weights
-//! - Source MaxCut: 3 vertices, 3 edges, max cut = 2
-//! - Target SpinGlass: 3 spins
+//! - Instance: Petersen graph (10 vertices, 15 edges) with unit edge weights
+//! - Source MaxCut: 10 vertices, 15 edges
+//! - Target SpinGlass: 10 spins
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maxcut_to_spinglass.json` and `maxcut_to_spinglass.result.json`.
@@ -17,10 +17,12 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    let maxcut = MaxCut::<SimpleGraph, i32>::new(3, vec![(0, 1, 1), (1, 2, 1), (0, 2, 1)]);
+    let (num_vertices, edges) = petersen();
+    let maxcut = MaxCut::<SimpleGraph, i32>::unweighted(num_vertices, edges.clone());
 
     let reduction = ReduceTo::<SpinGlass<SimpleGraph, i32>>::reduce_to(&maxcut);
     let sg = reduction.target_problem();
diff --git a/examples/reduction_maximumclique_to_ilp.rs b/examples/reduction_maximumclique_to_ilp.rs
index 4cbd929e..e75700b8 100644
--- a/examples/reduction_maximumclique_to_ilp.rs
+++ b/examples/reduction_maximumclique_to_ilp.rs
@@ -6,11 +6,10 @@
 //! Objective: maximize sum of w_v * x_v.
 //!
 //! ## This Example
-//! - Instance: 4-vertex graph with a triangle subgraph on {0,1,2} plus vertex 3
-//!   connected only to vertex 2. Edges: 0-1, 0-2, 1-2, 2-3.
-//! - Source MaximumClique: max clique is {0,1,2} (size 3)
-//! - Target ILP: 4 binary variables, 3 non-edge constraints
-//!   (non-edges: (0,3), (1,3))
+//! - Instance: Octahedron graph (K_{2,2,2}) with 6 vertices and 12 edges.
+//! - Source MaximumClique: max clique is size 3
+//! - Target ILP: 6 binary variables, 3 non-edge constraints
+//!   (non-edges: opposite vertex pairs (0,5), (1,4), (2,3))
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumclique_to_ilp.json` and `maximumclique_to_ilp.result.json`.
@@ -18,11 +17,13 @@
 use problemreductions::export::*;
 use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
+use problemreductions::topology::small_graphs::octahedral;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create MaximumClique instance: 4 vertices, triangle {0,1,2} plus vertex 3 connected to 2
-    let clique = MaximumClique::<SimpleGraph, i32>::new(4, vec![(0, 1), (0, 2), (1, 2), (2, 3)]);
+    // 1. Create MaximumClique instance: Octahedron (K_{2,2,2}), 6 vertices, 12 edges, clique number 3
+    let (num_vertices, edges) = octahedral();
+    let clique = MaximumClique::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&clique);
diff --git a/examples/reduction_maximumindependentset_to_ilp.rs b/examples/reduction_maximumindependentset_to_ilp.rs
index c6eed4ac..1d461944 100644
--- a/examples/reduction_maximumindependentset_to_ilp.rs
+++ b/examples/reduction_maximumindependentset_to_ilp.rs
@@ -6,9 +6,9 @@
 //! Objective: maximize sum of w_v * x_v.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges: 0-1-2-3)
-//! - Source IS: max size 2 (e.g., {0,2}, {0,3}, {1,3})
-//! - Target ILP: 4 binary variables, 3 constraints
+//! - Instance: Petersen graph (10 vertices, 15 edges, 3-regular)
+//! - Source IS: max size 4
+//! - Target ILP: 10 binary variables, 15 constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumindependentset_to_ilp.json` and `maximumindependentset_to_ilp.result.json`.
@@ -16,12 +16,13 @@
 use problemreductions::export::*;
 use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create IS instance: path graph P4
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(4, edges.clone());
+    // 1. Create IS instance: Petersen graph
+    let (num_vertices, edges) = petersen();
+    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&is);
diff --git a/examples/reduction_maximumindependentset_to_maximumsetpacking.rs b/examples/reduction_maximumindependentset_to_maximumsetpacking.rs
index 066bc6fd..8d073300 100644
--- a/examples/reduction_maximumindependentset_to_maximumsetpacking.rs
+++ b/examples/reduction_maximumindependentset_to_maximumsetpacking.rs
@@ -6,9 +6,9 @@
 //! incident edge sets, so IS maps to set packing with identical optimal value.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges: (0,1), (1,2), (2,3))
-//! - Source IS: max size 2
-//! - Target MaximumSetPacking: max packing 2
+//! - Instance: Petersen graph (10 vertices, 15 edges, 3-regular)
+//! - Source IS: max size 4
+//! - Target MaximumSetPacking: max packing 4
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumindependentset_to_maximumsetpacking.json` and `maximumindependentset_to_maximumsetpacking.result.json`.
@@ -17,17 +17,18 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("\n=== Independent Set -> Set Packing Reduction ===\n");
 
-    // Path graph P4: 0-1-2-3
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let source = MaximumIndependentSet::<SimpleGraph, i32>::new(4, edges.clone());
+    // Petersen graph: 10 vertices, 15 edges, 3-regular
+    let (num_vertices, edges) = petersen();
+    let source = MaximumIndependentSet::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
-    println!("Source: MaximumIndependentSet on P4");
-    println!("  Vertices: 4");
+    println!("Source: MaximumIndependentSet on Petersen graph");
+    println!("  Vertices: {}", num_vertices);
     println!("  Edges: {:?}", edges);
 
     // Reduce to MaximumSetPacking
@@ -75,8 +76,8 @@ fn main() {
     let source_size = source.solution_size(&source_sol);
     let target_size = target.solution_size(target_sol);
 
-    assert_eq!(source_size.size, 2, "IS on P4 has optimal size 2");
-    assert_eq!(target_size.size, 2, "MaximumSetPacking should also have size 2");
+    assert_eq!(source_size.size, 4, "IS on Petersen graph has optimal size 4");
+    assert_eq!(target_size.size, 4, "MaximumSetPacking should also have size 4");
 
     // Export JSON
     let overhead = lookup_overhead("MaximumIndependentSet", "MaximumSetPacking")
@@ -107,5 +108,5 @@ fn main() {
     let name = env!("CARGO_BIN_NAME").strip_prefix("reduction_").unwrap();
     write_example(name, &data, &results);
 
-    println!("\nDone: IS(P4) optimal=2 maps to MaximumSetPacking optimal=2");
+    println!("\nDone: IS(Petersen) optimal=4 maps to MaximumSetPacking optimal=4");
 }
diff --git a/examples/reduction_maximumindependentset_to_minimumvertexcover.rs b/examples/reduction_maximumindependentset_to_minimumvertexcover.rs
index 39539dc1..19f23f0f 100644
--- a/examples/reduction_maximumindependentset_to_minimumvertexcover.rs
+++ b/examples/reduction_maximumindependentset_to_minimumvertexcover.rs
@@ -5,9 +5,9 @@
 //! operation preserves optimality since |IS| + |VC| = |V| is constant.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges)
-//! - Source IS: max size 2 (e.g., {0, 2} or {0, 3} or {1, 3})
-//! - Target VC: min size 2
+//! - Instance: Petersen graph (10 vertices, 15 edges, 3-regular)
+//! - Source IS: max size 4
+//! - Target VC: min size 6
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumindependentset_to_minimumvertexcover.json` and `maximumindependentset_to_minimumvertexcover.result.json`.
@@ -16,12 +16,13 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create IS instance: path graph P4
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(4, edges.clone());
+    // 1. Create IS instance: Petersen graph
+    let (num_vertices, edges) = petersen();
+    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to VC
     let reduction = ReduceTo::<MinimumVertexCover<SimpleGraph, i32>>::reduce_to(&is);
diff --git a/examples/reduction_maximumindependentset_to_qubo.rs b/examples/reduction_maximumindependentset_to_qubo.rs
index 1ef3518b..f830e9c6 100644
--- a/examples/reduction_maximumindependentset_to_qubo.rs
+++ b/examples/reduction_maximumindependentset_to_qubo.rs
@@ -11,10 +11,10 @@
 //! independent set size while respecting adjacency constraints.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 with 4 vertices and 3 edges (0-1-2-3)
-//! - Source: MaximumIndependentSet with maximum size 2 (e.g., {0,2} or {1,3})
-//! - QUBO variables: 4 (one per vertex)
-//! - Expected: Two optimal solutions of size 2: vertices {0,2} and {1,3}
+//! - Instance: Petersen graph (10 vertices, 15 edges, 3-regular)
+//! - Source: MaximumIndependentSet with maximum size 4
+//! - QUBO variables: 10 (one per vertex)
+//! - Expected: Optimal solutions of size 4
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumindependentset_to_qubo.json` and `maximumindependentset_to_qubo.result.json`.
@@ -26,20 +26,21 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("=== Independent Set -> QUBO Reduction ===\n");
 
-    // Path graph P4: 0-1-2-3
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(4, edges.clone());
+    // Petersen graph: 10 vertices, 15 edges, 3-regular
+    let (num_vertices, edges) = petersen();
+    let is = MaximumIndependentSet::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&is);
     let qubo = reduction.target_problem();
 
-    println!("Source: MaximumIndependentSet on path P4 (4 vertices, 3 edges)");
+    println!("Source: MaximumIndependentSet on Petersen graph (10 vertices, 15 edges)");
     println!("Target: QUBO with {} variables", qubo.num_variables());
     println!("Q matrix:");
     for row in qubo.matrix() {
diff --git a/examples/reduction_maximummatching_to_ilp.rs b/examples/reduction_maximummatching_to_ilp.rs
index 757a7036..f9ad063d 100644
--- a/examples/reduction_maximummatching_to_ilp.rs
+++ b/examples/reduction_maximummatching_to_ilp.rs
@@ -6,9 +6,9 @@
 //! Objective: maximize sum of w_e * x_e.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges: 0-1, 1-2, 2-3)
-//! - Source MaximumMatching: max matching size 2 (e.g., {0-1, 2-3})
-//! - Target ILP: 3 binary variables (one per edge), 4 vertex constraints
+//! - Instance: Petersen graph (10 vertices, 15 edges), perfect matching of size 5
+//! - Source MaximumMatching: max matching size 5
+//! - Target ILP: 15 binary variables (one per edge), 10 vertex constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximummatching_to_ilp.json` and `maximummatching_to_ilp.result.json`.
@@ -16,12 +16,13 @@
 use problemreductions::export::*;
 use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create MaximumMatching instance: path graph P4 with unit weights
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let matching = MaximumMatching::<SimpleGraph, i32>::unweighted(4, edges.clone());
+    // 1. Create MaximumMatching instance: Petersen graph with unit weights
+    let (num_vertices, edges) = petersen();
+    let matching = MaximumMatching::<SimpleGraph, i32>::unweighted(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&matching);
diff --git a/examples/reduction_maximummatching_to_maximumsetpacking.rs b/examples/reduction_maximummatching_to_maximumsetpacking.rs
index 6b4105f9..644c27a9 100644
--- a/examples/reduction_maximummatching_to_maximumsetpacking.rs
+++ b/examples/reduction_maximummatching_to_maximumsetpacking.rs
@@ -6,9 +6,9 @@
 //! no shared elements) with the same weight.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges) with unit weights
-//! - Source matching: max size 2 (e.g., edges {(0,1), (2,3)})
-//! - Target MaximumSetPacking: max packing 2
+//! - Instance: Petersen graph (10 vertices, 15 edges), perfect matching of size 5
+//! - Source matching: max size 5
+//! - Target MaximumSetPacking: max packing 5
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximummatching_to_maximumsetpacking.json` and `maximummatching_to_maximumsetpacking.result.json`.
@@ -17,17 +17,18 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("\n=== MaximumMatching -> Set Packing Reduction ===\n");
 
-    // Path graph P4: 0-1-2-3 with unit weights
-    let edges = vec![(0, 1), (1, 2), (2, 3)];
-    let source = MaximumMatching::<SimpleGraph, i32>::unweighted(4, edges.clone());
+    // Petersen graph with unit weights
+    let (num_vertices, edges) = petersen();
+    let source = MaximumMatching::<SimpleGraph, i32>::unweighted(num_vertices, edges.clone());
 
-    println!("Source: MaximumMatching on P4");
-    println!("  Vertices: 4");
+    println!("Source: MaximumMatching on Petersen graph");
+    println!("  Vertices: {}", num_vertices);
     println!("  Edges: {:?}", edges);
 
     // Reduce to MaximumSetPacking
@@ -94,5 +95,5 @@ fn main() {
     let name = env!("CARGO_BIN_NAME").strip_prefix("reduction_").unwrap();
     write_example(name, &data, &results);
 
-    println!("\nDone: MaximumMatching(P4) optimal=2 maps to MaximumSetPacking optimal=2");
+    println!("\nDone: MaximumMatching(Petersen) optimal=5 maps to MaximumSetPacking optimal=5");
 }
diff --git a/examples/reduction_maximumsetpacking_to_ilp.rs b/examples/reduction_maximumsetpacking_to_ilp.rs
index 71d99af0..a55b06f5 100644
--- a/examples/reduction_maximumsetpacking_to_ilp.rs
+++ b/examples/reduction_maximumsetpacking_to_ilp.rs
@@ -6,10 +6,10 @@
 //! Objective: maximize sum of w_i * x_i.
 //!
 //! ## This Example
-//! - Instance: 3 sets: S0={0,1}, S1={1,2}, S2={2,3,4}
-//!   Overlapping pairs: (S0,S1) share element 1, (S1,S2) share element 2
-//! - Source MaximumSetPacking: max packing size 2 (S0 and S2 are disjoint)
-//! - Target ILP: 3 binary variables, 2 overlap constraints
+//! - Instance: 6 sets over universe {0,...,7}
+//!   - S0={0,1,2}, S1={2,3,4}, S2={4,5,6}, S3={6,7,0}, S4={1,3,5}, S5={0,4,7}
+//! - Source MaximumSetPacking: max packing size 2 (e.g., S0 and S2, or S1 and S3)
+//! - Target ILP: 6 binary variables, overlap constraints for each pair sharing elements
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumsetpacking_to_ilp.json` and `maximumsetpacking_to_ilp.result.json`.
@@ -19,12 +19,16 @@ use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
 
 fn main() {
-    // 1. Create MaximumSetPacking instance: 3 sets
-    let sp = MaximumSetPacking::<i32>::new(vec![
-        vec![0, 1],
-        vec![1, 2],
-        vec![2, 3, 4],
-    ]);
+    // 1. Create MaximumSetPacking instance: 6 sets over universe {0,...,7}
+    let sets = vec![
+        vec![0, 1, 2],    // S0
+        vec![2, 3, 4],    // S1 (overlaps S0 at 2)
+        vec![4, 5, 6],    // S2 (overlaps S1 at 4)
+        vec![6, 7, 0],    // S3 (overlaps S2 at 6, S0 at 0)
+        vec![1, 3, 5],    // S4 (overlaps S0, S1, S2)
+        vec![0, 4, 7],    // S5 (overlaps S0, S1, S3)
+    ];
+    let sp = MaximumSetPacking::<i32>::new(sets.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&sp);
@@ -32,7 +36,10 @@ fn main() {
 
     // 3. Print transformation
     println!("\n=== Problem Transformation ===");
-    println!("Source: MaximumSetPacking with {} variables", sp.num_variables());
+    println!("Source: MaximumSetPacking with {} sets over universe {{0,...,7}}", sp.num_variables());
+    for (i, s) in sets.iter().enumerate() {
+        println!("  S{} = {:?}", i, s);
+    }
     println!("Target: ILP with {} variables, {} constraints", ilp.num_vars, ilp.constraints.len());
 
     // 4. Solve target ILP
diff --git a/examples/reduction_maximumsetpacking_to_qubo.rs b/examples/reduction_maximumsetpacking_to_qubo.rs
index 8f86567e..04fa5568 100644
--- a/examples/reduction_maximumsetpacking_to_qubo.rs
+++ b/examples/reduction_maximumsetpacking_to_qubo.rs
@@ -11,21 +11,22 @@
 //! non-overlapping sets.
 //!
 //! ## This Example
-//! - Instance: 3 sets over universe {0,1,2,3,4}
-//!   - Set A = {0, 1}
-//!   - Set B = {1, 2}
-//!   - Set C = {2, 3, 4}
-//! - Overlaps: A-B share element 1, B-C share element 2
-//! - QUBO variables: 3 (one per set)
-//! - Expected: Optimal packing selects {A, C} (size 2), since A and C
-//!   do not overlap
+//! - Instance: 6 sets over universe {0,...,7}
+//!   - S0 = {0, 1, 2}
+//!   - S1 = {2, 3, 4}   (overlaps S0 at 2)
+//!   - S2 = {4, 5, 6}   (overlaps S1 at 4)
+//!   - S3 = {6, 7, 0}   (overlaps S2 at 6, S0 at 0)
+//!   - S4 = {1, 3, 5}   (overlaps S0, S1, S2)
+//!   - S5 = {0, 4, 7}   (overlaps S0, S1, S3)
+//! - QUBO variables: 6 (one per set)
+//! - Expected: Optimal packing selects 2 disjoint sets (e.g., {S0, S2} or {S1, S3})
 //!
 //! ## Output
 //! Exports `docs/paper/examples/maximumsetpacking_to_qubo.json` and `maximumsetpacking_to_qubo.result.json`.
 //!
 //! ## Usage
 //! ```bash
-//! cargo run --example reduction_setpacking_to_qubo
+//! cargo run --example reduction_maximumsetpacking_to_qubo
 //! ```
 
 use problemreductions::export::*;
@@ -34,21 +35,25 @@ use problemreductions::prelude::*;
 fn main() {
     println!("=== Set Packing -> QUBO Reduction ===\n");
 
-    // 3 sets over universe {0,1,2,3,4}
+    // 6 sets over universe {0,...,7}
     let sets = vec![
-        vec![0, 1],    // Set A
-        vec![1, 2],    // Set B
-        vec![2, 3, 4], // Set C
+        vec![0, 1, 2],    // S0
+        vec![2, 3, 4],    // S1 (overlaps S0 at 2)
+        vec![4, 5, 6],    // S2 (overlaps S1 at 4)
+        vec![6, 7, 0],    // S3 (overlaps S2 at 6, S0 at 0)
+        vec![1, 3, 5],    // S4 (overlaps S0, S1, S2)
+        vec![0, 4, 7],    // S5 (overlaps S0, S1, S3)
     ];
-    let set_names = ["Set-A".to_string(), "Set-B".to_string(), "Set-C".to_string()];
     let sp = MaximumSetPacking::<i32>::new(sets.clone());
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&sp);
     let qubo = reduction.target_problem();
 
-    println!("Source: MaximumSetPacking with 3 sets over universe {{0,1,2,3,4}}");
-    println!("  Set A = {{0, 1}}, Set B = {{1, 2}}, Set C = {{2, 3, 4}}");
+    println!("Source: MaximumSetPacking with 6 sets over universe {{0,...,7}}");
+    for (i, s) in sets.iter().enumerate() {
+        println!("  S{} = {:?}", i, s);
+    }
     println!("Target: QUBO with {} variables", qubo.num_variables());
     println!("Q matrix:");
     for row in qubo.matrix() {
@@ -64,14 +69,17 @@ fn main() {
     let mut solutions = Vec::new();
     for sol in &qubo_solutions {
         let extracted = reduction.extract_solution(sol);
-        let selected: Vec<String> = extracted
+        let selected: Vec<usize> = extracted
             .iter()
             .enumerate()
             .filter(|(_, &x)| x == 1)
-            .map(|(i, _)| set_names[i].clone())
+            .map(|(i, _)| i)
             .collect();
         let packing_size = selected.len();
-        println!("  Selected: {:?} (packing size {})", selected, packing_size);
+        println!(
+            "  Selected sets: {:?} (packing size {})",
+            selected, packing_size
+        );
 
         // Closed-loop verification: check solution is valid in original problem
         let sol_size = sp.solution_size(&extracted);
diff --git a/examples/reduction_minimumdominatingset_to_ilp.rs b/examples/reduction_minimumdominatingset_to_ilp.rs
index 708dc17e..c7d02105 100644
--- a/examples/reduction_minimumdominatingset_to_ilp.rs
+++ b/examples/reduction_minimumdominatingset_to_ilp.rs
@@ -6,9 +6,9 @@
 //! Objective: minimize sum of w_v * x_v.
 //!
 //! ## This Example
-//! - Instance: Path graph P4 (4 vertices, 3 edges: 0-1-2-3)
-//! - Source MinimumDominatingSet: min dominating set size 2 (e.g., {1,2})
-//! - Target ILP: 4 binary variables, 4 domination constraints
+//! - Instance: Petersen graph (10 vertices, 15 edges), min dominating set size 3
+//! - Source MinimumDominatingSet: min dominating set size 3
+//! - Target ILP: 10 binary variables, 10 domination constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumdominatingset_to_ilp.json` and `minimumdominatingset_to_ilp.result.json`.
@@ -16,11 +16,13 @@
 use problemreductions::export::*;
 use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create MinimumDominatingSet instance: path graph P4
-    let ds = MinimumDominatingSet::<SimpleGraph, i32>::new(4, vec![(0, 1), (1, 2), (2, 3)]);
+    // 1. Create MinimumDominatingSet instance: Petersen graph
+    let (num_vertices, edges) = petersen();
+    let ds = MinimumDominatingSet::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&ds);
diff --git a/examples/reduction_minimumsetcovering_to_ilp.rs b/examples/reduction_minimumsetcovering_to_ilp.rs
index 0a5771e2..00e3f33f 100644
--- a/examples/reduction_minimumsetcovering_to_ilp.rs
+++ b/examples/reduction_minimumsetcovering_to_ilp.rs
@@ -6,9 +6,10 @@
 //! Objective: minimize sum of w_i * x_i.
 //!
 //! ## This Example
-//! - Instance: Universe size 3, sets: S0={0,1}, S1={1,2}, S2={0,2}
-//! - Source MinimumSetCovering: min cover size 2 (any two sets cover all elements)
-//! - Target ILP: 3 binary variables, 3 element-coverage constraints
+//! - Instance: Universe size 8, 6 sets
+//!   - S0={0,1,2}, S1={2,3,4}, S2={4,5,6}, S3={6,7,0}, S4={1,3,5}, S5={0,4,7}
+//! - Source MinimumSetCovering: every element in {0,...,7} must be covered
+//! - Target ILP: 6 binary variables, 8 element-coverage constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumsetcovering_to_ilp.json` and `minimumsetcovering_to_ilp.result.json`.
@@ -18,15 +19,16 @@ use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
 
 fn main() {
-    // 1. Create MinimumSetCovering instance: universe {0,1,2}, 3 sets
-    let sc = MinimumSetCovering::<i32>::new(
-        3,
-        vec![
-            vec![0, 1],
-            vec![1, 2],
-            vec![0, 2],
-        ],
-    );
+    // 1. Create MinimumSetCovering instance: universe {0,...,7}, 6 sets
+    let sets = vec![
+        vec![0, 1, 2],    // S0
+        vec![2, 3, 4],    // S1
+        vec![4, 5, 6],    // S2
+        vec![6, 7, 0],    // S3
+        vec![1, 3, 5],    // S4
+        vec![0, 4, 7],    // S5
+    ];
+    let sc = MinimumSetCovering::<i32>::new(8, sets.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&sc);
@@ -34,7 +36,10 @@ fn main() {
 
     // 3. Print transformation
     println!("\n=== Problem Transformation ===");
-    println!("Source: MinimumSetCovering with {} variables", sc.num_variables());
+    println!("Source: MinimumSetCovering with {} sets over universe {{0,...,7}}", sc.num_variables());
+    for (i, s) in sets.iter().enumerate() {
+        println!("  S{} = {:?}", i, s);
+    }
     println!("Target: ILP with {} variables, {} constraints", ilp.num_vars, ilp.constraints.len());
 
     // 4. Solve target ILP
diff --git a/examples/reduction_minimumvertexcover_to_ilp.rs b/examples/reduction_minimumvertexcover_to_ilp.rs
index d70164fc..5e6baa8f 100644
--- a/examples/reduction_minimumvertexcover_to_ilp.rs
+++ b/examples/reduction_minimumvertexcover_to_ilp.rs
@@ -6,9 +6,9 @@
 //! Objective: minimize sum of w_v * x_v.
 //!
 //! ## This Example
-//! - Instance: Cycle C4 (4 vertices, 4 edges: 0-1-2-3-0)
-//! - Source VC: min cover size 2
-//! - Target ILP: 4 binary variables, 4 constraints
+//! - Instance: Petersen graph (10 vertices, 15 edges), VC=6
+//! - Source VC: min cover size 6
+//! - Target ILP: 10 binary variables, 15 constraints
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumvertexcover_to_ilp.json` and `minimumvertexcover_to_ilp.result.json`.
@@ -16,11 +16,13 @@
 use problemreductions::export::*;
 use problemreductions::prelude::*;
 use problemreductions::solvers::BruteForceFloat;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create VC instance: cycle C4
-    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(4, vec![(0, 1), (1, 2), (2, 3), (3, 0)]);
+    // 1. Create VC instance: Petersen graph (10 vertices, 15 edges), VC=6
+    let (num_vertices, edges) = petersen();
+    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // 2. Reduce to ILP
     let reduction = ReduceTo::<ILP>::reduce_to(&vc);
diff --git a/examples/reduction_minimumvertexcover_to_maximumindependentset.rs b/examples/reduction_minimumvertexcover_to_maximumindependentset.rs
index fa4f6ed9..e31ca128 100644
--- a/examples/reduction_minimumvertexcover_to_maximumindependentset.rs
+++ b/examples/reduction_minimumvertexcover_to_maximumindependentset.rs
@@ -6,9 +6,9 @@
 //! computes the complement: IS = V \ VC.
 //!
 //! ## This Example
-//! - Instance: Cycle C4 (4 vertices, 4 edges)
-//! - Source VC: min size 2
-//! - Target IS: max size 2
+//! - Instance: Petersen graph (10 vertices, 15 edges), VC=6
+//! - Source VC: min size 6
+//! - Target IS: max size 4
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumvertexcover_to_maximumindependentset.json` and `minimumvertexcover_to_maximumindependentset.result.json`.
@@ -17,10 +17,13 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(4, vec![(0, 1), (1, 2), (2, 3), (0, 3)]);
+    // Petersen graph: 10 vertices, 15 edges, VC=6
+    let (num_vertices, edges) = petersen();
+    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     let reduction = ReduceTo::<MaximumIndependentSet<SimpleGraph, i32>>::reduce_to(&vc);
     let is = reduction.target_problem();
diff --git a/examples/reduction_minimumvertexcover_to_minimumsetcovering.rs b/examples/reduction_minimumvertexcover_to_minimumsetcovering.rs
index c4f46b4f..c435a83a 100644
--- a/examples/reduction_minimumvertexcover_to_minimumsetcovering.rs
+++ b/examples/reduction_minimumvertexcover_to_minimumsetcovering.rs
@@ -6,9 +6,9 @@
 //! in the cover) maps to a set cover (every universe element in some set).
 //!
 //! ## This Example
-//! - Instance: Triangle K3 (3 vertices, 3 edges)
-//! - Source VC: min size 2
-//! - Target MinimumSetCovering: min cover 2
+//! - Instance: Petersen graph (10 vertices, 15 edges), VC=6
+//! - Source VC: min size 6
+//! - Target MinimumSetCovering: min cover 6
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumvertexcover_to_minimumsetcovering.json` and `minimumvertexcover_to_minimumsetcovering.result.json`.
@@ -17,17 +17,18 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("\n=== Vertex Cover -> Set Covering Reduction ===\n");
 
-    // Triangle K3: 3 vertices, 3 edges
-    let edges = vec![(0, 1), (1, 2), (0, 2)];
-    let source = MinimumVertexCover::<SimpleGraph, i32>::new(3, edges.clone());
+    // Petersen graph: 10 vertices, 15 edges, VC=6
+    let (num_vertices, edges) = petersen();
+    let source = MinimumVertexCover::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
-    println!("Source: MinimumVertexCover on K3");
-    println!("  Vertices: 3");
+    println!("Source: MinimumVertexCover on Petersen graph");
+    println!("  Vertices: {}", num_vertices);
     println!("  Edges: {:?}", edges);
 
     // Reduce to MinimumSetCovering
@@ -76,8 +77,8 @@ fn main() {
     let source_size = source.solution_size(&source_sol);
     let target_size = target.solution_size(target_sol);
 
-    assert_eq!(source_size.size, 2, "VC on K3 has optimal size 2");
-    assert_eq!(target_size.size, 2, "MinimumSetCovering should also have size 2");
+    assert_eq!(source_size.size, 6, "VC on Petersen has optimal size 6");
+    assert_eq!(target_size.size, 6, "MinimumSetCovering should also have size 6");
 
     // Export JSON
     let overhead = lookup_overhead("MinimumVertexCover", "MinimumSetCovering")
@@ -109,5 +110,5 @@ fn main() {
     let name = env!("CARGO_BIN_NAME").strip_prefix("reduction_").unwrap();
     write_example(name, &data, &results);
 
-    println!("\nDone: VC(K3) optimal=2 maps to MinimumSetCovering optimal=2");
+    println!("\nDone: VC(Petersen) optimal=6 maps to MinimumSetCovering optimal=6");
 }
diff --git a/examples/reduction_minimumvertexcover_to_qubo.rs b/examples/reduction_minimumvertexcover_to_qubo.rs
index 2be17d63..c244f787 100644
--- a/examples/reduction_minimumvertexcover_to_qubo.rs
+++ b/examples/reduction_minimumvertexcover_to_qubo.rs
@@ -11,10 +11,10 @@
 //! number of selected vertices while covering all edges.
 //!
 //! ## This Example
-//! - Instance: Cycle graph C4 with 4 vertices and 4 edges (0-1-2-3-0)
-//! - Source: MinimumVertexCover with minimum size 2
-//! - QUBO variables: 4 (one per vertex)
-//! - Expected: Optimal vertex covers of size 2 (e.g., {0,2} or {1,3})
+//! - Instance: Petersen graph (10 vertices, 15 edges), VC=6
+//! - Source: MinimumVertexCover with minimum size 6
+//! - QUBO variables: 10 (one per vertex)
+//! - Expected: Optimal vertex covers of size 6
 //!
 //! ## Output
 //! Exports `docs/paper/examples/minimumvertexcover_to_qubo.json` and `minimumvertexcover_to_qubo.result.json`.
@@ -26,20 +26,21 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
     println!("=== Vertex Covering -> QUBO Reduction ===\n");
 
-    // Cycle C4: 0-1-2-3-0
-    let edges = vec![(0, 1), (1, 2), (2, 3), (0, 3)];
-    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(4, edges.clone());
+    // Petersen graph: 10 vertices, 15 edges, VC=6
+    let (num_vertices, edges) = petersen();
+    let vc = MinimumVertexCover::<SimpleGraph, i32>::new(num_vertices, edges.clone());
 
     // Reduce to QUBO
     let reduction = ReduceTo::<QUBO>::reduce_to(&vc);
     let qubo = reduction.target_problem();
 
-    println!("Source: MinimumVertexCover on cycle C4 (4 vertices, 4 edges)");
+    println!("Source: MinimumVertexCover on Petersen graph (10 vertices, 15 edges)");
     println!("Target: QUBO with {} variables", qubo.num_variables());
     println!("Q matrix:");
     for row in qubo.matrix() {
@@ -77,12 +78,12 @@ fn main() {
         });
     }
 
-    // All optimal solutions should have size 2
+    // All optimal solutions should have size 6
     assert!(
-        solutions.iter().all(|s| s.source_config.iter().filter(|&&x| x == 1).count() == 2),
-        "All optimal VC solutions on C4 should have size 2"
+        solutions.iter().all(|s| s.source_config.iter().filter(|&&x| x == 1).count() == 6),
+        "All optimal VC solutions on Petersen graph should have size 6"
     );
-    println!("\nVerification passed: all solutions are valid with size 2");
+    println!("\nVerification passed: all solutions are valid with size 6");
 
     // Export JSON
     let overhead = lookup_overhead("MinimumVertexCover", "QUBO")
diff --git a/examples/reduction_qubo_to_spinglass.rs b/examples/reduction_qubo_to_spinglass.rs
index 6aab182b..259c596b 100644
--- a/examples/reduction_qubo_to_spinglass.rs
+++ b/examples/reduction_qubo_to_spinglass.rs
@@ -6,9 +6,9 @@
 //! Q_{ij} = -4J_{ij} for off-diagonal and Q_{ii} = 2*sum_j J_{ij} - 2h_i for diagonal.
 //!
 //! ## This Example
-//! - Instance: 3-variable QUBO with diagonal [-1, -2, -1] and coupling Q_{01} = 3
-//! - Source QUBO: 3 binary variables
-//! - Target SpinGlass: 3 spins
+//! - Instance: 10-variable QUBO with Petersen connectivity
+//! - Source QUBO: 10 binary variables
+//! - Target SpinGlass: 10 spins
 //!
 //! ## Output
 //! Exports `docs/paper/examples/qubo_to_spinglass.json` and
@@ -18,14 +18,21 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    let matrix = vec![
-        vec![-1.0, 3.0, 0.0],
-        vec![0.0, -2.0, 0.0],
-        vec![0.0, 0.0, -1.0],
-    ];
+    let (n, edges) = petersen();
+    let mut matrix = vec![vec![0.0; n]; n];
+    // Diagonal: linear terms
+    for (i, row) in matrix.iter_mut().enumerate() {
+        row[i] = -1.0 + 0.2 * i as f64;
+    }
+    // Off-diagonal: quadratic terms on Petersen edges
+    for (idx, &(u, v)) in edges.iter().enumerate() {
+        let (i, j) = if u < v { (u, v) } else { (v, u) };
+        matrix[i][j] = if idx % 2 == 0 { 2.0 } else { -1.5 };
+    }
     let qubo = QUBO::from_matrix(matrix.clone());
 
     let reduction = ReduceTo::<SpinGlass<SimpleGraph, f64>>::reduce_to(&qubo);
diff --git a/examples/reduction_satisfiability_to_kcoloring.rs b/examples/reduction_satisfiability_to_kcoloring.rs
index e1b8523e..88cc520c 100644
--- a/examples/reduction_satisfiability_to_kcoloring.rs
+++ b/examples/reduction_satisfiability_to_kcoloring.rs
@@ -6,9 +6,11 @@
 //! phi is satisfiable iff the constructed graph is 3-colorable.
 //!
 //! ## This Example
-//! - Instance: phi = (x1 v x2), 2 vars, 1 clause
-//! - Source SAT: satisfiable
-//! - Target: 3-Coloring with larger graph
+//! - Instance: 5-variable, 3-clause SAT formula with unit clauses
+//!   (OR-gadgets add 5 vertices per extra literal per clause, making BruteForce
+//!   infeasible for multi-literal clauses; unit clauses keep it at 13 vertices)
+//! - Source SAT: satisfiable (x1=1, x3=0, x5=1, x2/x4 free)
+//! - Target: 3-Coloring with 13 vertices
 //!
 //! ## Output
 //! Exports `docs/paper/examples/satisfiability_to_kcoloring.json` and `satisfiability_to_kcoloring.result.json`.
@@ -18,18 +20,24 @@ use problemreductions::prelude::*;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create SAT instance: phi = (x1 v x2), 2 variables, 1 clause
+    // 1. Create SAT instance: 5-variable, 3-clause formula with unit clauses
+    //    The SAT→KColoring reduction creates OR-gadgets that add 5 vertices per literal
+    //    beyond the first in each clause. BruteForce on 3-coloring is O(3^n), so we use
+    //    unit clauses (1 literal each) to keep vertex count at 2*5+3 = 13 (3^13 ~ 1.6M).
     let sat = Satisfiability::<i32>::new(
-        2,
+        5,
         vec![
-            CNFClause::new(vec![1, 2]), // x1 OR x2
+            CNFClause::new(vec![1]),     // x1 (unit clause)
+            CNFClause::new(vec![-3]),    // ~x3 (unit clause)
+            CNFClause::new(vec![5]),     // x5 (unit clause)
         ],
     );
 
     println!("=== SAT to 3-Coloring Reduction (Garey & Johnson 1979) ===\n");
-    println!("Source SAT formula:");
-    println!("  (x1 v x2)");
+    println!("Source SAT formula: 5-variable, 3-clause SAT (unit clauses to fit BruteForce)");
+    println!("  (x1) ^ (~x3) ^ (x5)");
     println!("  {} variables, {} clauses", sat.num_vars(), sat.num_clauses());
+    println!("  (Unit clauses avoid OR-gadgets, keeping vertex count manageable for BruteForce)");
 
     // 2. Reduce to 3-Coloring
     //    SAT reduces to KColoring<3, SimpleGraph, i32>
@@ -57,8 +65,8 @@ fn main() {
     let sat_solution = reduction.extract_solution(&coloring_solutions[0]);
     println!("Extracted SAT solution: {:?}", sat_solution);
     println!(
-        "  Interpretation: x1={}, x2={}",
-        sat_solution[0], sat_solution[1]
+        "  Interpretation: x1={}, x2={}, x3={}, x4={}, x5={}",
+        sat_solution[0], sat_solution[1], sat_solution[2], sat_solution[3], sat_solution[4]
     );
 
     let size = sat.solution_size(&sat_solution);
diff --git a/examples/reduction_satisfiability_to_ksatisfiability.rs b/examples/reduction_satisfiability_to_ksatisfiability.rs
index 523206d1..75a6b7b7 100644
--- a/examples/reduction_satisfiability_to_ksatisfiability.rs
+++ b/examples/reduction_satisfiability_to_ksatisfiability.rs
@@ -6,9 +6,12 @@
 //! in a chain that preserves satisfiability.
 //!
 //! ## This Example
-//! - Instance: phi = (x1) ^ (x1 v x2 v x3 v x4), 4 vars, 2 clauses
-//!   - First clause has 1 literal (will be padded to 3)
-//!   - Second clause has 4 literals (will be split into two 3-literal clauses)
+//! - Instance: 5-variable, 6-clause SAT formula with mixed clause sizes (1, 2, 3, 3, 4, 5 literals)
+//!   - 1-literal clause: padded to 3
+//!   - 2-literal clause: padded to 3
+//!   - 3-literal clauses: no change
+//!   - 4-literal clause: split into two 3-literal clauses
+//!   - 5-literal clause: split into three 3-literal clauses
 //! - Source SAT: satisfiable
 //! - Target: 3-SAT with 3 literals per clause
 //!
@@ -19,23 +22,30 @@ use problemreductions::export::*;
 use problemreductions::prelude::*;
 
 fn main() {
-    // 1. Create SAT instance with varied clause sizes:
-    //    phi = (x1) ^ (x1 v x2 v x3 v x4)
-    //    - Clause 1 has 1 literal (needs padding to 3)
-    //    - Clause 2 has 4 literals (needs splitting into 3-literal clauses)
+    // 1. Create SAT instance with varied clause sizes to demonstrate padding and splitting:
+    //    - 1 literal: padded to 3
+    //    - 2 literals: padded to 3
+    //    - 3 literals: no change (already 3-SAT)
+    //    - 4 literals: split into two 3-literal clauses
+    //    - 5 literals: split into three 3-literal clauses
     let sat = Satisfiability::<i32>::new(
-        4,
+        5,
         vec![
-            CNFClause::new(vec![1]),          // x1 (1 literal - will be padded)
-            CNFClause::new(vec![1, 2, 3, 4]), // x1 v x2 v x3 v x4 (4 literals - will be split)
+            CNFClause::new(vec![1]),              // 1 literal - will be padded
+            CNFClause::new(vec![2, -3]),           // 2 literals - will be padded
+            CNFClause::new(vec![-1, 3, 4]),        // 3 literals - no change
+            CNFClause::new(vec![2, -4, 5]),        // 3 literals - no change
+            CNFClause::new(vec![1, -2, 3, -5]),    // 4 literals - will be split
+            CNFClause::new(vec![-1, 2, -3, 4, 5]), // 5 literals - will be split
         ],
     );
 
     println!("=== SAT to 3-SAT Reduction ===\n");
-    println!("Source SAT formula:");
-    println!("  (x1) ^ (x1 v x2 v x3 v x4)");
+    println!("Source SAT formula: 5-variable, 6-clause SAT with mixed clause sizes");
+    println!("  (x1) ^ (x2 v ~x3) ^ (~x1 v x3 v x4) ^ (x2 v ~x4 v x5) ^");
+    println!("  (x1 v ~x2 v x3 v ~x5) ^ (~x1 v x2 v ~x3 v x4 v x5)");
     println!("  {} variables, {} clauses", sat.num_vars(), sat.num_clauses());
-    println!("  Clause sizes: 1 and 4 (both need transformation for 3-SAT)");
+    println!("  Clause sizes: 1, 2, 3, 3, 4, 5 (demonstrates padding and splitting)");
 
     // 2. Reduce to 3-SAT (K=3)
     let reduction = ReduceTo::<KSatisfiability<3, i32>>::reduce_to(&sat);
@@ -49,8 +59,10 @@ fn main() {
         ksat.num_clauses()
     );
     println!("  Additional variables: {} (ancilla/auxiliary)", ksat.num_vars() - sat.num_vars());
-    println!("  Clause (x1) padded: (x1 v a v b) ^ (x1 v a v ~b) ^ ... ");
-    println!("  Clause (x1 v x2 v x3 v x4) split: (x1 v x2 v c) ^ (~c v x3 v x4)");
+    println!("  1-literal (x1) padded: (x1 v a v b) ^ (x1 v a v ~b) ^ ... ");
+    println!("  2-literal (x2 v ~x3) padded similarly with auxiliary variables");
+    println!("  4-literal (x1 v ~x2 v x3 v ~x5) split: two 3-literal clauses");
+    println!("  5-literal (~x1 v x2 v ~x3 v x4 v x5) split: three 3-literal clauses");
 
     // Print target clauses
     println!("\n  Target 3-SAT clauses:");
@@ -68,8 +80,8 @@ fn main() {
     let sat_solution = reduction.extract_solution(&ksat_solutions[0]);
     println!("Extracted SAT solution: {:?}", sat_solution);
     println!(
-        "  Interpretation: x1={}, x2={}, x3={}, x4={}",
-        sat_solution[0], sat_solution[1], sat_solution[2], sat_solution[3]
+        "  Interpretation: x1={}, x2={}, x3={}, x4={}, x5={}",
+        sat_solution[0], sat_solution[1], sat_solution[2], sat_solution[3], sat_solution[4]
     );
 
     let size = sat.solution_size(&sat_solution);
diff --git a/examples/reduction_satisfiability_to_maximumindependentset.rs b/examples/reduction_satisfiability_to_maximumindependentset.rs
index 6b9991a4..c50c2a14 100644
--- a/examples/reduction_satisfiability_to_maximumindependentset.rs
+++ b/examples/reduction_satisfiability_to_maximumindependentset.rs
@@ -6,9 +6,9 @@
 //! edges connect complementary literals. phi is satisfiable iff G has IS of size m.
 //!
 //! ## This Example
-//! - Instance: phi = (x1 v x2) ^ (~x1 v x3) ^ (x2 v ~x3), 3 vars, 3 clauses
-//! - Source SAT: satisfiable (e.g., x1=1, x2=1, x3=1)
-//! - Target IS: size 3 (one vertex per clause)
+//! - Instance: 5-variable, 7-clause 3-SAT formula
+//! - Source SAT: satisfiable
+//! - Target IS: size 7 (one vertex per clause), 21 vertices total
 //!
 //! ## Output
 //! Exports `docs/paper/examples/satisfiability_to_maximumindependentset.json` and `satisfiability_to_maximumindependentset.result.json`.
@@ -18,20 +18,24 @@ use problemreductions::prelude::*;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create SAT instance: phi = (x1 v x2) ^ (~x1 v x3) ^ (x2 v ~x3)
-    //    3 variables, 3 clauses
+    // 1. Create SAT instance: 5-variable, 7-clause 3-SAT formula
     let sat = Satisfiability::<i32>::new(
-        3,
+        5,
         vec![
-            CNFClause::new(vec![1, 2]),   // x1 OR x2
-            CNFClause::new(vec![-1, 3]),  // NOT x1 OR x3
-            CNFClause::new(vec![2, -3]),  // x2 OR NOT x3
+            CNFClause::new(vec![1, 2, -3]),    // x1 v x2 v ~x3
+            CNFClause::new(vec![-1, 3, 4]),    // ~x1 v x3 v x4
+            CNFClause::new(vec![2, -4, 5]),    // x2 v ~x4 v x5
+            CNFClause::new(vec![-2, 3, -5]),   // ~x2 v x3 v ~x5
+            CNFClause::new(vec![1, -3, 5]),    // x1 v ~x3 v x5
+            CNFClause::new(vec![-1, -2, 4]),   // ~x1 v ~x2 v x4
+            CNFClause::new(vec![3, -4, -5]),   // x3 v ~x4 v ~x5
         ],
     );
 
     println!("=== SAT to Independent Set Reduction (Karp 1972) ===\n");
-    println!("Source SAT formula:");
-    println!("  (x1 v x2) ^ (~x1 v x3) ^ (x2 v ~x3)");
+    println!("Source SAT formula: 5-variable, 7-clause 3-SAT");
+    println!("  (x1 v x2 v ~x3) ^ (~x1 v x3 v x4) ^ (x2 v ~x4 v x5) ^");
+    println!("  (~x2 v x3 v ~x5) ^ (x1 v ~x3 v x5) ^ (~x1 v ~x2 v x4) ^ (x3 v ~x4 v ~x5)");
     println!("  {} variables, {} clauses", sat.num_vars(), sat.num_clauses());
 
     // 2. Reduce to Independent Set
@@ -59,8 +63,8 @@ fn main() {
     let sat_solution = reduction.extract_solution(&is_solutions[0]);
     println!("Extracted SAT solution: {:?}", sat_solution);
     println!(
-        "  Interpretation: x1={}, x2={}, x3={}",
-        sat_solution[0], sat_solution[1], sat_solution[2]
+        "  Interpretation: x1={}, x2={}, x3={}, x4={}, x5={}",
+        sat_solution[0], sat_solution[1], sat_solution[2], sat_solution[3], sat_solution[4]
     );
 
     let size = sat.solution_size(&sat_solution);
diff --git a/examples/reduction_satisfiability_to_minimumdominatingset.rs b/examples/reduction_satisfiability_to_minimumdominatingset.rs
index f290d50f..821fbeaf 100644
--- a/examples/reduction_satisfiability_to_minimumdominatingset.rs
+++ b/examples/reduction_satisfiability_to_minimumdominatingset.rs
@@ -6,9 +6,9 @@
 //! satisfiable iff the graph has a dominating set of size n.
 //!
 //! ## This Example
-//! - Instance: phi = (x1 v x2) ^ (~x1 v x2), 2 vars, 2 clauses
-//! - Source SAT: satisfiable (e.g., x2=1)
-//! - Target: Dominating set
+//! - Instance: 5-variable, 7-clause 3-SAT formula
+//! - Source SAT: satisfiable
+//! - Target: Dominating set with 3*5 + 7 = 22 vertices
 //!
 //! ## Output
 //! Exports `docs/paper/examples/satisfiability_to_minimumdominatingset.json` and `satisfiability_to_minimumdominatingset.result.json`.
@@ -18,18 +18,24 @@ use problemreductions::prelude::*;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    // 1. Create SAT instance: phi = (x1 v x2) ^ (~x1 v x2), 2 vars, 2 clauses
+    // 1. Create SAT instance: 5-variable, 7-clause 3-SAT formula
     let sat = Satisfiability::<i32>::new(
-        2,
+        5,
         vec![
-            CNFClause::new(vec![1, 2]),  // x1 OR x2
-            CNFClause::new(vec![-1, 2]), // NOT x1 OR x2
+            CNFClause::new(vec![1, 2, -3]),    // x1 v x2 v ~x3
+            CNFClause::new(vec![-1, 3, 4]),    // ~x1 v x3 v x4
+            CNFClause::new(vec![2, -4, 5]),    // x2 v ~x4 v x5
+            CNFClause::new(vec![-2, 3, -5]),   // ~x2 v x3 v ~x5
+            CNFClause::new(vec![1, -3, 5]),    // x1 v ~x3 v x5
+            CNFClause::new(vec![-1, -2, 4]),   // ~x1 v ~x2 v x4
+            CNFClause::new(vec![3, -4, -5]),   // x3 v ~x4 v ~x5
         ],
     );
 
     println!("=== SAT to Dominating Set Reduction (Garey & Johnson 1979) ===\n");
-    println!("Source SAT formula:");
-    println!("  (x1 v x2) ^ (~x1 v x2)");
+    println!("Source SAT formula: 5-variable, 7-clause 3-SAT");
+    println!("  (x1 v x2 v ~x3) ^ (~x1 v x3 v x4) ^ (x2 v ~x4 v x5) ^");
+    println!("  (~x2 v x3 v ~x5) ^ (x1 v ~x3 v x5) ^ (~x1 v ~x2 v x4) ^ (x3 v ~x4 v ~x5)");
     println!("  {} variables, {} clauses", sat.num_vars(), sat.num_clauses());
 
     // 2. Reduce to Dominating Set
@@ -45,7 +51,7 @@ fn main() {
     );
     println!("  Variable gadgets: 3 vertices per variable (pos, neg, dummy) forming triangles");
     println!("  Clause vertices: 1 per clause, connected to relevant literal vertices");
-    println!("  Layout: vertices 0-5 are variable gadgets, vertices 6-7 are clause vertices");
+    println!("  Layout: vertices 0-14 are variable gadgets (5 triangles), vertices 15-21 are clause vertices");
 
     // 3. Solve the target DS problem
     let solver = BruteForce::new();
@@ -57,8 +63,8 @@ fn main() {
     let sat_solution = reduction.extract_solution(&ds_solutions[0]);
     println!("Extracted SAT solution: {:?}", sat_solution);
     println!(
-        "  Interpretation: x1={}, x2={}",
-        sat_solution[0], sat_solution[1]
+        "  Interpretation: x1={}, x2={}, x3={}, x4={}, x5={}",
+        sat_solution[0], sat_solution[1], sat_solution[2], sat_solution[3], sat_solution[4]
     );
 
     let size = sat.solution_size(&sat_solution);
diff --git a/examples/reduction_spinglass_to_maxcut.rs b/examples/reduction_spinglass_to_maxcut.rs
index 1380f071..074936eb 100644
--- a/examples/reduction_spinglass_to_maxcut.rs
+++ b/examples/reduction_spinglass_to_maxcut.rs
@@ -6,9 +6,9 @@
 //! the Ising energy. When h_i != 0, an ancilla spin is added with w_{i,a} = h_i.
 //!
 //! ## This Example
-//! - Instance: 3-spin frustrated triangle (J_{01} = 1, J_{12} = 1, J_{02} = 1, h = 0)
-//! - Source SpinGlass: 3 spins, no external fields
-//! - Target MaxCut: 3 vertices (direct mapping, no ancilla)
+//! - Instance: Petersen graph with 10 spins, ±1 couplings, no external fields
+//! - Source SpinGlass: 10 spins on Petersen topology
+//! - Target MaxCut: 10 vertices (direct mapping, no ancilla)
 //!
 //! ## Output
 //! Exports `docs/paper/examples/spinglass_to_maxcut.json` and `spinglass_to_maxcut.result.json`.
@@ -17,14 +17,17 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    let sg = SpinGlass::<SimpleGraph, i32>::new(
-        3,
-        vec![((0, 1), 1), ((1, 2), 1), ((0, 2), 1)],
-        vec![0, 0, 0],
-    );
+    let (n, edges) = petersen();
+    let couplings: Vec<((usize, usize), i32)> = edges
+        .iter()
+        .enumerate()
+        .map(|(i, &(u, v))| ((u, v), if i % 2 == 0 { 1 } else { -1 }))
+        .collect();
+    let sg = SpinGlass::<SimpleGraph, i32>::new(n, couplings, vec![0; n]);
 
     let reduction = ReduceTo::<MaxCut<SimpleGraph, i32>>::reduce_to(&sg);
     let maxcut = reduction.target_problem();
diff --git a/examples/reduction_spinglass_to_qubo.rs b/examples/reduction_spinglass_to_qubo.rs
index 12131e88..67bf169b 100644
--- a/examples/reduction_spinglass_to_qubo.rs
+++ b/examples/reduction_spinglass_to_qubo.rs
@@ -6,11 +6,9 @@
 //! yields a QUBO objective Q(x) plus a constant offset.
 //!
 //! ## This Example
-//! - Instance: 3-spin antiferromagnetic chain with fields
-//!   - Couplings: J_{01} = -1.0, J_{12} = -1.0
-//!   - Fields: h = [0.5, -0.5, 0.5]
-//! - Source SpinGlass: 3 spins
-//! - Target QUBO: 3 binary variables
+//! - Instance: Petersen graph with 10 spins, 15 frustrated ±1 couplings, zero fields
+//! - Source SpinGlass: 10 spins on Petersen topology
+//! - Target QUBO: 10 binary variables
 //!
 //! ## Output
 //! Exports `docs/paper/examples/spinglass_to_qubo.json` and
@@ -20,14 +18,18 @@
 
 use problemreductions::export::*;
 use problemreductions::prelude::*;
+use problemreductions::topology::small_graphs::petersen;
 use problemreductions::topology::SimpleGraph;
 
 fn main() {
-    let sg = SpinGlass::<SimpleGraph, f64>::new(
-        3,
-        vec![((0, 1), -1.0), ((1, 2), -1.0)],
-        vec![0.5, -0.5, 0.5],
-    );
+    let (n, edges) = petersen();
+    // Alternating +/-1 couplings create frustration on odd cycles
+    let couplings: Vec<((usize, usize), f64)> = edges
+        .iter()
+        .enumerate()
+        .map(|(i, &(u, v))| ((u, v), if i % 2 == 0 { 1.0 } else { -1.0 }))
+        .collect();
+    let sg = SpinGlass::<SimpleGraph, f64>::new(n, couplings, vec![0.0; n]);
 
     let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&sg);
     let qubo = reduction.target_problem();
diff --git a/src/rules/ksatisfiability_qubo.rs b/src/rules/ksatisfiability_qubo.rs
index 5667b73b..49c01249 100644
--- a/src/rules/ksatisfiability_qubo.rs
+++ b/src/rules/ksatisfiability_qubo.rs
@@ -1,11 +1,15 @@
-//! Reduction from KSatisfiability (K=2) to QUBO (Max-2-SAT).
+//! Reduction from KSatisfiability to QUBO (Max-K-SAT).
 //!
-//! Each clause contributes to Q based on literal signs:
+//! For K=2 (quadratic penalty), each clause contributes to Q based on literal signs:
 //! - (x_i ∨ x_j): penalty (1-x_i)(1-x_j) → Q[i][i]-=1, Q[j][j]-=1, Q[i][j]+=1, const+=1
 //! - (¬x_i ∨ x_j): penalty x_i(1-x_j) → Q[i][i]+=1, Q[i][j]-=1
 //! - (x_i ∨ ¬x_j): penalty (1-x_i)x_j → Q[j][j]+=1, Q[i][j]-=1
 //! - (¬x_i ∨ ¬x_j): penalty x_i·x_j → Q[i][j]+=1
 //!
+//! For K≥3, we use the Rosenberg quadratization to reduce degree-K penalty terms
+//! to quadratic form by introducing auxiliary variables. Each clause of K literals
+//! requires K−2 auxiliary variables.
+//!
 //! CNFClause uses 1-indexed signed integers: positive = variable, negative = negated.
 
 use crate::models::optimization::QUBO;
@@ -17,10 +21,11 @@ use crate::rules::traits::{ReduceTo, ReductionResult};
 use crate::traits::Problem;
 use crate::types::ProblemSize;
 
-/// Result of reducing KSatisfiability<2> to QUBO.
+/// Result of reducing KSatisfiability to QUBO.
 #[derive(Debug, Clone)]
 pub struct ReductionKSatToQUBO {
     target: QUBO<f64>,
+    source_num_vars: usize,
     source_size: ProblemSize,
 }
 
@@ -33,7 +38,7 @@ impl ReductionResult for ReductionKSatToQUBO {
     }
 
     fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
-        target_solution.to_vec()
+        target_solution[..self.source_num_vars].to_vec()
     }
 
     fn source_size(&self) -> ProblemSize {
@@ -45,6 +50,264 @@ impl ReductionResult for ReductionKSatToQUBO {
     }
 }
 
+/// Result of reducing KSatisfiability<3> to QUBO.
+#[derive(Debug, Clone)]
+pub struct Reduction3SATToQUBO {
+    target: QUBO<f64>,
+    source_num_vars: usize,
+    source_size: ProblemSize,
+}
+
+impl ReductionResult for Reduction3SATToQUBO {
+    type Source = KSatisfiability<3, i32>;
+    type Target = QUBO<f64>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        target_solution[..self.source_num_vars].to_vec()
+    }
+
+    fn source_size(&self) -> ProblemSize {
+        self.source_size.clone()
+    }
+
+    fn target_size(&self) -> ProblemSize {
+        self.target.problem_size()
+    }
+}
+
+/// Convert a signed literal to (0-indexed variable, is_negated).
+fn lit_to_var(lit: i32) -> (usize, bool) {
+    let var = (lit.unsigned_abs() as usize) - 1;
+    let neg = lit < 0;
+    (var, neg)
+}
+
+/// Add the quadratic penalty term for a 2-literal clause to the QUBO matrix.
+///
+/// For clause (l_i ∨ l_j), the penalty for the clause being unsatisfied is
+/// the product of the complemented literals.
+fn add_2sat_clause_penalty(matrix: &mut [Vec<f64>], lits: &[i32]) {
+    assert_eq!(lits.len(), 2, "Expected 2-literal clause");
+
+    let (var_i, neg_i) = lit_to_var(lits[0]);
+    let (var_j, neg_j) = lit_to_var(lits[1]);
+
+    // Ensure i <= j for upper-triangular form
+    let (i, j, ni, nj) = if var_i <= var_j {
+        (var_i, var_j, neg_i, neg_j)
+    } else {
+        (var_j, var_i, neg_j, neg_i)
+    };
+
+    match (ni, nj) {
+        (false, false) => {
+            // (x_i ∨ x_j): penalty = (1-x_i)(1-x_j) = 1 - x_i - x_j + x_i·x_j
+            matrix[i][i] -= 1.0;
+            matrix[j][j] -= 1.0;
+            matrix[i][j] += 1.0;
+        }
+        (true, false) => {
+            // (¬x_i ∨ x_j): penalty = x_i(1-x_j) = x_i - x_i·x_j
+            matrix[i][i] += 1.0;
+            matrix[i][j] -= 1.0;
+        }
+        (false, true) => {
+            // (x_i ∨ ¬x_j): penalty = (1-x_i)x_j = x_j - x_i·x_j
+            matrix[j][j] += 1.0;
+            matrix[i][j] -= 1.0;
+        }
+        (true, true) => {
+            // (¬x_i ∨ ¬x_j): penalty = x_i·x_j
+            matrix[i][j] += 1.0;
+        }
+    }
+}
+
+/// Add the QUBO terms for a 3-literal clause using Rosenberg quadratization.
+///
+/// For clause (l1 ∨ l2 ∨ l3), the penalty for not satisfying the clause is:
+///   P = (1-l1)(1-l2)(1-l3) = y1·y2·y3
+/// where yi = 1 - li (complement of literal).
+///
+/// We introduce one auxiliary variable `a` and quadratize using the substitution
+/// a = y1·y2, adding penalty M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a) where M is a
+/// sufficiently large penalty (M = 2 suffices for Max-SAT).
+///
+/// The resulting quadratic form is:
+///   H = a·y3 + M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a)
+///
+/// `aux_var` is the 0-indexed auxiliary variable.
+fn add_3sat_clause_penalty(matrix: &mut [Vec<f64>], lits: &[i32], aux_var: usize) {
+    assert_eq!(lits.len(), 3, "Expected 3-literal clause");
+    let penalty = 2.0; // Rosenberg penalty weight
+
+    let (v1, n1) = lit_to_var(lits[0]);
+    let (v2, n2) = lit_to_var(lits[1]);
+    let (v3, n3) = lit_to_var(lits[2]);
+    let a = aux_var;
+
+    // We need to express yi = (1 - li) in terms of xi:
+    //   If literal is positive (li = xi): yi = 1 - xi
+    //   If literal is negated (li = 1 - xi): yi = xi
+    //
+    // So yi = xi if negated, yi = 1 - xi if positive.
+    //
+    // We compute the QUBO terms for:
+    //   H = a·y3 + M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a)
+    //
+    // Each term is expanded using yi = xi (if negated) or yi = 1-xi (if positive).
+
+    // Helper: add coefficient * yi * yj to the matrix
+    // where yi depends on variable vi and negation ni
+    let add_yy = |matrix: &mut [Vec<f64>], vi: usize, ni: bool, vj: usize, nj: bool, coeff: f64| {
+        // yi = xi if ni (negated literal), yi = 1 - xi if !ni (positive literal)
+        // yi * yj expansion:
+        if vi == vj {
+            // Same variable: yi * yj
+            // Both complemented the same way means yi = yj, so yi*yj = yi (binary)
+            // If ni == nj: yi*yj = yi^2 = yi (binary), add coeff * yi
+            // If ni != nj: yi * yj = xi * (1-xi) = 0 (always), add nothing
+            if ni == nj {
+                // yi * yi = yi (binary)
+                if ni {
+                    // yi = xi, add coeff * xi
+                    matrix[vi][vi] += coeff;
+                } else {
+                    // yi = 1 - xi, add coeff * (1 - xi) = coeff - coeff * xi
+                    // constant term ignored in QUBO (offset), diagonal:
+                    matrix[vi][vi] -= coeff;
+                }
+            }
+            // else: xi * (1-xi) = 0, nothing to add
+            return;
+        }
+        // Different variables: yi * yj
+        let (lo, hi, lo_neg, hi_neg) = if vi < vj {
+            (vi, vj, ni, nj)
+        } else {
+            (vj, vi, nj, ni)
+        };
+        // yi = xi if neg, else 1-xi
+        // yj = xj if neg, else 1-xj
+        // yi*yj = (xi if neg_i else 1-xi) * (xj if neg_j else 1-xj)
+        match (lo_neg, hi_neg) {
+            (true, true) => {
+                // xi * xj
+                matrix[lo][hi] += coeff;
+            }
+            (true, false) => {
+                // xi * (1 - xj) = xi - xi*xj
+                matrix[lo][lo] += coeff;
+                matrix[lo][hi] -= coeff;
+            }
+            (false, true) => {
+                // (1 - xi) * xj = xj - xi*xj
+                matrix[hi][hi] += coeff;
+                matrix[lo][hi] -= coeff;
+            }
+            (false, false) => {
+                // (1-xi)(1-xj) = 1 - xi - xj + xi*xj
+                // constant 1 ignored (offset)
+                matrix[lo][lo] -= coeff;
+                matrix[hi][hi] -= coeff;
+                matrix[lo][hi] += coeff;
+            }
+        }
+    };
+
+    // Helper: add coefficient * yi * a to the matrix
+    // where yi depends on variable vi and negation ni, a is aux variable
+    let add_ya = |matrix: &mut [Vec<f64>], vi: usize, ni: bool, a: usize, coeff: f64| {
+        // yi = xi if ni (negated literal), yi = 1-xi if !ni (positive literal)
+        // yi * a:
+        let (lo, hi) = if vi < a { (vi, a) } else { (a, vi) };
+        if ni {
+            // yi = xi, so yi * a = xi * a
+            matrix[lo][hi] += coeff;
+        } else {
+            // yi = 1 - xi, so yi * a = a - xi * a
+            matrix[a][a] += coeff;
+            matrix[lo][hi] -= coeff;
+        }
+    };
+
+    // Helper: add coefficient * yi to the matrix (linear term)
+    let add_y = |matrix: &mut [Vec<f64>], vi: usize, ni: bool, coeff: f64| {
+        if ni {
+            // yi = xi
+            matrix[vi][vi] += coeff;
+        } else {
+            // yi = 1 - xi, linear part: -coeff * xi (constant coeff ignored)
+            matrix[vi][vi] -= coeff;
+        }
+    };
+
+    // Term 1: a * y3 (coefficient = 1.0)
+    add_ya(matrix, v3, n3, a, 1.0);
+
+    // Term 2: M * y1 * y2
+    add_yy(matrix, v1, n1, v2, n2, penalty);
+
+    // Term 3: -2M * y1 * a
+    add_ya(matrix, v1, n1, a, -2.0 * penalty);
+
+    // Term 4: -2M * y2 * a
+    add_ya(matrix, v2, n2, a, -2.0 * penalty);
+
+    // Term 5: 3M * a (linear)
+    // a is a binary variable, a^2 = a, so linear a → diagonal
+    matrix[a][a] += 3.0 * penalty;
+
+    // We also need to add linear terms that come from constant offsets in products
+    // Actually, let's verify: the full expansion of
+    //   H = a·y3 + M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a)
+    // All terms are handled above.
+    //
+    // However, we need to account for the case where "add_ya" with !ni adds
+    // a linear term in `a`. Let me verify the add_ya logic handles this correctly.
+    //
+    // add_ya with ni=false: yi = 1-xi, yi*a = a - xi*a
+    //   matrix[a][a] += coeff (linear in a)
+    //   matrix[min(vi,a)][max(vi,a)] -= coeff (quadratic xi*a)
+    // This is correct.
+
+    // Note: We ignore constant terms (don't affect QUBO optimization).
+    let _ = add_y; // suppress unused warning - linear terms in y handled via products
+}
+
+/// Build a QUBO matrix from a KSatisfiability instance.
+///
+/// For K=2, directly encodes quadratic penalties.
+/// For K=3, uses Rosenberg quadratization with one auxiliary variable per clause.
+///
+/// Returns (matrix, num_source_vars) where matrix is (n + aux) x (n + aux).
+fn build_qubo_matrix(num_vars: usize, clauses: &[crate::models::satisfiability::CNFClause], k: usize) -> Vec<Vec<f64>> {
+    match k {
+        2 => {
+            let mut matrix = vec![vec![0.0; num_vars]; num_vars];
+            for clause in clauses {
+                add_2sat_clause_penalty(&mut matrix, &clause.literals);
+            }
+            matrix
+        }
+        3 => {
+            let num_aux = clauses.len(); // one auxiliary per clause
+            let total = num_vars + num_aux;
+            let mut matrix = vec![vec![0.0; total]; total];
+            for (idx, clause) in clauses.iter().enumerate() {
+                let aux_var = num_vars + idx;
+                add_3sat_clause_penalty(&mut matrix, &clause.literals, aux_var);
+            }
+            matrix
+        }
+        _ => unimplemented!("KSatisfiability to QUBO only supports K=2 and K=3"),
+    }
+}
+
 #[reduction(
     overhead = { ReductionOverhead::new(vec![("num_vars", poly!(num_vars))]) }
 )]
@@ -53,52 +316,31 @@ impl ReduceTo<QUBO<f64>> for KSatisfiability<2, i32> {
 
     fn reduce_to(&self) -> Self::Result {
         let n = self.num_vars();
-        let mut matrix = vec![vec![0.0; n]; n];
-
-        for clause in self.clauses() {
-            let lits = &clause.literals;
-            assert_eq!(lits.len(), 2, "Expected 2-SAT clause");
-
-            // Convert 1-indexed signed literals to 0-indexed (var, negated) pairs
-            let var_i = (lits[0].unsigned_abs() as usize) - 1;
-            let neg_i = lits[0] < 0;
-            let var_j = (lits[1].unsigned_abs() as usize) - 1;
-            let neg_j = lits[1] < 0;
-
-            let (i, j, ni, nj) = if var_i <= var_j {
-                (var_i, var_j, neg_i, neg_j)
-            } else {
-                (var_j, var_i, neg_j, neg_i)
-            };
-
-            // Penalty for unsatisfied clause: minimize penalty
-            // (l_i ∨ l_j) unsatisfied when both literals false
-            match (ni, nj) {
-                (false, false) => {
-                    // (x_i ∨ x_j): penalty = (1-x_i)(1-x_j) = 1 - x_i - x_j + x_i·x_j
-                    matrix[i][i] -= 1.0;
-                    matrix[j][j] -= 1.0;
-                    matrix[i][j] += 1.0;
-                }
-                (true, false) => {
-                    // (¬x_i ∨ x_j): penalty = x_i(1-x_j) = x_i - x_i·x_j
-                    matrix[i][i] += 1.0;
-                    matrix[i][j] -= 1.0;
-                }
-                (false, true) => {
-                    // (x_i ∨ ¬x_j): penalty = (1-x_i)x_j = x_j - x_i·x_j
-                    matrix[j][j] += 1.0;
-                    matrix[i][j] -= 1.0;
-                }
-                (true, true) => {
-                    // (¬x_i ∨ ¬x_j): penalty = x_i·x_j
-                    matrix[i][j] += 1.0;
-                }
-            }
-        }
+        let matrix = build_qubo_matrix(n, self.clauses(), 2);
 
         ReductionKSatToQUBO {
             target: QUBO::from_matrix(matrix),
+            source_num_vars: n,
+            source_size: self.problem_size(),
+        }
+    }
+}
+
+#[reduction(
+    overhead = { ReductionOverhead::new(vec![
+        ("num_vars", poly!(num_vars) + poly!(num_clauses)),
+    ]) }
+)]
+impl ReduceTo<QUBO<f64>> for KSatisfiability<3, i32> {
+    type Result = Reduction3SATToQUBO;
+
+    fn reduce_to(&self) -> Self::Result {
+        let n = self.num_vars();
+        let matrix = build_qubo_matrix(n, self.clauses(), 3);
+
+        Reduction3SATToQUBO {
+            target: QUBO::from_matrix(matrix),
+            source_num_vars: n,
             source_size: self.problem_size(),
         }
     }
diff --git a/src/rules/mod.rs b/src/rules/mod.rs
index 19662875..447215ee 100644
--- a/src/rules/mod.rs
+++ b/src/rules/mod.rs
@@ -62,7 +62,7 @@ pub use graph::{
 pub use coloring_qubo::ReductionKColoringToQUBO;
 pub use maximumindependentset_qubo::ReductionISToQUBO;
 pub use maximumindependentset_maximumsetpacking::{ReductionISToSP, ReductionSPToIS};
-pub use ksatisfiability_qubo::ReductionKSatToQUBO;
+pub use ksatisfiability_qubo::{ReductionKSatToQUBO, Reduction3SATToQUBO};
 pub use maximummatching_maximumsetpacking::ReductionMatchingToSP;
 pub use sat_coloring::ReductionSATToColoring;
 pub use maximumsetpacking_qubo::ReductionSPToQUBO;
diff --git a/src/unit_tests/rules/ksatisfiability_qubo.rs b/src/unit_tests/rules/ksatisfiability_qubo.rs
index 641e9216..328622fe 100644
--- a/src/unit_tests/rules/ksatisfiability_qubo.rs
+++ b/src/unit_tests/rules/ksatisfiability_qubo.rs
@@ -102,3 +102,83 @@ fn test_ksatisfiability_to_qubo_sizes() {
     assert!(!source_size.components.is_empty());
     assert!(!target_size.components.is_empty());
 }
+
+#[test]
+fn test_k3satisfiability_to_qubo_closed_loop() {
+    // 3-SAT: 5 vars, 7 clauses
+    let ksat = KSatisfiability::<3, i32>::new(
+        5,
+        vec![
+            CNFClause::new(vec![1, 2, -3]),   // x1 ∨ x2 ∨ ¬x3
+            CNFClause::new(vec![-1, 3, 4]),   // ¬x1 ∨ x3 ∨ x4
+            CNFClause::new(vec![2, -4, 5]),   // x2 ∨ ¬x4 ∨ x5
+            CNFClause::new(vec![-2, 3, -5]),  // ¬x2 ∨ x3 ∨ ¬x5
+            CNFClause::new(vec![1, -3, 5]),   // x1 ∨ ¬x3 ∨ x5
+            CNFClause::new(vec![-1, -2, 4]),  // ¬x1 ∨ ¬x2 ∨ x4
+            CNFClause::new(vec![3, -4, -5]),  // x3 ∨ ¬x4 ∨ ¬x5
+        ],
+    );
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&ksat);
+    let qubo = reduction.target_problem();
+
+    // QUBO should have 5 + 7 = 12 variables
+    assert_eq!(qubo.num_variables(), 12);
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_best(qubo);
+
+    // Verify all extracted solutions maximize satisfied clauses
+    for sol in &qubo_solutions {
+        let extracted = reduction.extract_solution(sol);
+        assert_eq!(extracted.len(), 5);
+        let satisfied = ksat.solution_size(&extracted).size;
+        assert_eq!(satisfied, 7, "Expected all 7 clauses satisfied");
+    }
+}
+
+#[test]
+fn test_k3satisfiability_to_qubo_single_clause() {
+    // Single 3-SAT clause: (x1 ∨ x2 ∨ x3) — 7 satisfying assignments
+    let ksat = KSatisfiability::<3, i32>::new(
+        3,
+        vec![CNFClause::new(vec![1, 2, 3])],
+    );
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&ksat);
+    let qubo = reduction.target_problem();
+
+    // 3 vars + 1 auxiliary = 4 total
+    assert_eq!(qubo.num_variables(), 4);
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_best(qubo);
+
+    // All solutions should satisfy the single clause
+    for sol in &qubo_solutions {
+        let extracted = reduction.extract_solution(sol);
+        assert_eq!(extracted.len(), 3);
+        assert!(ksat.solution_size(&extracted).is_valid);
+    }
+    // 7 out of 8 assignments satisfy (x1 ∨ x2 ∨ x3)
+    assert_eq!(qubo_solutions.len(), 7);
+}
+
+#[test]
+fn test_k3satisfiability_to_qubo_all_negated() {
+    // All negated: (¬x1 ∨ ¬x2 ∨ ¬x3) — 7 satisfying assignments
+    let ksat = KSatisfiability::<3, i32>::new(
+        3,
+        vec![CNFClause::new(vec![-1, -2, -3])],
+    );
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&ksat);
+    let qubo = reduction.target_problem();
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_best(qubo);
+
+    for sol in &qubo_solutions {
+        let extracted = reduction.extract_solution(sol);
+        assert!(ksat.solution_size(&extracted).is_valid);
+    }
+    // 7 out of 8 assignments satisfy (¬x1 ∨ ¬x2 ∨ ¬x3)
+    assert_eq!(qubo_solutions.len(), 7);
+}