Constraint Satisfaction

Why Constraint Satisfaction Matters

Suppose several exams need time slots, and exams sharing students cannot overlap. Choosing each slot independently fails because one decision changes what remains possible elsewhere. The same pattern appears in map coloring, scheduling, configuration, and Sudoku: variables have domains, constraints restrict combinations, and a solution must satisfy every constraint at once.

A constraint satisfaction problem, or CSP, makes that structure explicit. Instead of generating complete assignments and testing them afterward, a solver reduces domains as decisions are made. Forward checking removes conflicts immediately adjacent to an assignment. Propagation carries the consequences farther. MRV, the degree heuristic, and LCV choose where to search and which value to try first.

This lesson builds a finite-domain binary CSP solver with both forward-checking and maintaining-arc-consistency modes. You will trace an LCV decision, watch AC-3 propagate beyond the latest variable's neighbors, and follow a failed branch whose domain deletions must be completely undone before a later branch can succeed.

Core Concepts and Terminology

Write a finite CSP as variables $X_0,\ldots,X_{n-1}$, initial domains $D_0,\ldots,D_{n-1}$, and constraints. The implementation handles binary constraints: a relation $R_{ij}$ lists allowed pairs $(a,b)$ for variables $X_i,X_j$. Unary restrictions are represented by initial domains. Higher-arity constraints require another representation or an appropriate encoding.

The constraint graph has one vertex per variable and an edge when two variables share a relation. Its graph edge is undirected, but its relation can be asymmetric: $X_i\lt X_j$ does not mean $X_j\lt X_i$. When processing the reverse arc, swap the coordinates of every allowed pair. The programs construct that transpose automatically.

For an arc $X_i\to X_j$, a value $a\in D_i$ has support when some $b\in D_j$ satisfies $(a,b)\in R_{ij}$. The arc is consistent when every current value of $X_i$ has such support. A value with no support can never appear in a complete solution extending the current branch, so removing it is sound. If any domain becomes empty, the branch is impossible.

Arc consistency checks local supports, not all constraints simultaneously. A triangle of pairwise-different variables with domains $\{0,1\}$ is arc consistent: either value at either endpoint has a different neighbor value as support. Nevertheless, the triangle cannot be colored with two values.

The output is either SAT with one complete assignment or UNSAT after all required alternatives are excluded. This is a feasibility solver, not a counting, uniqueness, or optimization solver. No node budget is imposed; exhausting a separate budget would require an explicit incomplete result.

How Constraint Satisfaction Works, Step by Step

The solver maintains current domain masks, an assignment array using $-1$ for unassigned variables, and a trail of pairs (variable, old_mask). Every domain change records its previous mask before deletion. A checkpoint is just the current trail length. Restoring records in reverse order reconstructs the parent state even if propagation changed the same variable several times.

Choose a variable. Among unassigned variables, minimize current domain size. Among ties, maximize the number of unassigned graph neighbors. Break any remaining tie by variable index. Already assigned neighbors do not contribute to this dynamic degree.

Order its values. For candidate $a$ at variable $X$, let $S_{XY}(a)$ be the values of neighbor $Y$ allowed with $a$. The LCV score is

$$\operatorname{score}(a)=\sum_{Y\text{ unassigned neighbor of }X}\lvert D(Y)\setminus S_{XY}(a)\rvert.$$

Try smaller scores first, breaking ties by numeric value. This score predicts direct neighbor deletions; it does not simulate all future AC-3 consequences or prove that a candidate can lead to a solution.

Assign and forward-check. Save a checkpoint, mark $X=a$, and restrict its domain to $\{a\}$. Check constraints against assigned neighbors and intersect each unassigned neighbor's domain with its compatible-value mask. A wipeout rejects the attempt immediately.

Propagate when using MAC. Initialize an AC-3 queue with both directed arcs of every constraint edge. Revising $X\to Y$ removes every value of $X$ without support in $D(Y)$. If $D(X)$ changes, queue each incoming arc $Z\to X$ with $Z\ne Y$: those variables may have depended on a removed support. Stop on an empty domain or when the queue is empty.

The following map-coloring snapshot shows the difference. After assigning South Australia SA=0, its mainland neighbors have domains $\{1,2\}$. Assigning NT=1 directly forces WA=2 and Q=2. AC-3 then uses Q=2 to force NSW=1, which forces V=2; those last reductions are beyond NT's direct forward-checking step.

Recurse and undo. Enter a child only when filtering succeeds. At a complete assignment, save a copy of the answer. On every return, restore the trail to the checkpoint and set the selected variable back to unassigned, including when a solution was found. Root preprocessing also uses the trail, so the function can restore every original domain before returning.

