Constraint propagation (constraint programming)

Deep dive

How it works

Parts:

• domain store — one u64 bitset per variable (values 0..63) plus the trail of (var, old_domain)
• propagators — one per constraint; run deletes values, false on an empty domain
• watch lists — variable → propagators to wake (here on any change; solvers split bounds/fixed events)
• queue — propagators waiting, each at most once
• search — DFS, first-fail branching, trail undo

1. Post constraints: queue each propagator once, register it on its variables.
2. Pop a propagator and run it; each narrowing is trailed and queues the other watchers.
3. Empty domain: clear the queue, fail. Queue empty: fixpoint.
4. Pick the unfixed variable with the smallest domain (first-fail); none left = solution.
5. For each value: note the trail length, set the domain, propagate, recurse.
6. Undo the trail to the noted length, try the next value.

 queue -> propagator -> narrow: trail push, wake watchers -> queue
 queue empty = fixpoint -> pick var, try value -> queue
 empty domain -> undo trail -> next value

Consistency levels. Arc consistency: every value of x has a supporting value in y; AC-3 revises arcs from a queue (Mackworth 1977). Bounds consistency: only min and max must be supported (cheap for x + y = z, holes remain). Domain consistency: every remaining value fits some solution of that constraint.

Global constraints see many variables at once. The all-different below only removes fixed values and counts pigeonholes; Régin's matching filter reaches domain consistency in polynomial time.

Trailing is an undo journal as in journaled-state-revert: backtracking costs the changes made, not the store size. Two watched literals avoid even that.

Branching. First-fail takes the smallest domain; dom/wdeg divides domain size by the failure-weighted constraint count of the variable.

Lazy clause generation (Ohrimenko, Stuckey, Codish 2009) records an explanation clause per propagated change; a failure then yields a learned clause as in CDCL.

CP vs MIP vs SAT. CP: rich constraints, strong on feasibility and scheduling, weak on linear cost bounds. MIP: LP bounds, strong on cost. SAT: plain clauses, fast learning, large encodings. CP-SAT combines CP, SAT and LP.

Code

use std::time::Instant;
enum Prop { Ne(usize, usize, i32), AllDiff(Vec<usize>), Sum(usize, usize, usize) } // x != y+c; alldiff; x+y=z
struct Store {
    dom: Vec<u64>,            // bit v set = value v still possible (values 0..63)
    trail: Vec<(usize, u64)>, // (var, old domain), pushed before every change
    props: Vec<Prop>,
    watch: Vec<Vec<usize>>,   // var -> propagators woken when its domain changes
    queue: Vec<usize>, queued: Vec<bool>, nodes: u64,
}
fn bits(lo: i32, hi: i32) -> u64 { // mask of values lo..=hi clipped to 0..63
    let (lo, hi) = (lo.max(0), hi.min(63));
    if lo > hi { 0 } else { (u64::MAX >> (63 - (hi - lo) as u32)) << lo as u32 }
}
fn min(d: u64) -> i32 { d.trailing_zeros() as i32 }
fn max(d: u64) -> i32 { 63 - d.leading_zeros() as i32 }
impl Store {
    fn new(dom: Vec<u64>) -> Store {
        Store { watch: vec![vec![]; dom.len()], dom, trail: vec![], props: vec![], queue: vec![], queued: vec![], nodes: 0 }
    }
    fn post(&mut self, p: Prop) {
        let vars = match &p { Prop::Ne(x, y, _) => vec![*x, *y], Prop::AllDiff(v) => v.clone(), Prop::Sum(x, y, z) => vec![*x, *y, *z] };
        for v in vars { self.watch[v].push(self.props.len()); }
        self.queue.push(self.props.len());
        self.queued.push(true); self.props.push(p);
    }
    /// Narrow a domain, trail the old value, wake watchers; false = wipe-out.
    fn set(&mut self, v: usize, new: u64) -> bool {
        if new != self.dom[v] {
            self.trail.push((v, self.dom[v]));
            self.dom[v] = new;
            for &p in &self.watch[v] { if !self.queued[p] { self.queued[p] = true; self.queue.push(p); } }
        }
        new != 0
    }
    fn undo(&mut self, mark: usize) { while self.trail.len() > mark { let (v, d) = self.trail.pop().unwrap(); self.dom[v] = d; } }
    fn fixpoint(&mut self) -> bool {
        while let Some(p) = self.queue.pop() {
            self.queued[p] = false;
            if !self.run(p) { for q in self.queue.drain(..) { self.queued[q] = false; } return false; }
        }
        true
    }
    fn run(&mut self, p: usize) -> bool {
        match self.props[p] {
            Prop::Ne(x, y, c) => {
                let (dx, dy) = (self.dom[x], self.dom[y]); // a fixed side removes its value from the other
                (dy.count_ones() != 1 || self.set(x, dx & !bits(min(dy) + c, min(dy) + c)))
                    && (dx.count_ones() != 1 || self.set(y, self.dom[y] & !bits(min(dx) - c, min(dx) - c)))
            }
            Prop::AllDiff(ref vs) => {
                let (vs, mut union) = (vs.clone(), 0u64);
                for &v in &vs {
                    let d = self.dom[v];
                    union |= d;
                    if d.count_ones() == 1 {
                        for &w in &vs { if w != v && !self.set(w, self.dom[w] & !d) { return false; } }
                    }
                }
                union.count_ones() as usize >= vs.len() // pigeonhole
            }
            Prop::Sum(x, y, z) => { // bounds consistency for x + y = z
                let (a, b) = (self.dom[x], self.dom[y]);
                let c = self.dom[z] & bits(min(a) + min(b), max(a) + max(b));
                self.set(z, c)
                    && self.set(x, a & bits(min(c) - max(b), max(c) - min(b)))
                    && self.set(y, b & bits(min(c) - max(self.dom[x]), max(c) - min(self.dom[x])))
            }
        }
    }
    fn search(&mut self, leaf: &mut dyn FnMut(&[u64])) {
        let pick = (0..self.dom.len()).filter(|&v| self.dom[v].count_ones() > 1).min_by_key(|&v| self.dom[v].count_ones());
        let Some(v) = pick else { return leaf(&self.dom) };
        let d = self.dom[v];
        for val in (0..64).filter(|b| d >> b & 1 == 1) {
            self.nodes += 1;
            let mark = self.trail.len();
            if self.set(v, 1 << val) && self.fixpoint() { self.search(leaf); }
            self.undo(mark);
        }
    }
}
fn queens(n: usize) -> Store { // var = column, value = row
    let mut s = Store::new(vec![bits(0, n as i32 - 1); n]);
    s.post(Prop::AllDiff((0..n).collect()));
    for i in 0..n { for j in i + 1..n { // diagonals: row_i != row_j +- (j-i)
        s.post(Prop::Ne(i, j, (j - i) as i32)); s.post(Prop::Ne(i, j, i as i32 - j as i32));
    } }
    s
}
fn plain(n: usize, rows: &mut Vec<usize>, nodes: &mut u64) -> u64 {
    if rows.len() == n { return 1; }
    let col = rows.len();
    (0..n).map(|r| {
        *nodes += 1;
        if !rows.iter().enumerate().all(|(c, &q)| q != r && q.abs_diff(r) != col - c) { return 0; }
        rows.push(r); let f = plain(n, rows, nodes); rows.pop(); f
    }).sum()
}
fn main() {
    for (n, expect) in [(8usize, 92u64), (10, 724)] {
        let mut s = queens(n);
        assert!(s.fixpoint());
        let (root, mut count, t) = (s.dom.clone(), 0, Instant::now());
        s.search(&mut |d| {
            let r: Vec<i32> = d.iter().map(|&x| min(x)).collect();
            assert!((0..n).all(|i| (i + 1..n).all(|j| r[i] != r[j] && (r[i] - r[j]).abs() != (j - i) as i32)));
            count += 1;
        }); let t1 = t.elapsed();
        assert!(count == expect && s.dom == root); // trail restored the root
        let (mut nodes, t) = (0, Instant::now());
        assert_eq!(plain(n, &mut vec![], &mut nodes), expect);
        println!("n={n}: {count} solutions | propagation {} nodes {t1:?} | plain {nodes} nodes {:?}", s.nodes, t.elapsed());
        assert!(s.nodes < nodes);
    }
    let mut s = Store::new(vec![bits(2, 5), bits(3, 4), bits(0, 20)]); // x + y = z
    s.post(Prop::Sum(0, 1, 2));
    assert!(s.fixpoint() && s.dom[2] == bits(5, 9));
    assert!(s.set(2, bits(8, 8)) && s.fixpoint() && s.dom[0] == bits(4, 5));
}

