Solving the Magic Squares with constraint logic programming. Examples in prolog and python.

magic_squares.pl

% Solving Magic Squares using Sicstus PROLOG

:- use_module(library(clpfd)).
:- use_module(library(lists)).

matrixDomain([], _).
matrixDomain([Row | Matrix], N) :-
  NSquared is N * N,
  length(Row, N),
  domain(Row, 1, NSquared),
  matrixDomain(Matrix, N).

equal_sum_along_rows([], _).
equal_sum_along_rows([Row | Matrix], Sum) :-
  sum(Row, #=, Sum),
  equal_sum_along_rows(Matrix, Sum).

constraint_diagonals(Matrix, Sum) :-
  get_downwards_diagonal(Matrix, Diag1, 1),
  sum(Diag1, #=, Sum),
  get_upwards_diagonal(Matrix, Diag2, 1),
  sum(Diag2, #=, Sum).

get_downwards_diagonal([], [], _).
get_downwards_diagonal([Row | Matrix], [El | Diag], Idx) :-
  nth1(Idx, Row, El),
  NIdx is Idx + 1,
  get_downwards_diagonal(Matrix, Diag, NIdx).

get_upwards_diagonal([], [], _).
get_upwards_diagonal(Matrix, [El | Diag], Idx) :-
  last(Matrix, Row),
  select(Row, Matrix, TopRows),
  nth1(Idx, Row, El),
  NIdx is Idx + 1,
  get_upwards_diagonal(TopRows, Diag, NIdx).


magic_squares(Matrix, N) :-
  length(Matrix, N),
  matrixDomain(Matrix, N),

  append(Matrix, FlatList),
  all_distinct(FlatList),

  Sum is N * (N*N + 1) // 2,
  equal_sum_along_rows(Matrix, Sum),

  transpose(Matrix, TransposedMatrix),
  equal_sum_along_rows(TransposedMatrix, Sum),

  constraint_diagonals(Matrix, Sum),

  labeling([ffc], FlatList).
%% magic_squares(M, 4), printMatrix(M).


%% PRINT MATRIX
printMatrix([]).
printMatrix([Row | Matrix]) :-
  printRow(Row), nl,
  printMatrix(Matrix).

printRow([]).
printRow([El | Row]) :-
  write(El), write(' '),
  printRow(Row).

magic_squares.py

# Solving Magic Squares using a python constraint solver
# https://github.com/python-constraint/python-constraint

import sys
from constraint import *

def magic_square(n):
  problem = Problem()
  problem.addVariables(range(0, n**2), range(1, n**2 + 1))
  problem.addConstraint(AllDifferentConstraint(), range(0, n**2))

  sumResult = n * (n**2 + 1) // 2
  for row in range(n):
    problem.addConstraint(ExactSumConstraint(sumResult), [row * n+i for i in range(n)])

  for col in range(n):
    problem.addConstraint(ExactSumConstraint(sumResult), [col + n*i for i in range(n)])

  # restrict diagonal top-left to bot-right
  problem.addConstraint(ExactSumConstraint(sumResult), [i for i in range(0, n**2, n+1)])

  # restrict diagonal bot-left to top-right
  problem.addConstraint(ExactSumConstraint(sumResult), [i for i in range(n**2-n, 0, -(n-1))])

  return problem.getSolutionIter()

def print_square(square, n):
  for i in range(len(square)):
    if (i % n) == 0:
      print()
    print(square[i], " ", end="")


if (len(sys.argv) <= 1):
  print("Usage: " + sys.argv[0] + " N")
  sys.exit()

n = int(sys.argv[1])
solutions = magic_square(n)

for sol in solutions:
  print_square(sol, n)
  print()

magic_squares_ortools.py

# Solving Magic Squares using Google's ortools
# https://github.com/google/or-tools

from ortools.constraint_solver import pywrapcp

def main(n=4):
    # Create the solver.
    solver = pywrapcp.Solver('Magic square')
    # declare variables
    x = {}
    for i in range(n):
        for j in range(n):
            x[(i, j)] = solver.IntVar(1, n*n, 'x(%i,%i)' % (i, j))
    x_flat = [x[(i,j)] for i in range(n) for j in range(n)]

    # the sum
    s = solver.IntVar(1, n*n*n,'s')
    
    # constraints  
    solver.Add(solver.AllDifferent(x_flat, True))
    [solver.Add(solver.Sum([x[(i,j)] for j in range(n)]) == s) for i in range(n)]    
    [solver.Add(solver.Sum([x[(i,j)] for i in range(n)]) == s) for j in range(n)]        
    solver.Add(solver.Sum([ x[(i,i)]     for i in range(n)]) == s) # diag 1
    solver.Add(solver.Sum([ x[(i,n-i-1)] for i in range(n)]) == s) # diag 2       

    # solution and search
    solution = solver.Assignment()
    solution.Add(x_flat)
    solution.Add(s)
    db = solver.Phase(x_flat,
                      solver.CHOOSE_FIRST_UNBOUND,                     
                      solver.ASSIGN_CENTER_VALUE
                      )    
    solver.NewSearch(db)

    # output
    num_solutions = 0
    while solver.NextSolution():
        print("s:", s.Value())
        for i in range(n):
            for j in range(n):
                print("%2i" % x[(i,j)].Value(), end="")
            print()
        print()
        num_solutions += 1        
    solver.EndSearch()

    print()
    print("num_solutions:", num_solutions)
    print("failures:", solver.Failures())
    print("branches:", solver.Branches())
    print("wall_time:", solver.WallTime())

if __name__ == '__main__':
  main()