flowchart TD
    A["Validate domains and relation tables"] --> B["In MAC mode, run root AC-3"]
    B --> C{"Any empty domain?"}
    C -->|Yes| D["Restore root trail; return UNSAT"]
    C -->|No| E["Enter DFS; count one node"]
    E --> F{"All variables assigned?"}
    F -->|Yes| G["Copy complete answer"]
    F -->|No| H["MRV, degree tie-break, then LCV order"]
    H --> I["Checkpoint; assign next value; forward-check"]
    I --> J{"Forward checking succeeds?"}
    J -->|No| L["Undo to checkpoint; unassign variable"]
    J -->|Yes| Q["In MAC mode, propagate with AC-3"]
    Q --> K{"Propagation succeeds?"}
    K -->|Yes| E
    K -->|No| L
    G --> L
    L --> M{"Solution saved?"}
    M -->|Yes| N["Unwind callers; restore root state; return SAT"]
    M -->|No| O{"More ordered values?"}
    O -->|Yes| I
    O -->|No| P["Return failure; caller restores its checkpoint"]

The flowchart summarizes recursive calls: after a child returns, its caller performs the checkpoint restoration. An empty AC-3 queue means all queued consistency requirements have been processed. It is not an assignment and is not, by itself, a proof of satisfiability.

Worked Example with a Full Trace

Use variables $A,B,C$, indexed $0,1,2$, and impose $A\ne B$, $A\ne C$, $B\ne C$. Let initial domains be $D(A)=\{0,1\}$, $D(B)=\{0,1,2\}$, and $D(C)=\{0,2\}$. With bit $2^a$ representing value $a$, their masks are $[3,7,5]$.

Root AC-3 revises all six arcs without deleting a value. MRV ties $A$ and $C$ at size $2$, and their unassigned degrees are both $2$, so index $0$ selects $A$. The two LCV scores differ:

The complete successful descent is:

The deletion counter includes restricting the chosen variable itself, so choosing $A=1$ removes one value from $A$ and one from $B$. Choosing $B=0$ removes one from $B$ and one from $C$. Choosing already-singleton $C=2$ deletes nothing. This implementation initializes all six arcs after each decision, even when no further pruning occurs.

Unwinding is part of the trace. Undoing the $C$ decision changes its assignment back to unassigned, but its domain remains $\{2\}$ because that restriction came from the parent. Undoing the $B$ decision restores $D(B)=D(C)=\{0,2\}$. Undoing $A$ restores all original masks $[3,7,5]$. The saved answer remains $[1,0,2]$, while historical work counters remain unchanged.

Now give all three triangle variables domains $\{0,1\}$. Root AC-3 still removes nothing. Trying $A=0$ forward-checks both neighbors to $\{1\}$; revising $B\to C$ then empties $D(B)$ because equality is forbidden. Undo and try $A=1$; the symmetric conflict follows. MAC reports UNSAT with $1$ DFS node, $2$ candidate attempts, $8$ cumulative value deletions, and $14$ arc revisions. Forward checking alone needs $3$ DFS nodes and $4$ attempts because it reaches the conflict after explicitly choosing the next variable.

Implementation in Python, C++, and Java

Each panel contains the same generic binary-relation solver and six sample runs. Variables use indices $0$ through $n-1$, with $0\le n\le12$. Values use $0$ through $d-1$, with $1\le d\le8$. Initial domains are bit masks; an empty mask is a valid unsatisfiable input. An empty variable set with no constraints is satisfiable and has an empty assignment.

Each constraint edge supplies its two variable indices and an allowed-pair list. An empty relation forbids all endpoint combinations. Self-constraints and repeated edges, including a repeated edge in reverse order, are rejected. Combine multiple restrictions on the same pair into their intersection before calling the solver. Repeated pairs inside one allowed-pair list are harmless duplicates. All constraints must be represented: the absence of an edge means there is no binary restriction between those variables.

FC mode uses forward checking after each assignment. MAC mode additionally runs root AC-3 and restarts AC-3 from all directed arcs after each successful forward-checking step. More selective queue initialization is possible, but the full queue makes this implementation and its counters easier to reproduce. Both modes use MRV, unassigned-neighbor degree, and direct-deletion LCV.

The returned counters are historical totals: nodes counts entered DFS states, attempts counts tried values, pruned counts deleted domain values including decision restrictions, and revisions counts dequeued AC-3 arcs, even when nothing changes. Repeated deletions in different branches are counted repeatedly; undo does not subtract work. A root filtering failure can return with zero DFS nodes.

The restore fixture has variables $A,B,C,D,E,F$. The last five begin with domain $\{0,1,2\}$, while $A$ begins with $\{0,1\}$. Variables $B,C,D$ form an inequality triangle. When $A=0$, its relations restrict each of them to $\{0,1\}$ and leave $E,F$ unrestricted. When $A=1$, they keep all three values but $E,F$ must be $0$. LCV scores are $3$ for $A=0$ and $4$ for $A=1$, so the first choice is globally impossible despite its smaller local score. Correct rollback lets the later choice find $[1,0,1,2,0,0]$.

The map fixture uses order [WA,NT,SA,Q,NSW,V,T] and the usual mainland adjacency edges listed directly in the program. Initial degree selects SA, which has five neighbors. Tasmania T is unconstrained. The code returns one valid coloring, not a geographically drawn map or a unique solution.

from collections import deque


