880 magic squares of order 4 by Frenicle and Fitting​
​by Oliver Seipel
The 880 magic squares of order 4 found by Bernard Frenicle de Bessy (before 1693) proof by Friedrich (and Hans) Fitting 1931.
References:
https://mathshistory.st-andrews.ac.uk/Biographies/Frenicle_de_Bessy/​
​https://en.wikipedia.org/wiki/Bernard_Fr% C3 % A9nicle_de _Bessy
"Rein mathematische Behandlung des Problems der lateinischen Quadraten von 16 und von 64 Feldern", Jahresbericht DMV 1931http://www.digizeitschriften.de/dms/img/?PID=GDZPPN002129337​
https://gallica.bnf.fr/ark:/12148/bpt6k5493994j​​
Reference: https://mathshistory.st-andrews.ac.uk/Biographies/Fitting/
Hans Fitting's father, Prof. Dr. Friedrich Fitting (1862-1945), was a graduate secondary school teacher and a research mathematician who published over twenty papers. His doctorate was from the Martin-Luther University of Halle-Wittenberg in 1888 for the thesis "Über eine Klasse von Berührungstransformationen"....
Friedrich Fitting is best known today for giving a proof, in 1931, that there are exactly 880 magic squares of order 4. This result appears in his paper "Rein mathematische Behandlung des Problems der magischen Quadrate von 16 und von 64 Feldern" and is remarkable since these 880 magic squares had been given by Frenicle de Bessy in 1693 but no proof was found until Friedrich Fitting's 1931 paper appeared.

The construction all 880 magic squares of order 4 by F. Fitting 1930, with contribution from H. Fitting.

Mathematica code:

In[]:=
SemiMagicSquareQ[m_?SquareMatrixQ]:=Equal@@Join[Apply[Plus,m,1],Apply[Plus,Transpose[m],1]]
In[]:=
MagicSquareQ[m_?SquareMatrixQ]:=Equal@@Join[Apply[Plus,m,1],Apply[Plus,Transpose[m],1],{Tr[m],Tr[Reverse/@m]}]
In[]:=
PanDiagonalMagicSquareQ[m_?SquareMatrixQ]:=With[{d=Length[m]},​​Equal@@Join[Apply[Plus,m,1],Apply[Plus,Transpose[m],1],Tr/@NestList[RotateLeft,m,d-1],Tr/@NestList[RotateLeft,Reverse/@m,d-1]]]
https://en.wikipedia.org/wiki/Fr% C3 % A9nicle_standard _form
In[]:=
FrenicleNormalFormQ[m_?SquareMatrixQ]:=m[[1,1]]==Min[m[[1,1]],m[[-1,1]],m[[-1,-1]],m[[1,-1]]]&&m[[2,1]]>m[[1,2]]
In[]:=
ToFrenicleNormalForm[m_?SquareMatrixQ]:=With[{a=m〚1,1〛,b=m〚-1,1〛,b12=m〚-1,2〛,b21=m〚-2,1〛,c=m〚1,-1〛,c12=m〚1,-2〛,c21=m〚2,-1〛,d=m〚-1,-1〛,d12=m〚-1,-2〛,d21=m〚-2,-1〛},With[{min=Min[a,b,c,d]},Which[​​mind&&d12>d21,Reverse[Map[Reverse,m]],​​mind&&d12≤d21,Reverse[Map[Reverse,m]],​​minc&&c12>c21,Reverse[m],​​minc&&c12≤c21,Map[Reverse,m],​​minb&&b12>b21,Map[Reverse,m],​​minb&&b12≤b21,Reverse[m],​​m〚1,2〛>m〚2,1〛,m,​​True,m]]]
In[]:=
EightSymmetries[m_?SquareMatrixQ]:={m,Reverse[Map[Reverse,m]],m,Map[Reverse,m],Reverse[Map[Reverse,m]],Reverse[m],Map[Reverse,m],Reverse[m]}

Decomposition of a magic square into binary components:

