🧩 Constraint Solving POTD:Problem of the Day: Magic Squares

[github-actions[bot]]

Problem Statement

A magic square of order n is an n × n grid filled with distinct integers (typically 1 to n2) such that every row, column, and main diagonal sum to the same value — the magic constant M = n(n2 + 1) / 2.

For a 4 × 4 magic square, the magic constant is M = 34. Here's a valid solution:

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

Each row, column, and both main diagonals sum to 34.

Input/Output Specification:

• Input: Order n (and optionally, constraints: some cells pre-filled)
• Output: An n × n grid where:
  • Each cell contains a distinct value from {1, 2, ..., n2}
  • All row sums = M
  • All column sums = M
  • Both main diagonal sums = M

---

Why It Matters

Historical Fascination & Recreational Mathematics: Magic squares appear in ancient Chinese texts and Islamic geometry, and remain a classic subject in recreational mathematics and puzzle design.

Constraint Propagation Benchmark: Magic squares are excellent teaching examples for constraint propagation, global constraints (e.g., all_different), and symmetry-breaking in constraint solvers — they're often used in solver tutorials and research papers.

Combinatorial Search: Computing all magic squares for a given order is a hard combinatorial problem; for n=4, there are 880 essentially different solutions, making exhaustive enumeration computationally interesting.

---

Modeling Approaches

Approach 1: Constraint Programming (All-Different + Linear Sums)

The most natural CP formulation uses array variables and global constraints.

Decision Variables:

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

Constraints:

1. all_different(x[0][0], x[0][1], ..., x[n-1][n-1]) — all cells are distinct
2. For each row i: sum(x[i][0], x[i][1], ..., x[i][n-1]) = M
3. For each column j: sum(x[0][j], x[1][j], ..., x[n-1][j]) = M
4. Main diagonal: sum(x[0][0], x[1][1], ..., x[n-1][n-1]) = M
5. Anti-diagonal: sum(x[0][n-1], x[1][n-2], ..., x[n-1][0]) = M

Trade-offs:

• Pros: Highly expressive, leverages strong domain propagation for all_different and linear sum constraints.
• Cons: The search space is large (n2! possibilities before constraint propagation), requiring effective variable ordering.

Example Model (MiniZinc)

int: n = 4;
int: m = n * (n * n + 1) div 2;  % magic constant

array[1..n, 1..n] of var 1..n*n: x;

% All cells are distinct
constraint all_different([x[i,j] | i,j in 1..n]);

% Row sums
constraint forall(i in 1..n) (sum(j in 1..n) (x[i,j]) = m);

% Column sums
constraint forall(j in 1..n) (sum(i in 1..n) (x[i,j]) = m);

% Diagonal sums
constraint sum(i in 1..n) (x[i,i]) = m;
constraint sum(i in 1..n) (x[i,n+1-i]) = m;

solve satisfy;

Approach 2: Integer Linear Programming (ILP) with 0–1 Assignment Variables

An alternative formulation treats this as a combinatorial assignment problem.

Decision Variables:

• y[i][j][v] ∈ {0, 1}: binary variable indicating whether cell (i, j) contains value v

Constraints:

1. Each cell contains exactly one value: ∀(i, j): sum_v y[i][j][v] = 1
2. Each value appears exactly once: ∀v: sum_{i,j} y[i][j][v] = 1
3. Row/column/diagonal sum constraints expressed via: sum_{i,j,v} v · y[i][j][v] = M for each constraint

Trade-offs:

• Pros: Amenable to ILP solvers (e.g., CPLEX, Gurobi); LP relaxations provide strong bounds.
• Cons: Requires n4 binary variables; less natural for human intuition; LP relaxations often loose.

---

Key Techniques

1. Global Constraint Propagation (All-Different + Sum)

The all_different constraint, combined with domain reduction from sum constraints, eliminates large portions of the search space without explicit branching.

• Arc consistency on sum constraints removes infeasible value combinations.
• Strong propagators (like table, element) further prune domains.

2. Symmetry Breaking

Magic squares have rotational and reflectional symmetries (except for distinct solutions). Symmetry-breaking constraints reduce search by a factor of up to 8:

• Fix cell x[0][0] = 1 (or the smallest value)
• Enforce x[0][1] < x[n-1][0] (diagonal ordering constraint)
  These eliminate redundant branches without removing valid solutions.

3. Variable & Value Ordering Heuristics

• Variable ordering: First assign "most-constrained" cells (those with smallest domain).
• Value ordering: Try values that leave the most flexibility for neighboring constraints (e.g., middle values before extremes).

---

Challenge Corner

Open Questions:

1. Symmetry Reduction: How can you express all non-trivial symmetries of a magic square as constraints? Can you enumerate symmetries and add constraints that keep only a canonical representative?

2. Incomplete/Partial Magic Squares: Suppose you have a partially filled magic square (some cells pre-assigned). How does this change the solving complexity? Can you design an algorithm that incrementally solves such puzzles?

3. Computational Complexity: What is the growth rate of the number of distinct n × n magic squares as n increases? How does this relate to the hardness of constraint solving?

4. Generalization: How would you extend this to magic cubes (3D grids)? What new challenges emerge?

---

References

1. Rossi, F., van Beek, P., & Walsh, T. (Eds.). Handbook of Constraint Programming. Elsevier, 2006. (Ch. 2–3 on foundations; Ch. 34 on applications)

2. Hooker, J. N. Logic-based Methods for Optimization: Combining Optimization and Constraint Satisfaction. Wiley, 2000. (Ch. 4: search and pruning)

3. Eppstein, D. "Magic Squares" – Mathematics of Magic Squares and Cubes. University of California course notes. (Accessible introduction to construction methods and symmetry analysis)

4. Cheng, K. C., & Yap, R. H. C. "Applying local search to the Sudoku problem." ICGA Journal, 2008. (Discusses constraint models for Latin square-like problems, which are closely related)

---

Next Problem Preview: Expect a scheduling problem next! 🗓️

Generated by 🧩 Constraint Solving — Problem of the Day · haiku45 · 5.91 AIC · ⊞ 5.8K · ◷

• expires on Aug 14, 2026, 3:41 AM UTC-08:00