Solving Sudoku with Constraint Programming

Over the past few months, I’ve shared a couple of Sudoku solvers here — some built with Dynamic Programming, others with Mixed-Integer Linear Programming (MILP). Both approaches work and are great learning exercises, but they each have their quirks when it comes to constraint formulation:

  • In MILP, expressing “all-different” constraints means either crafting a set of linear inequalities or introducing auxiliary binary variables with Big-M tricks.
  • Dynamic Programming is more algorithmic and procedural, which makes it less declarative.

This time, I wanted something cleaner — so I turned to Constraint Programming (CP) with Google OR-Tools.

The beauty is: With CP, Sudoku rules can be expressed almost exactly as you’d explain them to a human:

  • Numbers in each row: AllDifferent
  • Numbers in each column: AllDifferent
  • Numbers in each 3×3 block: AllDifferent

No manual encoding of constraints, no binary variables for every number — just the logical structure of the puzzle.

To push it further, I ran the model on the infamous “world’s hardest Sudoku”, and it solved it in fractions of a second. Best part? I didn’t need to tweak the model at all.

Why CP shines here:

  • Declarative modeling: You state the rules, the solver figures out the search.
  • Concise code: A handful of loops and AllDifferent constraints define the whole puzzle.
  • Solver efficiency: CP explores Sudoku search spaces very effectively.
  • Extensibility: Want Killer Sudoku, Diagonal Sudoku, or Samurai Sudoku? Just add more constraints.

I’ve attached the full code below — the core is really just variable declarations plus three loops of AllDifferent constraints.

I’m curious:

  • Who here has experimented with OR-Tools’ CP solver, and how do you think it stacks up against MILP for puzzles like this?
  • Has anyone tried combining CP with custom search heuristics for exotic Sudoku variants?
"""Sudoku solver using Google OR-Tools CP-SAT.
This module exposes a single `solve` function that accepts a 9x9 grid
with zeros for blanks and returns a completed Sudoku grid.
"""

from ortools.sat.python import cp_model

def solve(puzzle):
    """Solve a standard 9×9 Sudoku with OR-Tools CP-SAT.

    The puzzle is modeled as 81 integer decision variables with the
    usual Sudoku constraints:
    - Rows contain numbers 1..9 exactly once.
    - Columns contain numbers 1..9 exactly once.
    - Each 3×3 block contains numbers 1..9 exactly once.
    Provided clues (non-zero entries) are enforced as fixed equalities.

    Args:
        puzzle: A 9×9 list of lists of ints. Use 0 to denote empty cells;
            filled cells must be in the range 1..9.

    Returns:
        A 9×9 list of lists of ints representing a solved Sudoku grid.

    """
    model = cp_model.CpModel()
    # Decision variables: one IntVar per cell with domain 1..9
    # cells[(r, c)] corresponds to the value in row r, column c (0-indexed).
    cells = {}
    for r in range(9):
        for c in range(9):
            cells[(r, c)] = model.NewIntVar(1, 9, f'cell_{r}_{c}')
    
    # Constraints
    # Rows & columns must contain all digits 1..9 exactly once.
    for i in range(9):
        model.AddAllDifferent([cells[(i, j)] for j in range(9)])
        model.AddAllDifferent([cells[(j, i)] for j in range(9)])
    # 3×3 subgrid constraints:
    for br in range(3):
        for bc in range(3):
            block = [cells[(r, c)]
                     for r in range(br*3, (br+1)*3)
                     for c in range(bc*3, (bc+1)*3)]
            model.AddAllDifferent(block)
    
    # Fix given clues (non-zero cells) to their provided values.
    for r in range(9):
        for c in range(9):
            if puzzle[r][c] != 0:
                model.Add(cells[(r, c)] == puzzle[r][c])
    
    solver = cp_model.CpSolver()
    solver.Solve(model)
    return [[solver.Value(cells[(r, c)]) for c in range(9)] for r in range(9)]

def main():        
    # Hardest ever Sudoku
    initial_grid = [[8, 0, 0, 0, 0, 0, 0, 0, 0],
                         [0, 0, 3, 6, 0, 0, 0, 0, 0],
                         [0, 7, 0, 0, 9, 0, 2, 0, 0],
                         [0, 5, 0, 0, 0, 7, 0, 0, 0],
                         [0, 0, 0, 0, 4, 5, 7, 0, 0],
                         [0, 0, 0, 1, 0, 0, 0, 3, 0],
                         [0, 0, 1, 0, 0, 0, 0, 6, 8],
                         [0, 0, 8, 5, 0, 0, 0, 1, 0],
                         [0, 9, 0, 0, 0, 0, 4, 0, 0],
                     ]
    solution = solve(initial_grid)
    print(solution)

if __name__ == "__main__":
    main()