def solve(initial, edges, d=3, use_ac3=True):
    n = len(initial)
    if not 0 <= n <= 12 or not 1 <= d <= 8:
        raise ValueError("invalid variable or value count")
    full = (1 << d) - 1
    if any(mask < 0 or mask & ~full for mask in initial):
        raise ValueError("invalid domain mask")
    neighbors = [[] for _ in range(n)]
    table = [[[0] * d for _ in range(n)] for _ in range(n)]
    seen = set()
    for x, y, pairs in edges:
        key = (min(x, y), max(x, y))
        if not 0 <= x < n or not 0 <= y < n or x == y or key in seen:
            raise ValueError("invalid or repeated constraint edge")
        seen.add(key)
        neighbors[x].append(y)
        neighbors[y].append(x)
        for a, b in pairs:
            if not 0 <= a < d or not 0 <= b < d:
                raise ValueError("invalid relation value")
            table[x][y][a] |= 1 << b
            table[y][x][b] |= 1 << a
    for row in neighbors:
        row.sort()
    domains, assigned, trail = list(initial), [-1] * n, []
    nodes = attempts = pruned = revisions = 0
    answer = []

    def shrink(x, new):
        nonlocal pruned
        old = domains[x]
        assert not new & ~old
        if new != old:
            trail.append((x, old))
            domains[x] = new
            pruned += old.bit_count() - new.bit_count()

    def undo(mark):
        while len(trail) > mark:
            x, old = trail.pop()
            domains[x] = old

    def ac3():
        nonlocal revisions
        queue = deque((x, y) for x in range(n) for y in neighbors[x])
        while queue:
            x, y = queue.popleft()
            revisions += 1
            old, new = domains[x], 0
            for a in range(d):
                if old & (1 << a) and table[x][y][a] & domains[y]:
                    new |= 1 << a
            if new != old:
                shrink(x, new)
                if new == 0:
                    return False
                for z in neighbors[x]:
                    if z != y:
                        queue.append((z, x))
        return True

    def forward(x, a):
        for y in neighbors[x]:
            allowed = table[x][y][a]
            if assigned[y] != -1:
                if not allowed & (1 << assigned[y]):
                    return False
            else:
                shrink(y, domains[y] & allowed)
                if domains[y] == 0:
                    return False
        return True

    def dfs():
        nonlocal nodes, attempts, answer
        nodes += 1
        open_vars = [x for x in range(n) if assigned[x] == -1]
        if not open_vars:
            answer = assigned.copy()
            return True
        x = min(open_vars, key=lambda v: (domains[v].bit_count(), -sum(assigned[y] == -1 for y in neighbors[v]), v))
        values = [a for a in range(d) if domains[x] & (1 << a)]

        def score(a):
            return sum((domains[y] & ~table[x][y][a]).bit_count() for y in neighbors[x] if assigned[y] == -1)

        values.sort(key=lambda a: (score(a), a))
        for a in values:
            attempts += 1
            mark = len(trail)
            assigned[x] = a
            shrink(x, 1 << a)
            ok = forward(x, a)
            if ok and use_ac3:
                ok = ac3()
            found = ok and dfs()
            undo(mark)
            assigned[x] = -1
            if found:
                return True
        return False

    ok = all(domains)
    if ok and use_ac3:
        ok = ac3()
    found = ok and dfs()
    undo(0)
    assert domains == list(initial) and assigned == [-1] * n and not trail
    return found, answer, nodes, attempts, pruned, revisions


def neq(x, y, d):
    return x, y, [(a, b) for a in range(d) for b in range(d) if a != b]


def fixtures():
    triangle = [neq(0, 1, 3), neq(0, 2, 3), neq(1, 2, 3)]
    triangle2 = [neq(0, 1, 2), neq(0, 2, 2), neq(1, 2, 2)]
    restore = [(0, x, [(a, b) for a in range(2) for b in range(3) if a == 1 or b < 2]) for x in [1, 2, 3]]
    restore += [(0, x, [(a, b) for a in range(2) for b in range(3) if a == 0 or b == 0]) for x in [4, 5]]
    restore += [neq(1, 2, 3), neq(1, 3, 3), neq(2, 3, 3)]
    australia = [neq(x, y, 3) for x, y in [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3), (2, 4), (2, 5), (3, 4), (4, 5)]]
    return [("asymmetric / FC", [3, 7, 5], triangle, 3, False),
            ("asymmetric / MAC", [3, 7, 5], triangle, 3, True),
            ("triangle2 / FC", [3, 3, 3], triangle2, 2, False),
            ("triangle2 / MAC", [3, 3, 3], triangle2, 2, True),
            ("restore / MAC", [3, 7, 7, 7, 7, 7], restore, 3, True),
            ("Australia / MAC", [7] * 7, australia, 3, True)]


def fmt(values):
    return "[" + ",".join(map(str, values)) + "]"


