problemreductions/rules/
factoring_ilp.rs

1//! Reduction from Factoring to ILP (Integer Linear Programming).
2//!
3//! The Integer Factoring problem can be formulated as an ILP using
4//! McCormick linearization for binary products combined with carry propagation.
5//!
6//! Given target N and bit widths m, n, find factors p (m bits) and q (n bits)
7//! such that p × q = N.
8//!
9//! ## Variables
10//! - `p_i ∈ {0,1}` for i = 0..m-1 (first factor bits)
11//! - `q_j ∈ {0,1}` for j = 0..n-1 (second factor bits)
12//! - `z_ij ∈ {0,1}` for each (i,j) pair (product p_i × q_j)
13//! - `c_k ∈ ℤ≥0` for k = 0..m+n-1 (carry at each bit position)
14//!
15//! ## Constraints
16//! 1. Product linearization (McCormick): z_ij ≤ p_i, z_ij ≤ q_j, z_ij ≥ p_i + q_j - 1
17//! 2. Bit-position sums: Σ_{i+j=k} z_ij + c_{k-1} = N_k + 2·c_k
18//! 3. No overflow: c_{m+n-1} = 0
19//! 4. Binary bounds: p_i ≤ 1, q_j ≤ 1
20//! 5. Carry bounds: 0 ≤ c_k ≤ min(m, n)
21
22use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
23use crate::models::misc::Factoring;
24use crate::reduction;
25use crate::rules::ilp_helpers::mccormick_product;
26use crate::rules::traits::{ReduceTo, ReductionResult};
27use std::cmp::min;
28
29/// Result of reducing Factoring to ILP.
30///
31/// This reduction creates an ILP where:
32/// - Binary variables represent factor bits and their products
33/// - Integer variables represent carries at each bit position
34/// - Constraints enforce the multiplication equals the target
35#[derive(Debug, Clone)]
36pub struct ReductionFactoringToILP {
37    target: ILP<i64>,
38    m: usize, // bits for first factor
39    n: usize, // bits for second factor
40}
41
42impl ReductionFactoringToILP {
43    /// Get the variable index for p_i (first factor bit i).
44    fn p_var(&self, i: usize) -> usize {
45        i
46    }
47
48    /// Get the variable index for q_j (second factor bit j).
49    fn q_var(&self, j: usize) -> usize {
50        self.m + j
51    }
52
53    /// Get the variable index for z_ij (product p_i × q_j).
54    #[cfg(test)]
55    fn z_var(&self, i: usize, j: usize) -> usize {
56        self.m + self.n + i * self.n + j
57    }
58
59    /// Get the variable index for carry at position k.
60    #[cfg(test)]
61    fn carry_var(&self, k: usize) -> usize {
62        self.m + self.n + self.m * self.n + k
63    }
64}
65
66impl ReductionResult for ReductionFactoringToILP {
67    type Source = Factoring;
68    type Target = ILP<i64>;
69
70    fn target_problem(&self) -> &ILP<i64> {
71        &self.target
72    }
73
74    /// Extract solution from ILP back to Factoring.
75    ///
76    /// The first m variables are p_i (first factor bits).
77    /// The next n variables are q_j (second factor bits).
78    /// Returns the decoded factors in ascending order.
79    fn extract_solution(
80        &self,
81        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
82    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
83        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
84
85        Ok({
86            // Extract p bits (first factor)
87            let p = (0..self.m)
88                .filter(|&i| target_solution[self.p_var(i)] == 1)
89                .fold(num_bigint::BigUint::from(0u8), |value, index| {
90                    value + (num_bigint::BigUint::from(1u8) << index)
91                });
92
93            // Extract q bits (second factor)
94            let q = (0..self.n)
95                .filter(|&j| target_solution[self.q_var(j)] == 1)
96                .fold(num_bigint::BigUint::from(0u8), |value, index| {
97                    value + (num_bigint::BigUint::from(1u8) << index)
98                });
99            if p <= q {
100                (p, q)
101            } else {
102                (q, p)
103            }
104        })
105    }
106}
107
108#[reduction(transform = upper_bound {
109    num_vars = "num_bits_first * num_bits_second + 2 * num_bits_first + 2 * num_bits_second + target_bits",
110    num_constraints = "3 * num_bits_first * num_bits_second + 4 * num_bits_first + 4 * num_bits_second + 3 * target_bits + 1",
111},
112    unavailable = {
113        num_nonzeros = "the exact target parameter is not represented by this reduction's symbolic transform",
114    }
115)]
116impl ReduceTo<ILP<i64>> for Factoring {
117    type Result = ReductionFactoringToILP;
118
119    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
120        let m = self.m();
121        let n = self.n();
122        let target = self.target();
123
124        // Calculate the number of bits needed for the target
125        let target_bits = self.target_bits();
126
127        // Number of bit positions to check: max(m+n, target_bits)
128        // For feasible instances, target_bits <= m+n (product of m-bit × n-bit has at most m+n bits).
129        // When target_bits > m+n, the ILP will be infeasible (target too large for given bit widths).
130        // Using max() here ensures proper infeasibility detection through the bit equations.
131        let num_bit_positions = std::cmp::max(m + n, target_bits);
132
133        // Total variables: m + n + m*n + num_bit_positions
134        let num_p = m;
135        let num_q = n;
136        let num_z = m * n;
137        let num_carries = num_bit_positions;
138        let num_vars = num_p + num_q + num_z + num_carries;
139
140        // Helper functions for variable indices
141        let p_var = |i: usize| -> usize { i };
142        let q_var = |j: usize| -> usize { m + j };
143        let z_var = |i: usize, j: usize| -> usize { m + n + i * n + j };
144        let carry_var = |k: usize| -> usize { m + n + m * n + k };
145
146        let mut constraints = Vec::new();
147
148        // Constraint 1: Product linearization (McCormick constraints)
149        // For each z_ij = p_i × q_j:
150        //   z_ij ≤ p_i
151        //   z_ij ≤ q_j
152        //   z_ij ≥ p_i + q_j - 1
153        for i in 0..m {
154            for j in 0..n {
155                let z = z_var(i, j);
156                let p = p_var(i);
157                let q = q_var(j);
158
159                constraints.extend(mccormick_product(z, p, q));
160            }
161        }
162
163        // Constraint 2: Bit-position equations
164        // For each bit position k = 0..num_bit_positions-1:
165        //   Σ_{i+j=k} z_ij + c_{k-1} = N_k + 2·c_k
166        // Rearranged: Σ_{i+j=k} z_ij + c_{k-1} - 2·c_k = N_k
167        for k in 0..num_bit_positions {
168            let mut terms: Vec<(usize, i64)> = Vec::new();
169
170            // Collect all z_ij where i + j = k
171            for i in 0..m {
172                if k >= i && k - i < n {
173                    let j = k - i;
174                    terms.push((z_var(i, j), 1));
175                }
176            }
177
178            // Add carry_in (from position k-1)
179            if k > 0 {
180                terms.push((carry_var(k - 1), 1));
181            }
182
183            // Subtract 2 × carry_out
184            terms.push((carry_var(k), -2));
185
186            // RHS is N_k (k-th bit of target).
187            let n_k = i64::from(target.bit(u64::try_from(k).expect("bit index fits u64")));
188            constraints.push(LinearConstraint::eq(terms, n_k));
189        }
190
191        // Constraint 3: Final carry must be zero (no overflow)
192        constraints.push(LinearConstraint::eq(
193            vec![(carry_var(num_bit_positions - 1), 1)],
194            0,
195        ));
196
197        // Constraint 4: Binary bounds for p_i and q_j (enforce 0/1 in integer domain)
198        for i in 0..m {
199            constraints.push(LinearConstraint::le(vec![(p_var(i), 1)], 1));
200        }
201        for j in 0..n {
202            constraints.push(LinearConstraint::le(vec![(q_var(j), 1)], 1));
203        }
204
205        // Constraint 5: Carry bounds (0 ≤ c_k ≤ min(m, n))
206        let carry_upper =
207            <Self as ReduceTo<ILP<i64>>>::exact_i64(min(m, n), "encoding a carry bound")?;
208        for k in 0..num_carries {
209            let cv = carry_var(k);
210            constraints.push(LinearConstraint::ge(vec![(cv, 1)], 0));
211            constraints.push(LinearConstraint::le(vec![(cv, 1)], carry_upper));
212        }
213
214        // Objective: feasibility problem (minimize 0)
215        let objective: Vec<(usize, i64)> = vec![];
216
217        let ilp = ILP::<i64>::new(num_vars, constraints, objective, ObjectiveSense::Minimize)
218            .map_err(<Self as ReduceTo<ILP<i64>>>::target_construction)?;
219
220        Ok(ReductionFactoringToILP { target: ilp, m, n })
221    }
222}
223
224#[cfg(feature = "example-db")]
225pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
226    vec![crate::example_db::specs::RuleExampleSpec {
227        id: "factoring_to_ilp",
228        build: || {
229            let source = Factoring::with_factor_bits(35, 3, 3);
230            crate::example_db::specs::rule_example_via_ilp::<_, i64>(source)
231        },
232    }]
233}
234
235#[cfg(test)]
236#[path = "../unit_tests/rules/factoring_ilp.rs"]
237mod tests;
problemreductions/rules/factoring_ilp.rs

problemreductions/rules/
factoring_ilp.rs
