problemreductions/rules/

factoring_ilp.rs

//! Reduction from Factoring to ILP (Integer Linear Programming).
//!
//! The Integer Factoring problem can be formulated as an 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
//! 4. Binary bounds: p_i ≤ 1, q_j ≤ 1
//! 5. Carry bounds: 0 ≤ c_k ≤ min(m, n)

use crate::models::algebraic::{Bounded, IntegerVariable, LinearConstraint, ObjectiveSense, ILP};
use crate::models::misc::Factoring;
use crate::reduction;
use crate::rules::ilp_helpers::mccormick_product;
use crate::rules::traits::{ReduceTo, ReductionResult};
use std::cmp::min;

/// 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<i64, i64, Bounded>,
    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<i64, i64, Bounded>;

    fn target_problem(&self) -> &ILP<i64, i64, Bounded> {
        &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 the decoded factors in ascending order.
    fn extract_solution(
        &self,
        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
        crate::rules::traits::validate_target_witness(
            self.target_problem(),
            target_solution,
            |value| value.value.is_some(),
            "target ILP assignment is infeasible",
        )?;

        Ok({
            // Extract p bits (first factor)
            let p = (0..self.m)
                .filter(|&i| target_solution[self.p_var(i)] == 1)
                .fold(num_bigint::BigUint::from(0u8), |value, index| {
                    value + (num_bigint::BigUint::from(1u8) << index)
                });

            // Extract q bits (second factor)
            let q = (0..self.n)
                .filter(|&j| target_solution[self.q_var(j)] == 1)
                .fold(num_bigint::BigUint::from(0u8), |value, index| {
                    value + (num_bigint::BigUint::from(1u8) << index)
                });
            if p <= q {
                (p, q)
            } else {
                (q, p)
            }
        })
    }
}

#[crate::aggregate_reduction(ilp_feasibility)]
impl crate::rules::AggregateReductionResult for ReductionFactoringToILP {}

// With L=max(m+n,target_bits), McCormick rows use 7mn terms; bit equations,
// final carry, factor bounds and carry bounds use mn+m+n+4L more.
// Substitute L <= m+n+target_bits to obtain a bound using source parameters.
#[reduction(transform = upper_bound {
    max_constraint_magnitude_bits = "num_bits_first + num_bits_second + 2",
    num_vars = "num_bits_first * num_bits_second + 2 * num_bits_first + 2 * num_bits_second + target_bits",
    num_constraints = "3 * num_bits_first * num_bits_second + 4 * num_bits_first + 4 * num_bits_second + 3 * target_bits + 1",
    num_nonzeros = "8 * num_bits_first * num_bits_second + 5 * num_bits_first + 5 * num_bits_second + 4 * target_bits",
})]
impl ReduceTo<ILP<i64, i64, Bounded>> for Factoring {
    type Result = ReductionFactoringToILP;

    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
        let m = self.m();
        let n = self.n();
        let target = self.target();

        // Calculate the number of bits needed for the target
        let target_bits = self.target_bits();

        // 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 };

        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);

                constraints.extend(mccormick_product(z, p, q));
            }
        }

        // 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, i64)> = 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));
                }
            }

            // Add carry_in (from position k-1)
            if k > 0 {
                terms.push((carry_var(k - 1), 1));
            }

            // Subtract 2 × carry_out
            terms.push((carry_var(k), -2));

            // RHS is N_k (k-th bit of target).
            let n_k = i64::from(target.bit(u64::try_from(k).expect("bit index fits u64")));
            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,
        ));

        // Constraint 4: Binary bounds for p_i and q_j (enforce 0/1 in integer domain)
        for i in 0..m {
            constraints.push(LinearConstraint::le(vec![(p_var(i), 1)], 1));
        }
        for j in 0..n {
            constraints.push(LinearConstraint::le(vec![(q_var(j), 1)], 1));
        }

        // Constraint 5: Carry bounds (0 ≤ c_k ≤ min(m, n))
        let carry_upper = <Self as ReduceTo<ILP<i64, i64, Bounded>>>::exact_i64(
            min(m, n),
            "encoding a carry bound",
        )?;
        for k in 0..num_carries {
            let cv = carry_var(k);
            constraints.push(LinearConstraint::ge(vec![(cv, 1)], 0));
            constraints.push(LinearConstraint::le(vec![(cv, 1)], carry_upper));
        }

        // Objective: feasibility problem (minimize 0)
        let objective: Vec<(usize, i64)> = vec![];

        let mut variables = vec![IntegerVariable::binary(); num_vars];
        variables[num_p + num_q + num_z..].fill(
            IntegerVariable::new(Some(0), Some(carry_upper))
                .map_err(<Self as ReduceTo<ILP<i64, i64, Bounded>>>::target_construction)?,
        );

        let ilp = ILP::<i64, i64, Bounded>::with_variables(
            variables,
            constraints,
            objective,
            ObjectiveSense::Minimize,
        )
        .map_err(<Self as ReduceTo<ILP<i64, i64, Bounded>>>::target_construction)?;

        Ok(ReductionFactoringToILP { target: ilp, m, n })
    }
}

#[cfg(feature = "example-db")]
pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
    vec![crate::example_db::specs::RuleExampleSpec {
        id: "factoring_to_ilp",
        build: || {
            let source = Factoring::with_factor_bits(35, 3, 3);
            crate::example_db::specs::rule_example_via_bounded_ilp::<_>(source)
        },
    }]
}

#[cfg(test)]
#[path = "../unit_tests/rules/factoring_ilp.rs"]
mod tests;