def main():
    ok = True
    for name, domains, edges, d, mac in fixtures():
        original = domains.copy()
        found, answer, nodes, attempts, pruned, revisions = solve(domains, edges, d, mac)
        print(f"{name}: {'SAT' if found else 'UNSAT'}, nodes={nodes}, attempts={attempts}, pruned={pruned}, revisions={revisions}")
        print("assignment = " + fmt(answer))
        valid = found == (not name.startswith("triangle2"))
        if found:
            valid = valid and len(answer) == len(domains) and all(domains[x] & (1 << a) for x, a in enumerate(answer))
            valid = valid and all((answer[x], answer[y]) in pairs for x, y, pairs in edges)
        ok = ok and valid and domains == original
    print("checks = " + ("PASS" if ok else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <deque>
#include <functional>
#include <iostream>
#include <stdexcept>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
using Count = std::int64_t;
struct Edge { int x, y; std::vector<std::pair<int,int>> pairs; };
struct Result { bool found; std::vector<int> answer; Count nodes, attempts, pruned, revisions; };
int popcount(int mask) { int count = 0; while (mask) { mask &= mask - 1; ++count; } return count; }
Result solve(const std::vector<int>& initial, const std::vector<Edge>& edges, int d = 3, bool use_ac3 = true) {
    int n = static_cast<int>(initial.size());
    if (initial.size() > 12 || d < 1 || d > 8) throw std::invalid_argument("invalid variable or value count");
    int full = (1 << d) - 1;
    for (int mask : initial) if (mask < 0 || (mask & ~full)) throw std::invalid_argument("invalid domain mask");
    std::vector<std::vector<int>> neighbors(n);
    std::vector<std::vector<bool>> seen(n, std::vector<bool>(n));
    std::vector<std::vector<std::vector<int>>> table(n, std::vector<std::vector<int>>(n, std::vector<int>(d)));
    for (const auto& edge : edges) {
        int x = edge.x, y = edge.y;
        if (x < 0 || x >= n || y < 0 || y >= n || x == y || seen[x][y]) throw std::invalid_argument("invalid or repeated edge");
        seen[x][y] = seen[y][x] = true; neighbors[x].push_back(y); neighbors[y].push_back(x);
        for (auto p : edge.pairs) {
            int a = p.first, b = p.second;
            if (a < 0 || a >= d || b < 0 || b >= d) throw std::invalid_argument("invalid relation value");
            table[x][y][a] |= 1 << b; table[y][x][b] |= 1 << a;
        }
    }
    for (auto& row : neighbors) std::sort(row.begin(),row.end());
    std::vector<int> domains = initial, assigned(n,-1), answer;
    std::vector<std::pair<int,int>> trail;
    Count nodes = 0, attempts = 0, pruned = 0, revisions = 0;
    auto shrink = [&](int x, int fresh) {
        int old = domains[x]; assert((fresh & ~old) == 0);
        if (fresh != old) { trail.emplace_back(x,old); domains[x] = fresh; pruned += popcount(old) - popcount(fresh); }
    };
    auto undo = [&](std::size_t mark) {
        while (trail.size() > mark) { auto p = trail.back(); trail.pop_back(); domains[p.first] = p.second; }
    };
    auto ac3 = [&]() {
        std::deque<std::pair<int,int>> queue;
        for (int x = 0; x < n; ++x) for (int y : neighbors[x]) queue.emplace_back(x,y);
        while (!queue.empty()) {
            auto p = queue.front(); queue.pop_front(); int x = p.first, y = p.second; ++revisions;
            int old = domains[x], fresh = 0;
            for (int a = 0; a < d; ++a) if ((old & (1 << a)) && (table[x][y][a] & domains[y])) fresh |= 1 << a;
            if (fresh != old) {
                shrink(x,fresh); if (fresh == 0) return false;
                for (int z : neighbors[x]) if (z != y) queue.emplace_back(z,x);
            }
        }
        return true;
    };
    auto forward = [&](int x, int a) {
        for (int y : neighbors[x]) {
            int allowed = table[x][y][a];
            if (assigned[y] != -1) { if (!(allowed & (1 << assigned[y]))) return false; }
            else { shrink(y,domains[y] & allowed); if (domains[y] == 0) return false; }
        }
        return true;
    };
    std::function<bool()> dfs = [&]() {
        ++nodes; int x = -1; std::tuple<int,int,int> best;
        for (int v = 0; v < n; ++v) if (assigned[v] == -1) {
            int degree = 0; for (int y : neighbors[v]) if (assigned[y] == -1) ++degree;
            auto key = std::make_tuple(popcount(domains[v]),-degree,v);
            if (x == -1 || key < best) { x = v; best = key; }
        }
        if (x == -1) { answer = assigned; return true; }
        std::vector<std::pair<int,int>> values;
        for (int a = 0; a < d; ++a) if (domains[x] & (1 << a)) {
            int score = 0;
            for (int y : neighbors[x]) if (assigned[y] == -1) score += popcount(domains[y] & ~table[x][y][a]);
            values.emplace_back(score,a);
        }
        std::sort(values.begin(),values.end());
        for (auto p : values) {
            int a = p.second; ++attempts; std::size_t mark = trail.size(); assigned[x] = a; shrink(x,1 << a);
            bool ok = forward(x,a); if (ok && use_ac3) ok = ac3();
            bool found = ok && dfs();
            undo(mark); assigned[x] = -1; if (found) return true;
        }
        return false;
    };
    bool ok = std::all_of(domains.begin(),domains.end(),[](int mask) { return mask != 0; });
    if (ok && use_ac3) ok = ac3();
    bool found = ok && dfs(); undo(0);
    assert(domains == initial && trail.empty()); for (int value : assigned) assert(value == -1);
    return {found,answer,nodes,attempts,pruned,revisions};
}
Edge neq(int x, int y, int d) {
    Edge edge{x,y,{} };
    for (int a = 0; a < d; ++a) for (int b = 0; b < d; ++b) if (a != b) edge.pairs.emplace_back(a,b);
    return edge;
}
struct Fixture { std::string name; std::vector<int> domains; std::vector<Edge> edges; int d; bool mac; };
std::vector<Fixture> fixtures() {
    std::vector<Edge> triangle = {neq(0,1,3),neq(0,2,3),neq(1,2,3)};
    std::vector<Edge> triangle2 = {neq(0,1,2),neq(0,2,2),neq(1,2,2)}, restore, australia;
    for (int x : {1,2,3}) {
        Edge edge{0,x,{} };
        for (int a = 0; a < 2; ++a) for (int b = 0; b < 3; ++b) if (a == 1 || b < 2) edge.pairs.emplace_back(a,b);
        restore.push_back(edge);
    }
    for (int x : {4,5}) {
        Edge edge{0,x,{} };
        for (int a = 0; a < 2; ++a) for (int b = 0; b < 3; ++b) if (a == 0 || b == 0) edge.pairs.emplace_back(a,b);
        restore.push_back(edge);
    }
    restore.push_back(neq(1,2,3)); restore.push_back(neq(1,3,3)); restore.push_back(neq(2,3,3));
    for (auto p : std::vector<std::pair<int,int>>{ {0,1},{0,2},{1,2},{1,3},{2,3},{2,4},{2,5},{3,4},{4,5} }) australia.push_back(neq(p.first,p.second,3));
    return { {"asymmetric / FC",{3,7,5},triangle,3,false},
             {"asymmetric / MAC",{3,7,5},triangle,3,true},
             {"triangle2 / FC",{3,3,3},triangle2,2,false},
             {"triangle2 / MAC",{3,3,3},triangle2,2,true},
             {"restore / MAC",{3,7,7,7,7,7},restore,3,true},
             {"Australia / MAC",std::vector<int>(7,7),australia,3,true} };
}
std::string fmt(const std::vector<int>& values) {
    std::string out = "[";
    for (std::size_t i = 0; i < values.size(); ++i) { if (i) out += ","; out += std::to_string(values[i]); }
    return out + "]";
}
int main() {
    bool ok = true;
    for (const auto& f : fixtures()) {
        auto original = f.domains; auto r = solve(f.domains,f.edges,f.d,f.mac);
        std::cout << f.name << ": " << (r.found ? "SAT" : "UNSAT") << ", nodes=" << r.nodes << ", attempts=" << r.attempts << ", pruned=" << r.pruned << ", revisions=" << r.revisions << '\n';
        std::cout << "assignment = " << fmt(r.answer) << '\n';
        bool valid = r.found == (f.name.find("triangle2") != 0);
        if (r.found) {
            valid = valid && r.answer.size() == f.domains.size();
            for (std::size_t x = 0; x < r.answer.size(); ++x) valid = valid && (f.domains[x] & (1 << r.answer[x]));
            for (const auto& edge : f.edges) valid = valid && std::find(edge.pairs.begin(),edge.pairs.end(),std::make_pair(r.answer[edge.x],r.answer[edge.y])) != edge.pairs.end();
        }
        ok = ok && valid && f.domains == original;
    }
    std::cout << "checks = " << (ok ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.ArrayList;
import java.util.ArrayDeque;
import java.util.Arrays;
import java.util.Collections;

public class Main {
    static class Edge {
        final int x,y; final ArrayList<int[]> pairs = new ArrayList<>();
        Edge(int x,int y) { this.x = x; this.y = y; }
        boolean has(int a,int b) { for (int[] p : pairs) if (p[0] == a && p[1] == b) return true; return false; }
    }
    static class Result {
        final boolean found; final int[] answer; final long nodes,attempts,pruned,revisions;
        Result(boolean found,int[] answer,long nodes,long attempts,long pruned,long revisions) {
            this.found = found; this.answer = answer; this.nodes = nodes; this.attempts = attempts; this.pruned = pruned; this.revisions = revisions;
        }
    }
    static Result solve(int[] initial,ArrayList<Edge> edges,int d,boolean useAc3) {
        int n = initial.length;
        if (n > 12 || d < 1 || d > 8) throw new IllegalArgumentException("invalid variable or value count");
        int full = (1 << d) - 1;
        for (int mask : initial) if (mask < 0 || (mask & ~full) != 0) throw new IllegalArgumentException("invalid domain mask");
        ArrayList<ArrayList<Integer>> neighbors = new ArrayList<>();
        for (int x = 0; x < n; ++x) neighbors.add(new ArrayList<>());
        int[][][] table = new int[n][n][d]; boolean[][] seen = new boolean[n][n];
        for (Edge edge : edges) {
            int x = edge.x,y = edge.y;
            if (x < 0 || x >= n || y < 0 || y >= n || x == y || seen[x][y]) throw new IllegalArgumentException("invalid or repeated edge");
            seen[x][y] = seen[y][x] = true; neighbors.get(x).add(y); neighbors.get(y).add(x);
            for (int[] p : edge.pairs) {
                int a = p[0],b = p[1];
                if (a < 0 || a >= d || b < 0 || b >= d) throw new IllegalArgumentException("invalid relation value");
                table[x][y][a] |= 1 << b; table[y][x][b] |= 1 << a;
            }
        }
        for (ArrayList<Integer> row : neighbors) Collections.sort(row);
        class Search {
            final int[] domains = initial.clone(),assigned = new int[n];
            final ArrayList<int[]> trail = new ArrayList<>();
            int[] answer = new int[0]; long nodes,attempts,pruned,revisions;
            Search() { Arrays.fill(assigned,-1); }
            void shrink(int x,int fresh) {
                int old = domains[x];
                if ((fresh & ~old) != 0) throw new AssertionError("domain expansion");
                if (fresh != old) { trail.add(new int[] {x,old}); domains[x] = fresh; pruned += Integer.bitCount(old) - Integer.bitCount(fresh); }
            }
            void undo(int mark) {
                while (trail.size() > mark) { int[] p = trail.remove(trail.size()-1); domains[p[0]] = p[1]; }
            }
            boolean ac3() {
                ArrayDeque<int[]> queue = new ArrayDeque<>();
                for (int x = 0; x < n; ++x) for (int y : neighbors.get(x)) queue.addLast(new int[] {x,y});
                while (!queue.isEmpty()) {
                    int[] p = queue.removeFirst(); int x = p[0],y = p[1]; ++revisions;
                    int old = domains[x],fresh = 0;
                    for (int a = 0; a < d; ++a) if ((old & (1 << a)) != 0 && (table[x][y][a] & domains[y]) != 0) fresh |= 1 << a;
                    if (fresh != old) {
                        shrink(x,fresh); if (fresh == 0) return false;
                        for (int z : neighbors.get(x)) if (z != y) queue.addLast(new int[] {z,x});
                    }
                }
                return true;
            }
            boolean forward(int x,int a) {
                for (int y : neighbors.get(x)) {
                    int allowed = table[x][y][a];
                    if (assigned[y] != -1) { if ((allowed & (1 << assigned[y])) == 0) return false; }
                    else { shrink(y,domains[y] & allowed); if (domains[y] == 0) return false; }
                }
                return true;
            }
            boolean dfs() {
                ++nodes; int x = -1,bestSize = d + 1,bestDegree = -1;
                for (int v = 0; v < n; ++v) if (assigned[v] == -1) {
                    int degree = 0; for (int y : neighbors.get(v)) if (assigned[y] == -1) ++degree;
                    int size = Integer.bitCount(domains[v]);
                    if (size < bestSize || (size == bestSize && degree > bestDegree)) { x = v; bestSize = size; bestDegree = degree; }
                }
                if (x == -1) { answer = assigned.clone(); return true; }
                ArrayList<int[]> values = new ArrayList<>();
                for (int a = 0; a < d; ++a) if ((domains[x] & (1 << a)) != 0) {
                    int score = 0;
                    for (int y : neighbors.get(x)) if (assigned[y] == -1) score += Integer.bitCount(domains[y] & ~table[x][y][a]);
                    values.add(new int[] {score,a});
                }
                values.sort((p,q) -> p[0] != q[0] ? Integer.compare(p[0],q[0]) : Integer.compare(p[1],q[1]));
                for (int[] p : values) {
                    int a = p[1]; ++attempts; int mark = trail.size(); assigned[x] = a; shrink(x,1 << a);
                    boolean ok = forward(x,a); if (ok && useAc3) ok = ac3();
                    boolean found = ok && dfs();
                    undo(mark); assigned[x] = -1; if (found) return true;
                }
                return false;
            }
        }
        Search search = new Search(); boolean ok = true;
        for (int mask : initial) if (mask == 0) ok = false;
        if (ok && useAc3) ok = search.ac3();
        boolean found = ok && search.dfs(); search.undo(0);
        if (!Arrays.equals(search.domains,initial) || !search.trail.isEmpty()) throw new AssertionError("undo failed");
        for (int value : search.assigned) if (value != -1) throw new AssertionError("assignment undo failed");
        return new Result(found,search.answer,search.nodes,search.attempts,search.pruned,search.revisions);
    }
    static Edge neq(int x,int y,int d) {
        Edge edge = new Edge(x,y);
        for (int a = 0; a < d; ++a) for (int b = 0; b < d; ++b) if (a != b) edge.pairs.add(new int[] {a,b});
        return edge;
    }
    static class Fixture {
        final String name; final int[] domains; final ArrayList<Edge> edges; final int d; final boolean mac;
        Fixture(String name,int[] domains,ArrayList<Edge> edges,int d,boolean mac) { this.name = name; this.domains = domains; this.edges = edges; this.d = d; this.mac = mac; }
    }
    static ArrayList<Fixture> fixtures() {
        ArrayList<Edge> triangle = new ArrayList<>(Arrays.asList(neq(0,1,3),neq(0,2,3),neq(1,2,3)));
        ArrayList<Edge> triangle2 = new ArrayList<>(Arrays.asList(neq(0,1,2),neq(0,2,2),neq(1,2,2)));
        ArrayList<Edge> restore = new ArrayList<>(),australia = new ArrayList<>();
        for (int x : new int[] {1,2,3}) {
            Edge edge = new Edge(0,x);
            for (int a = 0; a < 2; ++a) for (int b = 0; b < 3; ++b) if (a == 1 || b < 2) edge.pairs.add(new int[] {a,b});
            restore.add(edge);
        }
        for (int x : new int[] {4,5}) {
            Edge edge = new Edge(0,x);
            for (int a = 0; a < 2; ++a) for (int b = 0; b < 3; ++b) if (a == 0 || b == 0) edge.pairs.add(new int[] {a,b});
            restore.add(edge);
        }
        restore.add(neq(1,2,3)); restore.add(neq(1,3,3)); restore.add(neq(2,3,3));
        for (int[] p : new int[][] { {0,1},{0,2},{1,2},{1,3},{2,3},{2,4},{2,5},{3,4},{4,5} }) australia.add(neq(p[0],p[1],3));
        ArrayList<Fixture> out = new ArrayList<>();
        out.add(new Fixture("asymmetric / FC",new int[] {3,7,5},triangle,3,false));
        out.add(new Fixture("asymmetric / MAC",new int[] {3,7,5},triangle,3,true));
        out.add(new Fixture("triangle2 / FC",new int[] {3,3,3},triangle2,2,false));
        out.add(new Fixture("triangle2 / MAC",new int[] {3,3,3},triangle2,2,true));
        out.add(new Fixture("restore / MAC",new int[] {3,7,7,7,7,7},restore,3,true));
        out.add(new Fixture("Australia / MAC",new int[] {7,7,7,7,7,7,7},australia,3,true));
        return out;
    }
    static String fmt(int[] values) {
        StringBuilder b = new StringBuilder("[");
        for (int i = 0; i < values.length; ++i) { if (i > 0) b.append(','); b.append(values[i]); }
        return b.append(']').toString();
    }
    public static void main(String[] args) {
        boolean ok = true;
        for (Fixture f : fixtures()) {
            int[] original = f.domains.clone(); Result r = solve(f.domains,f.edges,f.d,f.mac);
            System.out.println(f.name + ": " + (r.found ? "SAT" : "UNSAT") + ", nodes=" + r.nodes + ", attempts=" + r.attempts + ", pruned=" + r.pruned + ", revisions=" + r.revisions);
            System.out.println("assignment = " + fmt(r.answer));
            boolean valid = r.found == !f.name.startsWith("triangle2");
            if (r.found) {
                valid = valid && r.answer.length == f.domains.length;
                for (int x = 0; x < r.answer.length; ++x) valid = valid && (f.domains[x] & (1 << r.answer[x])) != 0;
                for (Edge edge : f.edges) valid = valid && edge.has(r.answer[edge.x],r.answer[edge.y]);
            }
            ok = ok && valid && Arrays.equals(f.domains,original);
        }
        System.out.println("checks = " + (ok ? "PASS" : "FAIL"));
    }
}

Use Python 3.10 or later for int.bit_count(), compile C++ with g++ -std=c++14 -O2 main.cpp -o csp, or compile Java 17 or later with javac Main.java and run java Main. All three programs print these thirteen lines:

• asymmetric / FC: SAT, nodes=4, attempts=3, pruned=4, revisions=0
• assignment = [1,0,2]
• asymmetric / MAC: SAT, nodes=4, attempts=3, pruned=4, revisions=24
• assignment = [1,0,2]
• triangle2 / FC: UNSAT, nodes=3, attempts=4, pruned=8, revisions=0
• assignment = []
• triangle2 / MAC: UNSAT, nodes=1, attempts=2, pruned=8, revisions=14
• assignment = []
• restore / MAC: SAT, nodes=8, attempts=9, pruned=23, revisions=150
• assignment = [1,0,1,2,0,0]
• Australia / MAC: SAT, nodes=8, attempts=7, pruned=14, revisions=147
• assignment = [2,1,0,2,1,2,0]
• checks = PASS

Complexity, Correctness, and Invariants

Let $n$ be the number of variables, $d$ the value-universe size, $e$ the number of binary edges, and $\Delta$ the maximum graph degree. The supported domains and support sets fit in one word. The following operation bounds treat their bit operations and bounded-word population counts as constant cost.

MRV with dynamic degree costs $O(n+e)$ per DFS node. Scoring all values for a selected variable costs $O(d\Delta)$, followed by $O(d\log(d+1))$ sorting. Forward checking one attempted value costs $O(1+\Delta)$ because each neighbor filter is one mask intersection plus reversible bookkeeping.

A pair-by-pair AC-3 implementation commonly has an $O(ed^3)$ bound. Here support masks make one arc revision cost $O(d)$: loop over the tail's candidate values and test each support mask against the head domain. During one AC-3 invocation, each domain can shrink strictly at most $d$ times. A shrink queues incoming arcs, so there are $O(ed)$ arc visits overall, plus initialization. Thus this implementation has an $O(n+ed^2)$ bound per AC-3 call.

Search still has a potentially exponential tree. At depth $k$, at most $d^k$ assignments can have been attempted along distinct branches, regardless of the variable-ordering policy. The full node bound is $\sum_{k=0}^n d^k$, which is $O(d^n)$ for $d\ge2$ and $O(n+1)$ when $d=1$. Propagation and heuristics can dramatically reduce actual work without making general binary CSP solving polynomial.

Let $T$ be the number of DFS entries and $K$ the number of tried candidates. An informative MAC work bound, excluding model construction, is

$$O\!\left(n+ed^2+T\bigl(n+e+d\Delta+d\log(d+1)\bigr)+K\bigl(\Delta+n+ed^2\bigr)\right).$$

The initial term covers root checks and propagation, and $K\le dT$ whenever search occurs. In FC mode, omit all AC-3 work and retain candidate filtering and undo costs. Model construction initializes the dense support table in $O(n^2d)$ time, reads $M$ supplied allowed-pair entries, and sorts neighbor lists in $O(e\log(n+1))$ time.

The actual representation uses $O(n^2d)$ words for its dense support masks and $O(n+e)$ adjacency storage. Along one active search path, at most $nd$ domain values can be removed, so its trail has $O(nd)$ records; recursive value-order lists also use $O(nd)$ space. The non-deduplicated AC-3 queue has an $O(ed)$ bound. Excluding the caller's relation lists, logically live auxiliary storage is therefore $O(n^2d+nd+ed+n+e)$. A sparse relation-table representation can reduce the dense-table term.

Best-case behavior depends on the input and objective. An initial empty domain or a root propagation contradiction needs no DFS entries, although validation and model construction still cost work. An easy satisfiable instance still uses $n+1$ DFS entries here because singleton domains are explicitly assigned. There is no general average-case bound without an input distribution; fewer nodes alone do not establish a faster solver.

Soundness at a leaf follows because each committed value belongs to its domain and every binary edge was checked when its later endpoint was assigned. Completeness follows because no sound deletion removes an extendable solution and every surviving value is considered unless a full answer has already been copied. Search depth is at most $n$; propagation terminates because domains shrink monotonically within an invocation and new queue work is introduced only after a shrink.

Undo must restore the parent domains, not blindly restore every variable's original domain. A parent may already have made legitimate reductions. Saving old masks handles this distinction automatically. Root propagation is also trailed, letting the implementation restore its own entire working state before returning either verdict.

Variants, Pitfalls, and Optimizations

The most damaging bug is incomplete rollback. In the restore fixture, the failed $A=0$ branch removes value $2$ from all of $B,C,D$. Restoring only $A$ leaves those deletions behind. The later $A=1$ branch would then wrongly face a two-color triangle, losing the valid completion that uses a third value. The correct checkpoint restores all three neighbor domains as well as $A$.

Other design choices need separate reasoning:

• Incremental propagation: once the parent is arc consistent, initialize only arcs affected by the decision and the preceding FC deletions. Seed every relevant changed domain; propagating only from the assigned variable can miss reductions made in its neighbors.
• Queue deduplication and support caching: avoid repeated pending arcs or remember useful support locations. Any cached support must be checked against current domains, and branch-dependent metadata must remain valid after rollback.
• Stronger constraints: a global all-different propagator can detect Hall-set effects that a decomposition into pairwise inequalities misses. Arc consistency of the binary decomposition is a different, generally weaker filtering condition.
• Constraint decomposition: higher-arity relations can be encoded with auxiliary variables or decomposed where justified. A poor encoding can enlarge domains, weaken propagation, or change the problem if it omits relationships.
• Disconnected components: solve independent graph components separately and combine their assignments. The example's Tasmania variable is independent of the mainland coloring; the basic solver keeps it in the common search loop for transparency.
• Learning and backjumping: record explanations for conflicts and avoid rediscovering equivalent failures. Learned constraints must be logically implied, and backjump targets must match the recorded dependencies.
• Alternative search: min-conflicts and other local methods can be useful for finding solutions, but a failure to find one within a limit is not a proof of UNSAT. This backtracking solver reports failure only after complete exclusion.
• Counting and optimization: enumerating solutions, proving uniqueness, and minimizing a cost require different stopping conditions. The first satisfying assignment does not answer those questions.

Practice Problems and Real-World Use

Compare small instances against exhaustive assignment enumeration, using the original domains and allowed-pair lists as the independent specification. Test asymmetric relations, reversed edge orientations, empty domains, empty relations, disconnected graphs, singleton domains, and the empty CSP. FC and MAC must agree on satisfiability, even when their first returned assignment or work differs.

Check each returned assignment directly, preserve every input, and audit current domains after each deletion and rollback. For a stronger oracle on small instances, enumerate every completion extending the current assignment and verify that propagation has not removed any value used by those completions. Use the restore fixture to ensure a failed early branch cannot poison the later successful branch.

For further study, see the AIMA CSP implementation and examples for backtracking, filtering, and map-coloring formulations. The solver here uses deterministic ordering, bit-mask support tables, explicit domain checkpoints, and the precise counter definitions above so its behavior can be checked across languages.