Multidimensional Convex Polyhedra and Searching for Magic Squares

Christopher J. Henrich

(Buffalo, N.Y.)

Motivation: Magic Squares

Needs["Convex`ConvexPolyhedra`"];

A magic square of order n is a square array of the numbers {1 ,..., (n^2)} suh that the rows,,the columns, and the diagonals all add up to the same sum. Example:

Grid[{{16,3,2,13},{5,10,11,8},{9,6,7,12},{4,15,14,1}},Frame->True]

Almost half the squares of order 4 follow an algebraic pattern. (See http://www.mathinteract.com). What about the others?

Goal: a list of all squares of order 4; shorten the list by exploiting symmetries

The Search Problem

16 values.

CellNames = {a,b,c,d,e,f,g,h,i,k,l,m,n,o,p,q};

CellList is a list of the cells in the order in which they are to be given trial values in the search. Choosing a different order is an easy way of reorganizing the search procedure.

CellList = CellNames;

NameGrid = Grid[Partition[CellNames,4], Frame->True]

Equality constraints

a + b + c + d = 34; nine others. (There is one linear relation between them.)

In Convex`ConvexPolyhedra, an equality constraint is represented by a linear expression which is required to be zero...

MagicEqualities = { a + b + c + d - 34, e + f + g + h -34,
i + k + l + m - 34, n + o + p + q - 34,
a + e + i + n - 34, b + f + k + o - 34,
c + g + l + p - 34, d + h + m + q - 34,
a + f + l + q - 34, d + g + k + n - 34};

Inequality constraints

... and an inequality constraint is represented by a linear expression that is required to be ≤ 0.

Values from 1 to 16

1 ≤ a; fifteeen more like this;
a ≤ 16; fifteen more like this.

LowerBounds = ((1-#)&) /@ CellNames;

UpperBounds = ((#-16)&) /@ CellNames;

Inequality Constraints for symmetries of the problem

Because the values are meant to be distinct integers, I changed inequalities such as a < d into a ≤ d-1, or a-d+1≤0.

SymmetryInequalities = { a-d+1, a-f+1, a-g+1, a-k+1, a-l+1, a-n+1, a-q+1} ~ Join ~
{b-c+1, b-e+1, b-i+1};

Inequalities = Join[LowerBounds,UpperBounds, SymmetryInequalities];

Symmetries of the problem

If the square is rotated through 90°, it becomes another magic square.
The smallest of {a,d, n, q} may be rotated into the upper left corner. Therefore we may assume
a <d,         a < n,       a<q.
A reflection around the diagonal changes the square into another magic square. It is convenient to take advantage of this by assuming b < e.
Another symmetry allows us to assume
b < c,       b < i.
Another allows us to assume 
a < f,    a<g,    a<k,    a<l.

Because the values are meant to be distinct integers, I changed inequalities such as a < d into a ≤ d-1, or a-d+1≤0.

SymmetryInequalities = { a-d+1, a-f+1, a-g+1, a-k+1, a-l+1, a-n+1, a-q+1} ~ Join ~
{b-c+1, b-e+1, b-i+1};

Inequalities = Join[LowerBounds,UpperBounds, SymmetryInequalities];

The values must be distinct integers

The Search Space

The magic inequalities define a 7-dimensional affine subspace of (ℝ^16). The inequalities cut this down to a convex polyhedron, which I will calll SearchPolyhedron.

The Search Strategy

Representing the solutions

A magic square is represented by the expression MagicSquare[assignmentList], where assignmentList is a list of Rules to replace the symbols in CellList by the integer values which are to make up the magic square.

MakeMagicSquare[symbolList, valuesList] assumes that symbolList is a list of symbols and valuesList is an equally long list of values. It returns a MagicSquare containing a list of Rules, to replace each symbol by the respective value.

MakeMagicSquare[symbolList_, valueList_] := With[{f=(#1->#2)&},MagicSquare[Thread[f[symbolList, valueList]]]];

A solution point is a vector ((Subscript[v, 1]) ,..., (Subscript[v, 16])) of distinct integers which lies in SearchPolyhedron.

A partial solution point, of dimension α, is a vector ((Subscript[v, 1]) ,..., (Subscript[v, α])) of distinct integers that satisfy all the equaiities and inequalities imposed by requiring them to be the initial α coordinates of a point in SearchPolyhedron.  That is, a partial solution point lies in a projection of SearchPolyhedron.

MagicSquareList will be a list of all magic squares corresponding to solution points. It will be convenient to construct this as a nested list, and flatten it only as the last operation.

A family of partial problems

A partial solution point, of dimension α, is a vector ((Subscript[v, 1]) ,..., (Subscript[v, α])) of distinct integers that satisfy all the equaiities and inequalities imposed by requiring them to be the initial α coordinates of a point in SearchPolyhedron.  That is, a partial solution point lies in a projection of SearchPolyhedron.

If partial is a partial solution point, then MagicSquarePartialList[partial] is a nested list of magic squares constructed from solutions points whose initial coordinates agree with partial.

If partial is a partial solution point ((Subscript[v, 1]) ,..., (Subscript[v, α])) of dimension α≤16, then Extensions[partial] will be a list of all partial solution points of the from ((Subscript[v, 1]) ,..., (Subscript[v, α]), w).

MagicSquarePartialList[partial_] := 
If[Length[partial] == 16, 
MakeMagicSquare[CellList, partial], (*else*)
MagicSquarePartialList /@ Extensions[partial]];

NextCoordinates[partial] will be a list of all the integers w such that Append[partial, w] is a partialsolution.

Extensions[partial_] := Append[partial, #]& /@ NextCoordinates[partial];

The inequalities and equalities which characterize SearchPolyhedron imply that the members of NextCoordinates[partial] lie in the range between two values, which need not be integers and need not be distinct. NextCoordinateInterval[partial] will be a list of these two values.

PolyhedronIntersection is the function that duilds a description of a polyhedron from sets of constraints.

In Convex`ConvexPolyhedra, a "full constraint" is an equality, that restricts a point to lie on a hyperplane; a "half constraint" is an inequality, that restricts a point to lie in a half space.

IntegersInInterval[{lower_,upper_}] := Range[Ceiling[lower], Floor[upper]];

NextCoordinates[partial_] := IntegersInInterval[NextCoordinateInterval[partial]] ~ Complement ~ partial;

Before trying to define NextCoordinateInterval, let is look more closely at convex polyhedroa.

Convex Polyhedra

Defined as intersection of halfspaces...

Graphics[{{LightGreen,Polygon[{{0,2},{2,4},{6,4},{8,2},{4,0}}]},Line[{{-1,1},{3,5}}],Arrow[{{0.75,3.25},{1.25,2.75}}],
Line[{{1,4},{7,4}}],Arrow[{{4,4.35},{4,3.65}}],
Line[{{5,5},{9,1}}],Arrow[{{7.25,3.25},{6.75,2.75}}],
Line[{{3,-0.5},{9,2.5}}],Arrow[{{6.25,0.5},{5.75,1.5}}],
Line[{{-1,2.5},{5,-0.5}}],Arrow[{{1.75,0.5},{2.25,1.5}}]}
]

... or as convex hull of generators.

Graphics[{{EdgeForm[Thick],LightGreen,Polygon[{{0,2},{2,4},{6,4},{8,2},{4,0}}]},
Disk[{0,2},0.15],Disk[{2,4},0.15],Disk[{6,4},0.15],Disk[{8,2},0.15],Disk[{4,0},0.15]}]

Double description: maintain both generators and hyperplanes.

Operations include intersection with hyperplane...

... and  projection along a line.

With the double description, these operations are similar.

Using Convex`ConvexPolyhedra

PolyhedronIntersection is the function that builds a description of a polyhedron from sets of constraints.

In Convex`ConvexPolyhedra, a "full constraint" is an equality, that restricts a point to lie on a hyperplane; a "half constraint" is an inequality, that restricts a point to lie in a half space.

Convert from symbolic form to list of coefficients

Convex`Polyhedra represents a vector in (ℝ^n) by a list of n components. A constraint is a linear function with a constant term, represented by a list of n+1 components.

Row1Constraint = a+b+c+d-34;

{1,1,1,1,  0,0,0,0,  0,0,0,0,  0,0,0,0,  -34};

The conversion from the first form to the second needs to be automated. The strategy is to define a set of rules, each of which transforms one of the names of a cell into a list gererated by UnitVector.

ListVectorRules[variableList_] := With[{dd=Length[variableList]},MapIndexed[(#1->UnitVector[dd+1,First[#2]])&,variableList]];

CellRules = ListVectorRules[CellList]

How about extracting the constant term from a constraint? I chose to make a set of rules that transform each cell name to 0, and apply it to the constraint.

SetZeroRules[variableList_] := ((#1->0)&)/@variableList;

ConstantRules = SetZeroRules[CellNames]

Row1Constraint/.ConstantRules

MakeListVector[variableList_, linearExpr_] := With[{vr = ListVectorRules[variableList], constantTerm=(linearExpr /. SetZeroRules[variableList]), d=Length[variableList]},((linearExpr-constantTerm) /.vr) +constantTerm*UnitVector[d+1,d+1] ];

MakeListVector[CellList, Row1Constraint]

Constructing the Search Space

Constructing the search space

FullConstraints = MakeListVector[CellList, #]& /@ MagicEqualities;

HalfConstraints = MakeListVector[CellList, #]& /@ Inequalities;

{compTime, SearchPolyhedron} = Timing[ PolyhedronIntersection[FullConstraints, HalfConstraints]];compTime

To solve the family of problems for different dimensions, we need the projections of SearchPolyhedron on coordinate subspaces corresponding to subsets of cell names, such as {a,b,c,d,e}. Convex`ConvexPolyhedra makes such projections one axis at a time, with the function ProjectPolyhedronAlongAxis. We project successively along the axes running from 16 down to 2.

{compTime,PartialSearchPolyhedra} = Timing[Reverse[FoldList[ProjectPolyhedronAlongAxis,SearchPolyhedron,Reverse[Range[2,16]]]]]; compTime

An example

PartialSearchPolyhedra[[5]]

Reformatting the description for NextCoordinateInterval

If partial is of length α.NextCoordinateInterval[partial] is a description of the intersection of a partial search polyhedron with a line parallel to the α+1-st coordinate axis.

The information about the polyhedron, needed to compute this, will be called IntersectionData.

The polyhedron might or might not be confined to a hyperplane which intersects the d+1-st coordinate axis in a single point. 

If it is thus confined, then the intersection of the polyhedron with another line parallel to that axis is also a single point, whose d+1-st coordinate is a linear function of the first d coordinates. In this case, the intersection data are given by a PointIntersection.

If the polyhedron is not confined to a hyperplane, then the intersection is an interval. Its lower endpoint is the maximum of a collection of linear functions, and its upper endpoint is the minimum of another such collection. The intersection data are given by an IntervalInterection.

PointSection

An object  PointIntersection[coefficientList] contains a vector of dimension α+1, giving the coefficients and the constant term of a linear function.

IntervalSection

An object IntervalSection[{lowerBounds, upperBounds}] contains two lists of α+1-dimensional vectors. Each vector gives the coefficients and constant term of a linear function. The first list is used to determine the lower endpoint of an interval, and the second to determine the upper endpoint.

Constructing "IntersectionData"

MakeIntersectionData converts a Polyhedron description into the form of an IntervalIntersection or a PointIntersection.

FormVector is a convenience function.

FormVector[vector_] := With[{divisor=Part[vector,-2]},(-1/divisor)*Drop[vector,{-2}]];

MakeIntersectionData[search_] := Module[{d,hyp,nzx,hsp},
d = AmbientDimension[search];
hyp = Hyperplanes[search];
hsp = Halfspaces[search];
nzx = Select[Range@Length@hyp,(hyp[[#,d]] ≠ 0)&];
If[Length[nzx]==0,
(*No hyperplane constraint affects the d-th coordinate. *)
IntervalIntersection[FormVector /@ Select[hsp,(Part[#,d]<0)&], FormVector /@ Select[hsp,(Part[#,d]>0)&]],
(*else... A hyperplane constraint determines the d-th coordinate. *)
PointIntersection[FormVector[Part[hyp,First[nzx]]]]
]];

IntersectionData = MakeIntersectionData /@ PartialSearchPolyhedra;

An example:

IntersectionData[[5]]

How NextCoordinateInterval works

NextCoordinateInterval[partial_] := MakeSection[IntersectionData[[1+Length[partial]]], partial];

MakeSection[PointIntersection[vector_],base_] := With[{v = vector . Append[base,1]},{v,v}];

MakeSection[IntervalIntersection[lowbounds_, highbounds_],base_] :=With[{base1=Append[base,1]},{Max@@(lowbounds.base1),Min@@(highbounds.base1)}];

Let's Try It

{compTime,MagicSquareList} = Timing[MagicSquarePartialList[{}]//Flatten]; compTime

Length[MagicSquareList]

M1 =First[MagicSquareList]

NameGrid/.M1[[1]]

Does a list of rules actually define a magic square? The function PolyhedronContainsPoint may be applied to the list of values.

IsItMagic[aSquare_]:= PolyhedronContainsPoint[SearchPolyhedron, CellList/.First[aSquare]]===True;

Here is a check that all the members of MagicSquareList are magic squares.

And @@ (IsItMagic /@ MagicSquareList)