This is Duerer’s famous magic square of order 4 from his 1514 painting “Melancholia I”
​https://en.wikipedia.org/wiki/Melencolia_I
In[]:=
MagicSquareQ[Duerer1514={{15,2,1,12},{4,9,10,7},{8,5,6,11},{3,14,13,0}}]
Out[]=
True
The magic square starts here with 0, but adding 1 to all numbers gives the normal magic square:
In[]:=
Grid[{{Grid[Duerer1514,FrameAll,ItemSize{1,1.3},BackgroundLightOrange],Grid[Duerer1514+1,FrameAll,ItemSize{1,1.3},BackgroundLightOrange]}}]
Out[]=
Decompose an order 4 magic square (starting with 0) into binary components:
In[]:=
ToBinaryComponents[m_?SquareMatrixQ]:=Table[Map[IntegerDigits[#,2,4][[k]]&,m,{2}],{k,1,4}]
In[]:=
Grid[{Map[Grid,ToBinaryComponents[Duerer1514]]},FrameAll,BackgroundLightOrange]
Out[]=
1
0
0
1
0
1
1
0
1
0
0
1
0
1
1
0
1
0
0
1
1
0
0
1
0
1
1
0
0
1
1
0
1
1
0
0
0
0
1
1
0
0
1
1
1
1
0
0
1
0
1
0
0
1
0
1
0
1
0
1
1
0
1
0
Here every component is itself a binary magic square:
In[]:=
Map[MagicSquareQ,ToBinaryComponents[Duerer1514]]
Out[]=
{True,True,True,True}
And reconstruction from the binary components:
In[]:=
FromBinaryComponents[{a_,b_,c_,d_}]:=8a+4b+2c+d
In[]:=
Grid[1+FromBinaryComponents[ToBinaryComponents[Duerer1514]],FrameAll,ItemSize{1,1.3},BackgroundLightOrange]

I. All components are magic (528):

There can be 8 different magic components:
In[]:=
A1={{1,0,1,0},{0,1,0,1},{0,1,0,1},{1,0,1,0}};​​A2={{1,0,0,1},{0,1,1,0},{1,0,0,1},{0,1,1,0}};​​A3={{1,0,0,1},{1,0,0,1},{0,1,1,0},{0,1,1,0}};​​A4={{0,0,1,1},{1,1,0,0},{1,1,0,0},{0,0,1,1}};​​A5={{1,0,1,0},{1,0,1,0},{0,1,0,1},{0,1,0,1}};​​A6={{0,0,1,1},{1,1,0,0},{0,0,1,1},{1,1,0,0}};​​A7={{1,1,0,0},{1,0,1,0},{0,1,0,1},{0,0,1,1}};​​A8={{0,1,0,1},{1,1,0,0},{0,0,1,1},{1,0,1,0}};
In[]:=
Grid[Partition[MapIndexed[Labeled[Grid[#,FrameTrue,BackgroundLightOrange],
A
#2[[1]]
]&,{A1,A2,A3,A4,A5,A6,A7,A8}],4]]
Out[]=
1
0
1
0
0
1
0
1
0
1
0
1
1
0
1
0
A
1
1
0
0
1
0
1
1
0
1
0
0
1
0
1
1
0
A
2
1
0
0
1
1
0
0
1
0
1
1
0
0
1
1
0
A
3
0
0
1
1
1
1
0
0
1
1
0
0
0
0
1
1
A
4
1
0
1
0
1
0
1
0
0
1
0
1
0
1
0
1
A
5
0
0
1
1
1
1
0
0
0
0
1
1
1
1
0
0
A
6
1
1
0
0
1
0
1
0
0
1
0
1
0
0
1
1
A
7
0
1
0
1
1
1
0
0
0
0
1
1
1
0
1
0
A
8
They are indeed magic:
In[]:=
Map[MagicSquareQ,{A1,A2,A3,A4,A5,A6,A7,A8}]
Out[]=
{True,True,True,True,True,True,True,True}
And four of them are panmagic (A1, A2, A5, A6):
In[]:=
Map[PanDiagonalMagicSquareQ,{A1,A2,A3,A4,A5,A6,A7,A8}]
Out[]=
{True,True,False,False,True,True,False,False}
The complements are also magic:
In[]:=
Map[MagicSquareQ,1-{A1,A2,A3,A4,A5,A6,A7,A8}]
Out[]=
{True,True,True,True,True,True,True,True}
The complements are also pan magic if the original is panmagic:
In[]:=
Map[PanDiagonalMagicSquareQ,1-{A1,A2,A3,A4,A5,A6,A7,A8}]
Out[]=
{True,True,False,False,True,True,False,False}
In[]:=
Grid[Partition[MapIndexed[Labeled[Grid[1-#,FrameTrue,BackgroundLightOrange],1-
A
#2[[1]]
]&,{A1,A2,A3,A4,A5,A6,A7,A8}],4]]
Out[]=
0
1
0
1
1
0
1
0
1
0
1
0
0
1
0
1
1-
A
1
0
1
1
0
1
0
0
1
0
1
1
0
1
0
0
1
1-
A
2
0
1
1
0
0
1
1
0
1
0
0
1
1
0
0
1
1-
A
3
1
1
0
0
0
0
1
1
0
0
1
1
1
1
0
0
1-
A
4
0
1
0
1
0
1
0
1
1
0
1
0
1
0
1
0
1-
A
5
1
1
0
0
0
0
1
1
1
1
0
0
0
0
1
1
1-
A
6
0
0
1
1
0
1
0
1
1
0
1
0
1
1
0
0
1-
A
7
1
0
1
0
0
0
1
1
1
1
0
0
0
1
0
1
1-
A
8
Fitting shows then that there are 11 combinations with all four components magic:
In[]:=
eleven={{a1,a2,a3,a4},{a1,a2,a5,a6},{a1,a4,a7,a8},{a2,a3,a7,a8},{a1,a2,a3,a5},{a1,a3,a4,a5},{a2,a3,a4,a6},{a3,a4,a5,a6},{a1,a2,a4,a6},{a1,a4,a5,a6},{a2,a3,a5,a6}};
And there are sixteen possibilities to replace a component A with its complement 1-A:
In[]:=
sixteen[{a_,b_,c_,d_}]:=Tuples[{a1-a,b1-b,c1-c,d1-d}]
This leads to 528 magic squares (with magic components)
In[]:=
magic528=Select[Flatten[Table[Map[FromBinaryComponents[#]+1&,Permutations[k]],{k,Flatten[Map[sixteen,eleven],1]/.{a1A1,a2A2,a3A3,a4A4,a5A5,a6A6,a7A7,a8A8}}],1],FrenicleNormalFormQ];
There are 528 such squares:
In[]:=
Length[magic528]
Out[]=
528
They are all magic:
In[]:=
And@@Map[MagicSquareQ,magic528]
Out[]=
True
They are all in Frenicle normal form:
In[]:=
And@@Map[FrenicleNormalFormQ,magic528]
Out[]=
True
48 of them are pandiagonal magic squares:
In[]:=
Grid[Map[Grid[#,FrameAll,ItemSize{1,1.3},BackgroundLightBlue]&,Partition[Select[magic528,PanDiagonalMagicSquareQ],6],{2}]]

II. There are also semimagic components possible

Fitting gives two examples of magic squares of order 4 with semimagic components. I give here only the first example:
In[]:=
MagicSquareQ[hassemimagiccomponents1={{15,0,7,8},{2,11,12,5},{9,6,1,14},{4,13,10,3}}]
Out[]=
True
In[]:=
Grid[hassemimagiccomponents1,FrameAll,ItemSize{1,1.3},BackgroundLightOrange]
In[]:=
Grid[{Map[Grid,ToBinaryComponents[hassemimagiccomponents1]]},FrameAll,BackgroundLightOrange]
Out[]=
1
0
0
1
0
1
1
0
1
0
0
1
0
1
1
0
1
0
1
0
0
0
1
1
0
1
0
1
1
1
0
0
1
0
1
0
1
1
0
0
0
1
0
1
0
0
1
1
1
0
1
0
0
1
0
1
1
0
1
0
0
1
0
1
Only the first component is magic:
In[]:=
Map[MagicSquareQ,ToBinaryComponents[hassemimagiccomponents1]]
Out[]=
{True,False,False,False}
But all components are semi magic squares:
In[]:=
Map[SemiMagicSquareQ,ToBinaryComponents[hassemimagiccomponents1]]
Out[]=
{True,True,True,True}
Fitting introduced therefore new semimagic binary squares
B
1
,
C
1
,
C
2
:
In[]:=
B1={{1,0,0,1},{0,1,0,1},{1,0,1,0},{0,1,1,0}};​​C1={{0,0,1,1},{0,0,1,1},{1,1,0,0},{1,1,0,0}};​​C2={{0,1,0,1},{0,0,1,1},{1,1,0,0},{1,0,1,0}};
In[]:=
Map[SemiMagicSquareQ,{B1,C1,C2}]
Out[]=
{True,True,True}
In[]:=
Map[MagicSquareQ,{B1,C1,C2}]
Out[]=
{False,False,False}
In[]:=
Grid[MapIndexed[Labeled[Grid[#,BackgroundLightOrange,FrameTrue],If[#2[[1]]1,
B
1
,
C
#2[[1]]-1
]]&,{B1,C1,C2}]]
Out[]=
1
0
0
1
0
1
0
1
1
0
1
0
0
1
1
0
B
1
0
0
1
1
0
0
1
1
1
1
0
0
1
1
0
0
C
1
0
1
0
1
0
0
1
1
1
1
0
0
1
0
1
0
C
2
And three transformations for the type B:
In[]:=
Transform1B[m_?SquareMatrixQ]:=m[[{1,3,2,4}]][[{1,3,2,4}]]​​Transform2B[m_?SquareMatrixQ]:=m[[{2,1,4,3}]][[{2,1,4,3}]]​​Transform3B[m_?SquareMatrixQ]:=m[[{3,1,4,2}]][[{3,1,4,2}]]
The components (cf. above) can then be constructed from
A
2
,
B
1
,
C
1
,
C
2
In[]:=
Grid[{Map[Grid,{A2,EightSymmetries[Transform1B[B1]][[6]],EightSymmetries[Transform3B[B1]][[3]],EightSymmetries[Transform3B[C1]][[4]]}]},FrameAll,BackgroundLightOrange]
Out[]=
1
0
0
1
0
1
1
0
1
0
0
1
0
1
1
0
1
0
1
0
0
0
1
1
0
1
0
1
1
1
0
0
1
0
1
0
1
1
0
0
0
1
0
1
0
0
1
1
1
0
1
0
0
1
0
1
1
0
1
0
0
1
0
1
Compare to:
In[]:=
Grid[{Map[Grid,ToBinaryComponents[hassemimagiccomponents1]]},FrameAll,BackgroundLightOrange]
Out[]=
1
0
0
1
0
1
1
0
1
0
0
1
0
1
1
0
1
0
1
0
0
0
1
1
0
1
0
1
1
1
0
0
1
0
1
0
1
1
0
0
0
1
0
1
0
0
1
1
1
0
1
0
0
1
0
1
1
0
1
0
0
1
0
1
a) only one component is a B-Form: (192)
Fitting shows, that only some combinations result in a magic square:
Exact one is
B
1
and one is either
C
1
or
C
2
and two are
A
4
and
A
6
or
1-
A
4
and
1-
A
6
In[]:=
oneBForm={{a4,a6,b1,c1},​​{a4,b1,c1,a6},​​{b1,c1,a4,a6}};​​oneBForm=Join[oneBForm,oneBForm/.{c1c2}];​​oneBForm=Join[oneBForm,oneBForm/.{a4a6,a6a4}];​​oneBForm=Join[oneBForm,oneBForm/.{a41-a4,a61-a6},oneBForm/.{a41-a4},oneBForm/.{a61-a6}];​​oneBForm=oneBForm/.{a4A4,a6A6,b1B1,c1C1,c2C2};
In[]:=
Length[oneBForm]
Out[]=
48
Here is the n-th form with a single B-form: (there are 48, so
1≤n≤48
)
In[]:=
With[{n=6},Grid[{Map[Grid[#,ItemSizeAll]&,oneBForm[[n]]]},FrameAll,BackgroundLightOrange]]
Out[]=
1
0
0
1
0
1
0
1
1
0
1
0
0
1
1
0
0
1
0
1
0
0
1
1
1
1
0
0
1
0
1
0
0
0
1
1
1
1
0
0
1
1
0
0
0
0
1
1
0
0
1
1
1
1
0
0
0
0
1
1
1
1
0
0
This leads to 48 magic squares:
In[]:=
b48=Map[ToFrenicleNormalForm,Flatten[Table[Map[FromBinaryComponents[#]+1&,{k}],{k,oneBForm}],1]];
In[]:=
And@@Map[MagicSquareQ,b48]
Out[]=
True
And now he uses the three transforms from above to quadrupple this number:
In[]:=
Length[b192=Join[b48,Map[ToFrenicleNormalForm,Map[Transform1B,b48]],Map[ToFrenicleNormalForm,Map[Transform2B,b48]],Map[ToFrenicleNormalForm,Map[Transform3B,b48]]]]
Out[]=
192
In[]:=
And@@Map[MagicSquareQ,b192]
Out[]=
True
b) Two components are B-Forms: (112)
Fitting shows, that there can be 2 semimagic B forms
(
B
2
to
B
6
)
In[]:=
B2={{1,0,0,1},{1,0,1,0},{0,1,0,1},{0,1,1,0}};​​B3={{0,1,0,1},{1,0,0,1},{0,1,1,0},{1,0,1,0}};​​B4={{0,1,0,1},{0,1,1,0},{1,0,0,1},{1,0,1,0}};​​B5={{0,0,1,1},{1,0,0,1},{0,1,1,0},{1,1,0,0}};​​B6={{0,0,1,1},{0,1,1,0},{1,0,0,1},{1,1,0,0}};
In[]:=
Map[SemiMagicSquareQ,{B2,B3,B4,B5,B6}]
Out[]=
{True,True,True,True,True}
In[]:=
Map[MagicSquareQ,{B2,B3,B4,B5,B6}]
Out[]=
{False,False,False,False,False}
In[]:=
Grid[MapIndexed[Labeled[Grid[#,BackgroundLightOrange,FrameTrue],
B
#2[[1]]+1
]&,{B2,B3,B4,B5,B6}]]
Out[]=
1
0
0
1
1
0
1
0
0
1
0
1
0
1
1
0
B
2
0
1
0
1
1
0
0
1
0
1
1
0
1
0
1
0
B
3
0
1
0
1
0
1
1
0
1
0
0
1
1
0
1
0
B
4
0
0
1
1
1
0
0
1
0
1
1
0
1
1
0
0
B
5
0
0
1
1
0
1
1
0
1
0
0
1
1
1
0
0
B
6
Fitting shows, if there are two B forms, then there must be one C form and one A form (the details are quite tricky):
In[]:=
twoBForms={{b1,c1,b2,a4},​​{b1,c1,b3,a4},{b1,c1,b3,a6},{b1,c1,b4,a6},{b1,c2,b2,a1},{b1,c2,b5,a8},{b1,c2,b6,a8}};​​twoBForms=Join[twoBForms,twoBForms/.{a11-a1,a41-a4,a61-a6,a81-a8}];​​twoBForms=twoBForms/.{a1A1,a4A4,a6A6,a8A8,b1B1,c1C1,c2C2,b2B2,b3B3,b4B4,b5B5,b6B6};
In[]:=
Length[twoBForms]
Out[]=
14
In[]:=
bb1=Map[ToFrenicleNormalForm,Flatten[Table[Map[(*{d,c,b,a}*){8,2,4,1}.#+1&,{k}],{k,twoBForms}],1]];
In[]:=
Length[bb1]
Out[]=
14
In[]:=
bb2=Map[ToFrenicleNormalForm,Flatten[Table[Map[(*{d,c,b,a}*){4,1,2,8}.#+1&,{k}],{k,twoBForms}],1]];
In[]:=
Length[bb2]
Out[]=
14
In[]:=
bb28=Join[bb1,bb2];
In[]:=
Length[bb28]
Out[]=
28
And now he uses the three transforms from above to quadrupple this number:
In[]:=
428
Out[]=
112
In[]:=
bb112=Map[ToFrenicleNormalForm,Join[bb28,Map[Transform1B,bb28],Map[Transform2B,bb28],Map[Transform3B,bb28]]];​​Length[bb112]
Out[]=
112
In[]:=
And@@Map[MagicSquareQ,bb112]
Out[]=
True
c) Three components are B-Forms: (48)
Fitting shows, if there are three B forms, then there must be one C form (the details are quite tricky):
In[]:=
threeBForms={{b1,c1,b2,b3},​​{b1,c1,b2,b4},{b1,c1,b3,b4},{b1,c2,b2,b5},{b1,c2,b2,b6},{b1,c2,b5,b6}};​​threeBForms=threeBForms/.{b1B1,c1C1,c2C2,b2B2,b3B3,b4B4,b5B5,b6B6};
In[]:=
bbb1=Map[ToFrenicleNormalForm,Flatten[Table[Map[(*d,c,b,a*){8,1,4,2}.#+1&,{k}],{k,threeBForms}],1]];Length[bbb1]
Out[]=
6
In[]:=
bbb2=Map[ToFrenicleNormalForm,Flatten[Table[Map[(*d,c,b,a*){8,1,2,4}.#+1&,{k}],{k,threeBForms}],1]];Length[bbb2]
Out[]=
6
In[]:=
bbb12=Join[bbb1,bbb2];
In[]:=
Length[bbb12]
Out[]=
12
And now he uses the three transforms from above to quadrupple this number:
In[]:=
412
Out[]=
48
In[]:=
bbb48=Map[ToFrenicleNormalForm,Join[bbb12,Map[Transform1B,bbb12],Map[Transform2B,bbb12],Map[Transform3B,bbb12]]];
In[]:=
Length[bbb48]
Out[]=
48
In[]:=
And@@Map[MagicSquareQ,bbb48]
Out[]=
True
Fitting then shows that there are no more magic squares of order 4 possible. So there are 880 different magic squares of order 4 in Frenicle normal form, as was conjectured by Frenicle before 1693.
In[]:=
fitting880=Join[bbb48,bb112,b192,magic528];
In[]:=
Length[fitting880]
Out[]=
880
In[]:=
And@@Map[MagicSquareQ,fitting880]
Out[]=
True
Here are all 880 magic squares of order 4 in Frenicle normal form in groups of 40:
The pandiagonal magic squares are with lightblue background (401 to 448) the others with lightorange background:
In[]:=
Manipulate[With[{m=40},Labeled[Grid[Map[Grid[#,BackgroundIf[PanDiagonalMagicSquareQ[#],LightBlue,LightOrange],FrameAll,ItemSize{1,1.3}]&,Partition[Take[fitting880,{mi+1,m(i+1)}],5],{2}]],Row[{"Frénicle's and Fitting's magic squares of order 4, from ",mi+1," to ",m(i+1)}],Top]],{{i,0,"40 squares each"},0,21,1},SaveDefinitionsTrue]
Out[]=
If the condition that a magic square is in Frenicle normal form is released, then there are 8 times as much magic squares.
Here are all 7040 possible magic squares of order 4 in 880 groups of 8 (only the first is in Frenicle normal form), the pandiagonal magic squares are with lightblue background (401 to 448), the others with lightorange background:
In[]:=
8880
Out[]=
7040
In[]:=
Manipulate[Labeled[Grid[Partition[Map[Grid[#,BackgroundIf[PanDiagonalMagicSquareQ[#],LightBlue,LightOrange],FrameAll,ItemSize{1,1.3}]&,EightSymmetries[fitting880[[i]]]],4]],Row[{i,". magic square (from 880) in 8 similar forms."}],Top],{{i,1,"n-th magic square of order 4 in 8 similar forms."},1,880,1},SaveDefinitionsTrue]
Out[]=