Compiled and run (Rust stable, release). Asserts: counts 92 and 724 (OEIS A000170), every solution re-checked, root domains restored exactly after the search, fewer nodes with propagation. "Nodes" are value assignments attempted.

Complexity and constants

• Propagator run: O(arity) for Ne and Sum, O(arity^2) for this all-different (it also clones its list per run).
• Measured on my Mac (noisy, about 10 runs): n=8, 662 nodes in 0.43-0.95 ms vs 15,720 plain nodes in 0.10-0.30 ms; n=10, 10,014 nodes in 8-18 ms vs 348,150 plain in 2.2-8 ms. Plain won in wall time every run; per node about 1 µs vs about 10 ns (derived, order of).

Efficiency

• Cost: per node, the triggered propagator runs plus one trail push per narrowing; bound by memory access, not arithmetic.
• Levers:
  • Bitset domains — a removal is one AND, a singleton test one popcount (this page) — small domains.
  • Event filtering (bounds / fixed / any) — skips propagators a change cannot affect (not measured here) — expensive propagators.
  • Clause learning — avoids repeating failures (Ohrimenko et al. 2009) — hard instances.
• Naive → standard → state of the art: plain backtracking · 348,150 nodes at n=10, about 10 ns each; propagation with first-fail · 10,014 nodes, about 1 µs each (measured); matching globals plus learning, as in CP-SAT and Chuffed.

Implementation notes

• Route every change through set (16-byte trail entry); a direct write to dom breaks undo.
• A non-idempotent propagator must requeue itself.

Pitfalls

• queued flags left set after a failure: those propagators never run again.
• Arc consistent does not mean solvable.

Variants and alternatives

In production

CP-SAT, Gecode, Chuffed and Choco-solver are general solvers; MiniZinc is a solver-independent modelling language. Fields: timetable planning, constrained-random hardware stimulus, tile synthesis.

Implementations

Sources

• Alan K. Mackworth, Consistency in networks of relations, Artificial Intelligence 8(1), 99-118, 1977. paper
• Olga Ohrimenko, Peter J. Stuckey, Michael Codish, Propagation via lazy clause generation, Constraints 14, 357-391, 2009. paper
• Jean-Charles Régin, A filtering algorithm for constraints of difference in CSPs, AAAI 1994 (cited from memory, no link).
• N. J. A. Sloane (ed.), OEIS A000170, nonattacking queens counts. entry