diff --git a/docs/paper/reduction_graph.json b/docs/paper/reduction_graph.json index d9a6c40a..cf26ef10 100644 --- a/docs/paper/reduction_graph.json +++ b/docs/paper/reduction_graph.json @@ -1,8 +1,13 @@ { "nodes": [ { - "id": "MaxCut", - "label": "MaxCut", + "id": "CircuitSAT", + "label": "CircuitSAT", + "category": "satisfiability" + }, + { + "id": "Coloring", + "label": "Coloring", "category": "graph" }, { @@ -11,13 +16,13 @@ "category": "graph" }, { - "id": "SpinGlass", - "label": "SpinGlass", - "category": "optimization" + "id": "Factoring", + "label": "Factoring", + "category": "specialized" }, { - "id": "QUBO", - "label": "QUBO", + "id": "ILP", + "label": "ILP", "category": "optimization" }, { @@ -26,49 +31,49 @@ "category": "graph" }, { - "id": "SetPacking", - "label": "SetPacking", - "category": "set" - }, - { - "id": "Satisfiability", - "label": "Satisfiability", + "id": "KSatisfiability", + "label": "KSatisfiability", "category": "satisfiability" }, { - "id": "Coloring", - "label": "Coloring", + "id": "Matching", + "label": "Matching", "category": "graph" }, { - "id": "VertexCovering", - "label": "VertexCovering", + "id": "MaxCut", + "label": "MaxCut", "category": "graph" }, { - "id": "KSatisfiability", - "label": "KSatisfiability", - "category": "satisfiability" + "id": "QUBO", + "label": "QUBO", + "category": "optimization" }, { - "id": "CircuitSAT", - "label": "CircuitSAT", + "id": "Satisfiability", + "label": "Satisfiability", "category": "satisfiability" }, { - "id": "Factoring", - "label": "Factoring", - "category": "specialized" + "id": "SetCovering", + "label": "SetCovering", + "category": "set" }, { - "id": "Matching", - "label": "Matching", - "category": "graph" + "id": "SetPacking", + "label": "SetPacking", + "category": "set" }, { - "id": "SetCovering", - "label": "SetCovering", - "category": "set" + "id": "SpinGlass", + "label": "SpinGlass", + "category": "optimization" + }, + { + "id": "VertexCovering", + "label": "VertexCovering", + "category": "graph" } ], "edges": [ @@ -78,38 +83,33 @@ "bidirectional": false }, { - "source": "Matching", - "target": "SetPacking", + "source": "Factoring", + "target": "CircuitSAT", "bidirectional": false }, - { - "source": "Satisfiability", - "target": "KSatisfiability", - "bidirectional": true - }, { "source": "Factoring", - "target": "CircuitSAT", + "target": "ILP", "bidirectional": false }, { - "source": "IndependentSet", - "target": "SetPacking", + "source": "KSatisfiability", + "target": "Satisfiability", "bidirectional": true }, { - "source": "SpinGlass", - "target": "QUBO", - "bidirectional": true + "source": "Matching", + "target": "SetPacking", + "bidirectional": false }, { - "source": "MaxCut", - "target": "SpinGlass", - "bidirectional": true + "source": "Satisfiability", + "target": "Coloring", + "bidirectional": false }, { - "source": "VertexCovering", - "target": "SetCovering", + "source": "Satisfiability", + "target": "DominatingSet", "bidirectional": false }, { @@ -118,18 +118,28 @@ "bidirectional": false }, { - "source": "Satisfiability", - "target": "Coloring", - "bidirectional": false + "source": "SetPacking", + "target": "IndependentSet", + "bidirectional": true + }, + { + "source": "SpinGlass", + "target": "MaxCut", + "bidirectional": true }, { - "source": "IndependentSet", - "target": "VertexCovering", + "source": "SpinGlass", + "target": "QUBO", "bidirectional": true }, { - "source": "Satisfiability", - "target": "DominatingSet", + "source": "VertexCovering", + "target": "IndependentSet", + "bidirectional": true + }, + { + "source": "VertexCovering", + "target": "SetCovering", "bidirectional": false } ] diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index ccd6cdb5..601289af 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -64,6 +64,7 @@ "VertexCovering": (-0.5, 2), "Matching": (-2, 2), "SpinGlass": (2.5, 2), + "ILP": (3.5, 1), // Row 3: Leaf nodes "SetPacking": (-1.5, 3), "SetCovering": (0.5, 3), @@ -330,6 +331,53 @@ In all graph problems below, $G = (V, E)$ denotes an undirected graph with $|V| _Solution extraction._ Without ancilla: identity. With ancilla: if $sigma_a = 1$, flip all spins before removing ancilla. ] +#theorem[ + *(Factoring $arrow.r$ ILP)* Integer factorization reduces to binary ILP using McCormick linearization with $O(m n)$ variables and constraints. +] + +#proof[ + _Construction._ For target $N$ with $m$-bit factor $p$ and $n$-bit factor $q$: + + _Variables:_ Binary $p_i, q_j in {0,1}$ for factor bits; binary $z_(i j) in {0,1}$ for products $p_i dot q_j$; integer $c_k >= 0$ for carries at each bit position. + + _Product linearization (McCormick):_ For each $z_(i j) = p_i dot q_j$: + $ z_(i j) <= p_i, quad z_(i j) <= q_j, quad z_(i j) >= p_i + q_j - 1 $ + + _Bit-position equations:_ For each bit position $k$: + $ sum_(i+j=k) z_(i j) + c_(k-1) = N_k + 2 c_k $ + where $N_k$ is the $k$-th bit of $N$ and $c_(-1) = 0$. + + _No overflow:_ $c_(m+n-1) = 0$. + + _Correctness._ The McCormick constraints enforce $z_(i j) = p_i dot q_j$ for binary variables. The bit equations encode $p times q = N$ via carry propagation, matching array multiplier semantics. + + _Solution extraction._ Read $p = sum_i p_i 2^i$ and $q = sum_j q_j 2^j$ from the binary variables. +] + +_Example: Factoring 15._ The following Rust code demonstrates the closed-loop reduction (requires `ilp` feature: `cargo add problemreductions --features ilp`): + +```rust +use problemreductions::prelude::*; + +// 1. Create factoring instance: find p (4-bit) × q (4-bit) = 15 +let problem = Factoring::new(4, 4, 15); + +// 2. Reduce to ILP +let reduction = ReduceTo::::reduce_to(&problem); +let ilp = reduction.target_problem(); + +// 3. Solve ILP +let solver = ILPSolver::new(); +let ilp_solution = solver.solve(ilp).unwrap(); + +// 4. Extract factoring solution +let extracted = reduction.extract_solution(&ilp_solution); + +// 5. Verify: reads factors and confirms p × q = 15 +let (p, q) = problem.read_factors(&extracted); +assert_eq!(p * q, 15); // e.g., (3, 5) or (5, 3) +``` + = Summary #let gray = rgb("#e8e8e8") @@ -352,6 +400,7 @@ In all graph problems below, $G = (V, E)$ denotes an undirected graph with $|V| [CircuitSAT $arrow.r$ SpinGlass], [$O(|"gates"|)$], [@whitfield2012 @lucas2014], [Factoring $arrow.r$ CircuitSAT], [$O(m n)$], [Folklore], [SpinGlass $arrow.l.r$ MaxCut], [$O(n + |J|)$], [@barahona1982 @lucas2014], + table.cell(fill: gray)[Factoring $arrow.r$ ILP], table.cell(fill: gray)[$O(m n)$], table.cell(fill: gray)[—], ), caption: [Summary of reductions. Gray rows indicate trivial reductions.] ) diff --git a/docs/plans/2026-01-31-factoring-ilp-design.md b/docs/plans/2026-01-31-factoring-ilp-design.md new file mode 100644 index 00000000..2c049628 --- /dev/null +++ b/docs/plans/2026-01-31-factoring-ilp-design.md @@ -0,0 +1,169 @@ +# Factoring to ILP Reduction Design + +**Issue:** #21 - Implement integer programming based Factoring problem solver +**Date:** 2026-01-31 +**Status:** Approved + +## Overview + +Implement a reduction from the Factoring problem to Integer Linear Programming (ILP). Given a target number N and bit widths m and n, find binary factors p and q such that p × q = N. + +The key challenge is linearizing the multiplication constraint, which we solve using McCormick relaxation for binary products combined with carry-propagation to handle bit-position sums. + +## ILP Formulation + +### Variables + +| Variable | Range | Count | Description | +|----------|-------|-------|-------------| +| `p_i` | {0,1} | m | First factor bits (i = 0..m-1) | +| `q_j` | {0,1} | n | Second factor bits (j = 0..n-1) | +| `z_ij` | {0,1} | m×n | Product terms `p_i × q_j` | +| `c_k` | [0, min(m,n)] | m+n | Carry at bit position k | + +**Total variables:** m + n + m×n + (m+n) = O(m×n) + +### Constraints + +**1. Product linearization** (for each i ∈ [0,m), j ∈ [0,n)): +``` +z_ij ≤ p_i +z_ij ≤ q_j +z_ij ≥ p_i + q_j - 1 +``` +These McCormick constraints ensure z_ij = p_i × q_j for binary variables. + +**2. Bit-position equations** (for each k ∈ [0, m+n)): +``` +Σ_{i+j=k} z_ij + c_{k-1} = N_k + 2·c_k +``` +Where: +- N_k is the k-th bit of target N +- c_{-1} = 0 (no incoming carry at position 0) + +**3. No overflow:** +``` +c_{m+n-1} = 0 +``` + +**Total constraints:** 3×m×n + (m+n) + 1 = O(m×n) + +### Objective + +Feasibility problem: minimize 0. + +## Implementation + +### File Structure + +``` +src/rules/factoring_ilp.rs +``` + +### Data Structures + +```rust +pub struct ReductionFactoringToILP { + target: ILP, + source_size: ProblemSize, + m: usize, // bits for first factor + n: usize, // bits for second factor +} +``` + +### Variable Layout + +Contiguous indices in ILP: +``` +[0, m) → p_i (first factor bits) +[m, m+n) → q_j (second factor bits) +[m+n, m+n+m*n) → z_ij (products, row-major order) +[m+n+m*n, 2(m+n)+m*n) → c_k (carries) +``` + +### Helper Methods + +```rust +impl ReductionFactoringToILP { + fn p_var(&self, i: usize) -> usize { i } + fn q_var(&self, j: usize) -> usize { self.m + j } + fn z_var(&self, i: usize, j: usize) -> usize { + self.m + self.n + i * self.n + j + } + fn carry_var(&self, k: usize) -> usize { + self.m + self.n + self.m * self.n + k + } +} +``` + +### Solution Extraction + +1. Read p_i values from indices [0, m) → first factor bits +2. Read q_j values from indices [m, m+n) → second factor bits +3. Return concatenated bit vector matching Factoring problem format + +## Testing Strategy + +### Correctness Tests + +| Test | Target | Bits | Expected Factors | +|------|--------|------|------------------| +| factor_6 | 6 | 2,2 | 2×3 or 3×2 | +| factor_15 | 15 | 3,3 | 3×5, 5×3, 1×15, 15×1 | +| factor_35 | 35 | 3,3 | 5×7 or 7×5 | + +### Edge Cases + +| Test | Description | +|------|-------------| +| factor_one | 1 = 1×1 | +| factor_prime | 7 = 1×7 or 7×1 | +| factor_square | 9 = 3×3 | +| infeasible | Target 100 with 2-bit factors (max 9) | + +### Validation Tests + +- Compare ILP solutions with brute force solver +- Verify solution extraction produces valid factorization +- Check variable and constraint counts match formulas + +### Closed-Loop Pattern (Issue #3) + +```rust +let problem = Factoring::new(m, n, target); +let reduction = ReduceTo::::reduce_to(&problem); +let ilp = reduction.target_problem(); +let ilp_solution = ILPSolver::new().solve(ilp)?; +let extracted = reduction.extract_solution(&ilp_solution); +assert!(problem.is_valid_factorization(&extracted)); +``` + +## Documentation Updates + +Per issue #3 requirements: + +1. Add Factoring → ILP to `docs/paper/reductions.typ`: + - Theorem with proof sketch + - Add to summary table with overhead O(m×n) + +2. Update reduction graph: + - Register in inventory + - Add ILP node (already exists from Coloring PR) + - Add edge Factoring → ILP + +3. Regenerate `docs/paper/reduction_graph.json` + +## Complexity Analysis + +| Metric | Value | +|--------|-------| +| Variables | m + n + m×n + (m+n) | +| Constraints | 3×m×n + (m+n) + 1 | +| Time complexity | O(m×n) for reduction | +| Space complexity | O(m×n) | + +## References + +- McCormick, G. P. (1976). Computability of global solutions to factorable nonconvex programs. Mathematical Programming. +- Issue #3: Coding rules for AI agents +- Existing implementation: `src/rules/vertexcovering_ilp.rs` (pattern reference) diff --git a/docs/src/reductions/graph.md b/docs/src/reductions/graph.md index a737470e..00044cef 100644 --- a/docs/src/reductions/graph.md +++ b/docs/src/reductions/graph.md @@ -37,10 +37,12 @@ flowchart TD SG_f64["SpinGlass<f64>"] QUBO["QUBO<f64>"] MaxCut["MaxCut<i32>"] + ILP["ILP"] end %% Factoring chain Factoring --> CircuitSAT + Factoring --> ILP CircuitSAT --> SG_i32 %% SAT reductions @@ -75,6 +77,7 @@ flowchart TD style SG_f64 fill:#ff9,stroke:#333 style QUBO fill:#ff9,stroke:#333 style MaxCut fill:#ff9,stroke:#333 + style ILP fill:#ff9,stroke:#333 ``` ## Legend @@ -126,6 +129,7 @@ println!("Types: {}, Reductions: {}", graph.num_types(), graph.num_reductions()) | Satisfiability | DominatingSet | No | | CircuitSAT | SpinGlass<i32> | No | | Factoring | CircuitSAT | No | +| Factoring | ILP | No | ## API diff --git a/src/rules/factoring_ilp.rs b/src/rules/factoring_ilp.rs new file mode 100644 index 00000000..055cb95f --- /dev/null +++ b/src/rules/factoring_ilp.rs @@ -0,0 +1,587 @@ +//! Reduction from Factoring to ILP (Integer Linear Programming). +//! +//! The Integer Factoring problem can be formulated as a binary ILP using +//! McCormick linearization for binary products combined with carry propagation. +//! +//! Given target N and bit widths m, n, find factors p (m bits) and q (n bits) +//! such that p × q = N. +//! +//! ## Variables +//! - `p_i ∈ {0,1}` for i = 0..m-1 (first factor bits) +//! - `q_j ∈ {0,1}` for j = 0..n-1 (second factor bits) +//! - `z_ij ∈ {0,1}` for each (i,j) pair (product p_i × q_j) +//! - `c_k ∈ ℤ≥0` for k = 0..m+n-1 (carry at each bit position) +//! +//! ## Constraints +//! 1. Product linearization (McCormick): z_ij ≤ p_i, z_ij ≤ q_j, z_ij ≥ p_i + q_j - 1 +//! 2. Bit-position sums: Σ_{i+j=k} z_ij + c_{k-1} = N_k + 2·c_k +//! 3. No overflow: c_{m+n-1} = 0 + +use crate::models::optimization::{ILP, LinearConstraint, ObjectiveSense, VarBounds}; +use crate::models::specialized::Factoring; +use crate::polynomial::{Monomial, Polynomial}; +use crate::rules::registry::{ReductionEntry, ReductionOverhead}; +use crate::rules::traits::{ReduceTo, ReductionResult}; +use crate::traits::Problem; +use crate::types::ProblemSize; +use std::cmp::min; + +// Register reduction in the inventory for automatic discovery +inventory::submit! { + ReductionEntry { + source_name: "Factoring", + target_name: "ILP", + source_graph: "Factoring", + target_graph: "ILPMatrix", + overhead_fn: || ReductionOverhead::new(vec![ + // num_vars = m + n + m*n + num_carries where num_carries = max(m+n, target_bits) + // For feasible instances, target_bits <= m+n, so this is 2(m+n) + m*n + ("num_vars", Polynomial { + terms: vec![ + Monomial::var("num_bits_first").scale(2.0), + Monomial::var("num_bits_second").scale(2.0), + Monomial { + coefficient: 1.0, + variables: vec![("num_bits_first", 1), ("num_bits_second", 1)], + }, + ] + }), + // num_constraints = 3*m*n + num_bit_positions + 1 + // For feasible instances (target_bits <= m+n), this is 3*m*n + (m+n) + 1 + ("num_constraints", Polynomial { + terms: vec![ + Monomial { + coefficient: 3.0, + variables: vec![("num_bits_first", 1), ("num_bits_second", 1)], + }, + Monomial::var("num_bits_first"), + Monomial::var("num_bits_second"), + Monomial::constant(1.0), + ] + }), + ]), + } +} + +/// Result of reducing Factoring to ILP. +/// +/// This reduction creates an ILP where: +/// - Binary variables represent factor bits and their products +/// - Integer variables represent carries at each bit position +/// - Constraints enforce the multiplication equals the target +#[derive(Debug, Clone)] +pub struct ReductionFactoringToILP { + target: ILP, + source_size: ProblemSize, + m: usize, // bits for first factor + n: usize, // bits for second factor +} + +impl ReductionFactoringToILP { + /// Get the variable index for p_i (first factor bit i). + fn p_var(&self, i: usize) -> usize { + i + } + + /// Get the variable index for q_j (second factor bit j). + fn q_var(&self, j: usize) -> usize { + self.m + j + } + + /// Get the variable index for z_ij (product p_i × q_j). + #[cfg(test)] + fn z_var(&self, i: usize, j: usize) -> usize { + self.m + self.n + i * self.n + j + } + + /// Get the variable index for carry at position k. + #[cfg(test)] + fn carry_var(&self, k: usize) -> usize { + self.m + self.n + self.m * self.n + k + } +} + +impl ReductionResult for ReductionFactoringToILP { + type Source = Factoring; + type Target = ILP; + + fn target_problem(&self) -> &ILP { + &self.target + } + + /// Extract solution from ILP back to Factoring. + /// + /// The first m variables are p_i (first factor bits). + /// The next n variables are q_j (second factor bits). + /// Returns concatenated bit vector [p_0, ..., p_{m-1}, q_0, ..., q_{n-1}]. + fn extract_solution(&self, target_solution: &[usize]) -> Vec { + // Extract p bits (first factor) + let p_bits: Vec = (0..self.m) + .map(|i| target_solution.get(self.p_var(i)).copied().unwrap_or(0)) + .collect(); + + // Extract q bits (second factor) + let q_bits: Vec = (0..self.n) + .map(|j| target_solution.get(self.q_var(j)).copied().unwrap_or(0)) + .collect(); + + // Concatenate p and q bits + let mut result = p_bits; + result.extend(q_bits); + result + } + + fn source_size(&self) -> ProblemSize { + self.source_size.clone() + } + + fn target_size(&self) -> ProblemSize { + self.target.problem_size() + } +} + +impl ReduceTo for Factoring { + type Result = ReductionFactoringToILP; + + fn reduce_to(&self) -> Self::Result { + let m = self.m(); + let n = self.n(); + let target = self.target(); + + // Calculate the number of bits needed for the target + let target_bits = if target == 0 { + 1 + } else { + (64 - target.leading_zeros()) as usize + }; + + // Number of bit positions to check: max(m+n, target_bits) + // For feasible instances, target_bits <= m+n (product of m-bit × n-bit has at most m+n bits). + // When target_bits > m+n, the ILP will be infeasible (target too large for given bit widths). + // Using max() here ensures proper infeasibility detection through the bit equations. + let num_bit_positions = std::cmp::max(m + n, target_bits); + + // Total variables: m + n + m*n + num_bit_positions + let num_p = m; + let num_q = n; + let num_z = m * n; + let num_carries = num_bit_positions; + let num_vars = num_p + num_q + num_z + num_carries; + + // Helper functions for variable indices + let p_var = |i: usize| -> usize { i }; + let q_var = |j: usize| -> usize { m + j }; + let z_var = |i: usize, j: usize| -> usize { m + n + i * n + j }; + let carry_var = |k: usize| -> usize { m + n + m * n + k }; + + // Variable bounds + let mut bounds = Vec::with_capacity(num_vars); + + // p_i, q_j, z_ij are binary + for _ in 0..(num_p + num_q + num_z) { + bounds.push(VarBounds::binary()); + } + + // c_k are non-negative integers with upper bound min(m, n) + // (at most min(m, n) products can contribute to any position) + let carry_upper = min(m, n) as i64; + for _ in 0..num_carries { + bounds.push(VarBounds::bounded(0, carry_upper)); + } + + let mut constraints = Vec::new(); + + // Constraint 1: Product linearization (McCormick constraints) + // For each z_ij = p_i × q_j: + // z_ij ≤ p_i + // z_ij ≤ q_j + // z_ij ≥ p_i + q_j - 1 + for i in 0..m { + for j in 0..n { + let z = z_var(i, j); + let p = p_var(i); + let q = q_var(j); + + // z_ij - p_i ≤ 0 + constraints.push(LinearConstraint::le( + vec![(z, 1.0), (p, -1.0)], + 0.0, + )); + + // z_ij - q_j ≤ 0 + constraints.push(LinearConstraint::le( + vec![(z, 1.0), (q, -1.0)], + 0.0, + )); + + // z_ij - p_i - q_j ≥ -1 + constraints.push(LinearConstraint::ge( + vec![(z, 1.0), (p, -1.0), (q, -1.0)], + -1.0, + )); + } + } + + // Constraint 2: Bit-position equations + // For each bit position k = 0..num_bit_positions-1: + // Σ_{i+j=k} z_ij + c_{k-1} = N_k + 2·c_k + // Rearranged: Σ_{i+j=k} z_ij + c_{k-1} - 2·c_k = N_k + for k in 0..num_bit_positions { + let mut terms: Vec<(usize, f64)> = Vec::new(); + + // Collect all z_ij where i + j = k + for i in 0..m { + if k >= i && k - i < n { + let j = k - i; + terms.push((z_var(i, j), 1.0)); + } + } + + // Add carry_in (from position k-1) + if k > 0 { + terms.push((carry_var(k - 1), 1.0)); + } + + // Subtract 2 × carry_out + terms.push((carry_var(k), -2.0)); + + // RHS is N_k (k-th bit of target). For k >= 64, the bit is 0 for u64. + let n_k = if k < 64 { + ((target >> k) & 1) as f64 + } else { + 0.0 + }; + constraints.push(LinearConstraint::eq(terms, n_k)); + } + + // Constraint 3: Final carry must be zero (no overflow) + constraints.push(LinearConstraint::eq( + vec![(carry_var(num_bit_positions - 1), 1.0)], + 0.0, + )); + + // Objective: feasibility problem (minimize 0) + let objective: Vec<(usize, f64)> = vec![]; + + let ilp = ILP::new( + num_vars, + bounds, + constraints, + objective, + ObjectiveSense::Minimize, + ); + + ReductionFactoringToILP { + target: ilp, + source_size: self.problem_size(), + m, + n, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::solvers::{BruteForce, ILPSolver, Solver}; + + #[test] + fn test_reduction_creates_valid_ilp() { + // Factor 6 with 2-bit factors + let problem = Factoring::new(2, 2, 6); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + // Check variable count: m + n + m*n + (m+n) = 2 + 2 + 4 + 4 = 12 + assert_eq!(ilp.num_vars, 12); + + // Check constraint count: 3*m*n + (m+n) + 1 = 12 + 4 + 1 = 17 + assert_eq!(ilp.constraints.len(), 17); + + assert_eq!(ilp.sense, ObjectiveSense::Minimize); + } + + #[test] + fn test_variable_layout() { + let problem = Factoring::new(3, 2, 6); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + + // p variables: [0, 1, 2] + assert_eq!(reduction.p_var(0), 0); + assert_eq!(reduction.p_var(2), 2); + + // q variables: [3, 4] + assert_eq!(reduction.q_var(0), 3); + assert_eq!(reduction.q_var(1), 4); + + // z variables: [5, 6, 7, 8, 9, 10] (3x2 = 6) + assert_eq!(reduction.z_var(0, 0), 5); + assert_eq!(reduction.z_var(0, 1), 6); + assert_eq!(reduction.z_var(1, 0), 7); + assert_eq!(reduction.z_var(2, 1), 10); + + // carry variables: [11, 12, 13, 14, 15] (m+n = 5) + assert_eq!(reduction.carry_var(0), 11); + assert_eq!(reduction.carry_var(4), 15); + } + + #[test] + fn test_factor_6() { + // 6 = 2 × 3 or 3 × 2 + let problem = Factoring::new(2, 2, 6); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + // Verify it's a valid factorization + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 6); + } + + #[test] + fn test_factor_15() { + // Closed-loop test for factoring 15 = 3 × 5 (or 5 × 3, 1 × 15, 15 × 1) + + // 1. Create factoring instance: find p (4-bit) × q (4-bit) = 15 + let problem = Factoring::new(4, 4, 15); + + // 2. Reduce to ILP + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + // 3. Solve ILP + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + + // 4. Extract factoring solution + let extracted = reduction.extract_solution(&ilp_solution); + + // 5. Verify: solution is valid and p × q = 15 + assert!(problem.is_valid_factorization(&extracted)); + let (p, q) = problem.read_factors(&extracted); + assert_eq!(p * q, 15); // e.g., (3, 5) or (5, 3) + } + + #[test] + fn test_factor_35() { + // 35 = 5 × 7 or 7 × 5 + let problem = Factoring::new(3, 3, 35); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 35); + } + + #[test] + fn test_factor_one() { + // 1 = 1 × 1 + let problem = Factoring::new(2, 2, 1); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 1); + } + + #[test] + fn test_factor_prime() { + // 7 is prime: 7 = 1 × 7 or 7 × 1 + let problem = Factoring::new(3, 3, 7); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 7); + } + + #[test] + fn test_factor_square() { + // 9 = 3 × 3 + let problem = Factoring::new(3, 3, 9); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 9); + } + + #[test] + fn test_infeasible_target_too_large() { + // Target 100 with 2-bit factors (max product is 3 × 3 = 9) + let problem = Factoring::new(2, 2, 100); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let result = ilp_solver.solve(ilp); + + assert!(result.is_none(), "Should be infeasible"); + } + + #[test] + fn test_ilp_matches_brute_force() { + let problem = Factoring::new(2, 2, 6); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + // Get ILP solution + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let ilp_factors = reduction.extract_solution(&ilp_solution); + + // Get brute force solutions + let bf = BruteForce::new(); + let bf_solutions = bf.find_best(&problem); + + // ILP solution should be among brute force solutions + let (a, b) = problem.read_factors(&ilp_factors); + let bf_pairs: Vec<(u64, u64)> = bf_solutions + .iter() + .map(|s| problem.read_factors(s)) + .collect(); + + assert!( + bf_pairs.contains(&(a, b)), + "ILP solution ({}, {}) should be in brute force solutions {:?}", + a, + b, + bf_pairs + ); + } + + #[test] + fn test_solution_extraction() { + let problem = Factoring::new(2, 2, 6); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + + // Manually construct ILP solution for 2 × 3 = 6 + // p = 2 = binary 10 -> p_0=0, p_1=1 + // q = 3 = binary 11 -> q_0=1, q_1=1 + // z_00 = p_0 * q_0 = 0, z_01 = p_0 * q_1 = 0 + // z_10 = p_1 * q_0 = 1, z_11 = p_1 * q_1 = 1 + // Variables: [p0, p1, q0, q1, z00, z01, z10, z11, c0, c1, c2, c3] + let ilp_solution = vec![0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0]; + let extracted = reduction.extract_solution(&ilp_solution); + + // Should extract [p0, p1, q0, q1] = [0, 1, 1, 1] + assert_eq!(extracted, vec![0, 1, 1, 1]); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a, 2); + assert_eq!(b, 3); + assert_eq!(a * b, 6); + } + + #[test] + fn test_source_and_target_size() { + let problem = Factoring::new(3, 4, 12); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + + let source_size = reduction.source_size(); + let target_size = reduction.target_size(); + + assert_eq!(source_size.get("num_bits_first"), Some(3)); + assert_eq!(source_size.get("num_bits_second"), Some(4)); + + // num_vars = 3 + 4 + 12 + 7 = 26 + assert_eq!(target_size.get("num_vars"), Some(26)); + + // num_constraints = 3*12 + 7 + 1 = 44 + assert_eq!(target_size.get("num_constraints"), Some(44)); + } + + #[test] + fn test_solve_reduced() { + let problem = Factoring::new(2, 2, 6); + + let ilp_solver = ILPSolver::new(); + let solution = ilp_solver + .solve_reduced(&problem) + .expect("solve_reduced should work"); + + assert!(problem.is_valid_factorization(&solution)); + } + + #[test] + fn test_asymmetric_bit_widths() { + // 12 = 3 × 4 or 4 × 3 or 2 × 6 or 6 × 2 or 1 × 12 or 12 × 1 + let problem = Factoring::new(2, 4, 12); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let ilp_solver = ILPSolver::new(); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert!(problem.is_valid_factorization(&extracted)); + + let (a, b) = problem.read_factors(&extracted); + assert_eq!(a * b, 12); + } + + #[test] + fn test_constraint_count_formula() { + // Verify constraint count matches formula: 3*m*n + (m+n) + 1 + for (m, n) in [(2, 2), (3, 3), (2, 4), (4, 2)] { + let problem = Factoring::new(m, n, 1); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let expected = 3 * m * n + (m + n) + 1; + assert_eq!( + ilp.constraints.len(), + expected, + "Constraint count mismatch for m={}, n={}", + m, + n + ); + } + } + + #[test] + fn test_variable_count_formula() { + // Verify variable count matches formula: m + n + m*n + (m+n) + for (m, n) in [(2, 2), (3, 3), (2, 4), (4, 2)] { + let problem = Factoring::new(m, n, 1); + let reduction: ReductionFactoringToILP = ReduceTo::::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let expected = m + n + m * n + (m + n); + assert_eq!( + ilp.num_vars, expected, + "Variable count mismatch for m={}, n={}", + m, n + ); + } + } +} diff --git a/src/rules/graph.rs b/src/rules/graph.rs index d96c53e0..2cb09ce9 100644 --- a/src/rules/graph.rs +++ b/src/rules/graph.rs @@ -162,9 +162,6 @@ impl ReductionGraph { } } - // Also register manual reductions for backward compatibility - Self::register_reductions(&mut graph, &name_indices); - Self { graph, name_indices, @@ -260,6 +257,7 @@ impl ReductionGraph { SpinGlass => "SpinGlass", SpinGlass => "SpinGlass", QUBO => "QUBO", + ILP => "ILP", // Satisfiability problems Satisfiability => "Satisfiability", KSatisfiability<3, i32> => "KSatisfiability", @@ -269,57 +267,6 @@ impl ReductionGraph { } } - fn register_reductions( - graph: &mut DiGraph<&'static str, ReductionEdge>, - name_indices: &HashMap<&'static str, NodeIndex>, - ) { - // Add an edge between two problem types by name (with default overhead). - // This is for backward compatibility with manually registered reductions. - macro_rules! add_edge { - ($src:expr => $dst:expr) => { - if let (Some(&src), Some(&dst)) = (name_indices.get($src), name_indices.get($dst)) { - // Avoid duplicate edges - if graph.find_edge(src, dst).is_none() { - graph.add_edge( - src, - dst, - ReductionEdge { - source_graph: "SimpleGraph", - target_graph: "SimpleGraph", - overhead: ReductionOverhead::default(), - }, - ); - } - } - }; - } - - // Graph problem reductions - add_edge!("IndependentSet" => "VertexCovering"); - add_edge!("VertexCovering" => "IndependentSet"); - add_edge!("IndependentSet" => "SetPacking"); - add_edge!("SetPacking" => "IndependentSet"); - add_edge!("VertexCovering" => "SetCovering"); - add_edge!("Matching" => "SetPacking"); - - // Optimization reductions - add_edge!("SpinGlass" => "QUBO"); - add_edge!("QUBO" => "SpinGlass"); - add_edge!("MaxCut" => "SpinGlass"); - add_edge!("SpinGlass" => "MaxCut"); - - // SAT-based reductions - add_edge!("Satisfiability" => "KSatisfiability"); - add_edge!("KSatisfiability" => "Satisfiability"); - add_edge!("Satisfiability" => "IndependentSet"); - add_edge!("Satisfiability" => "Coloring"); - add_edge!("Satisfiability" => "DominatingSet"); - - // Circuit reductions - add_edge!("CircuitSAT" => "SpinGlass"); - add_edge!("Factoring" => "CircuitSAT"); - } - /// Check if `sub` is a subtype of `sup` (or equal). pub fn is_graph_subtype(&self, sub: &str, sup: &str) -> bool { sub == sup @@ -571,8 +518,8 @@ impl ReductionGraph { } } - // Build nodes with categories - let nodes: Vec = self + // Build nodes with categories, sorted by id for deterministic output + let mut nodes: Vec = self .name_indices .keys() .map(|&name| { @@ -584,9 +531,11 @@ impl ReductionGraph { } }) .collect(); + nodes.sort_by(|a, b| a.id.cmp(&b.id)); // Build edges (only include one direction for bidirectional edges) - let edges: Vec = edge_set + // Sort by (source, target) for deterministic output + let mut edges: Vec = edge_set .into_iter() .map(|((src, dst), bidirectional)| EdgeJson { source: src.to_string(), @@ -594,6 +543,7 @@ impl ReductionGraph { bidirectional, }) .collect(); + edges.sort_by(|a, b| (&a.source, &a.target).cmp(&(&b.source, &b.target))); ReductionGraphJson { nodes, edges } } @@ -624,7 +574,7 @@ impl ReductionGraph { "graph" } else if name.contains("SetPacking") || name.contains("SetCover") { "set" - } else if name.contains("SpinGlass") || name.contains("QUBO") { + } else if name.contains("SpinGlass") || name.contains("QUBO") || name.contains("ILP") { "optimization" } else if name.contains("Satisfiability") || name.contains("SAT") { "satisfiability" diff --git a/src/rules/mod.rs b/src/rules/mod.rs index ac502adc..fecbfd01 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -34,6 +34,8 @@ mod setcovering_ilp; mod setpacking_ilp; #[cfg(feature = "ilp")] mod vertexcovering_ilp; +#[cfg(feature = "ilp")] +mod factoring_ilp; pub use graph::{EdgeJson, NodeJson, ReductionEdge, ReductionGraph, ReductionGraphJson, ReductionPath}; pub use traits::{ReduceTo, ReductionResult}; @@ -67,3 +69,5 @@ pub use setcovering_ilp::ReductionSCToILP; pub use setpacking_ilp::ReductionSPToILP; #[cfg(feature = "ilp")] pub use vertexcovering_ilp::ReductionVCToILP; +#[cfg(feature = "ilp")] +pub use factoring_ilp::ReductionFactoringToILP;