🧩 Constraint Solving POTD:Magic Squares

[github-actions[bot]]

Problem of the Day: Magic Squares

Problem Statement

A magic square is an n×n grid filled with distinct integers from 1 to n2 such that every row, column, and main diagonal sums to the same value, called the magic constant.

For a 4×4 magic square, the magic constant is (1 + 2 + ... + 16) / 4 = 136 / 4 = 34.

Concrete Instance (4×4):

 ?  ?  ?  ?
 ?  ?  ?  ?
 ?  ?  ?  ?
 ?  ?  ?  ?

Each cell must contain a unique integer from 1 to 16. Each row, column, and the two main diagonals must sum to 34.

One valid solution:

 1  14  15   4
12   7   6   9
 8  11  10   5
13   2   3  16

Input/Output Specification:

• Input: Integer n (grid size)
• Output: An n×n grid of integers 1 to n2 where all rows, columns, and diagonals sum to n(n2 + 1)/2

Why It Matters

Recreation & Puzzle Design: Magic squares appear in recreational mathematics, puzzle books, and educational settings worldwide. Generating magic squares automatically is essential for creating puzzle benchmarks.

Combinatorial Explosion Test Case: Magic squares provide a challenging benchmark for constraint solvers because the search space is enormous (16! possible arrangements for 4×4) yet highly constrained—there are only 880 distinct 4×4 magic squares (symmetries aside).

Symmetry & Structure Learning: Solving magic squares efficiently teaches solvers about exploiting symmetry, using global constraints, and combining multiple propagation strategies.

Modeling Approaches

Approach 1: Pure Constraint Programming (CP)

Paradigm: Constraint Programming with global constraints
Decision variables:

• grid[i][j] ∈ {1..n2} for each cell (i, j)

Constraints:

AllDifferent(all grid variables)
∀i: sum(grid[i][0..n-1]) = magic_constant
∀j: sum(grid[0..n][j]) = magic_constant
sum(grid[0][0], grid[1][1], ..., grid[n-1][n-1]) = magic_constant
sum(grid[0][n-1], grid[1][n-2], ..., grid[n-1][0]) = magic_constant

Trade-offs:

• ✓ Expressive and natural modeling
• ✓ Strong propagation via AllDifferent and sum constraints
• ✓ Scales well with dedicated CP solvers (Choco, OR-Tools)
• ✗ Requires a dedicated CP solver

Approach 2: SAT/MIP Encoding

Paradigm: Boolean SAT or Mixed-Integer Programming
Decision variables:

• Binary variables x[i][j][v] = 1 if cell (i, j) contains value v, else 0
• Continuous variables row_sum[i] in an MIP formulation

Constraints:

∀(i,j): sum_v(x[i][j][v]) = 1           [each cell has exactly one value]
∀v: sum_{i,j}(x[i][j][v]) = 1          [each value used exactly once]
∀i: sum_{j,v}(j * x[i][j][v]) = magic_constant  [row sum constraint]
∀j: sum_{i,v}(i * x[i][j][v]) = magic_constant  [column sum constraint]
[similar for diagonals]

Trade-offs:

• ✓ Directly encodes as SAT/MIP with standard solvers
• ✓ Leverages mature SAT/MIP technology (CaDiCaL, Gurobi, CPLEX)
• ✗ Requires O(n4) variables (vs. O(n2) in CP)
• ✗ Weaker propagation without specialized global constraints

Key Techniques

1. Global Constraints & Strong Propagation

The AllDifferent constraint combined with sum constraints creates tight domain reduction. Arc consistency algorithms can eliminate many values early, dramatically pruning the search tree. Specialized algorithms like range-consistency propagation are particularly effective.

2. Symmetry Breaking

Magic squares exhibit rotational and reflectional symmetry. Fixing the position of certain values (e.g., placing 1 in the top-left corner, or enforcing an ordering on the top row) eliminates equivalent solutions without losing valid ones. This can reduce search time by orders of magnitude.

3. Smart Variable & Value Ordering

• Variable ordering: Prioritize cells in constrained positions (e.g., the center cell in odd-sized squares is involved in more constraints).
• Value ordering: Try values that sum closer to the magic constant first in high-degree cells.
• Restart strategies: Periodic restarts with different orderings can escape local dead-ends in large instances.

Challenge Corner

1. Symmetry Elimination: Can you formulate symmetry-breaking constraints that eliminate rotations and reflections without eliminating any fundamentally different magic squares? (Hint: Consider canonical orderings or structural properties.)

2. Scaling to Larger Grids: The 5×5 case has over 275 million magic squares. How would you efficiently sample or enumerate a representative subset?

3. Hybrid Approach: Could you combine constraint programming (for tight propagation) with a local search phase (to escape local optima in partial solutions) to speed up search?

References

• Rossi, F., van Beek, P., & Walsh, T. (Eds.). Handbook of Constraint Programming. Elsevier, 2006. [Chapters on CSP modeling and solving techniques]
• Bouvier, J., & Saux, Y. "Magic Squares: A Survey of Known Results and Trends." Journal of Recreational Mathematics, 2003. [Comprehensive overview of mathematical properties]
• Hendrix, H., et al. "Symmetry Breaking in Constraint Programming." ACM Computing Surveys, 2015. [Detailed techniques for symmetry elimination]
• Google OR-Tools Documentation: (developers.google.com/redacted) [Practical example and solver guide]

---

Posted on 2026-06-04

Generated by 🧩 Constraint Solving — Problem of the Day · haiku45 30.8K · ◷

• expires on Jun 11, 2026, 12:28 PM UTC

[github-actions[bot]]

This discussion has been marked as outdated by Constraint Solving — Problem of the Day.

A newer discussion is available at Discussion #37106.