Order 4 allows to compute easily the set of magic (or semimagic) squares. Then you can handle this set for examination and you have a lot of results on the structure. For orders superior to 4, it is practically impossible to handle the set: the squares are too numerous or can't be computed.
The 32 transformations of a given magic square of order 4 can be obtained as follows:
 you take first the 8 classical transformations of the octic group of the square: I, V, H, G, D, R1, R2, R3 (the notations are arbitrary: I=Identity, V=Vertical, H=Horizontal, G=left or Gauche in French, D=right or Droite in French, R1 = rotation of 90° anticlockwise, R2=symetry from the centre, R3=rotation of 90° clockwise)
 then you multiply these 8 transformations by the transformation IT (see hereafter; IT=permutation of interior rows and columns, that is to say in the order 1324): you get 16 transformations,
 after that, youmultiply these 16 transformations by the transformation A (see hereafter; A=permutation in the order 2143 of rows and columns).
With the matricial notation, these 32 transformations are:
A1  A2  A3  A4 
B1  B2  B3  B4 
C1  C2  C3  C4 
D1  D2  D3  D4 
I
 
A4  A3  A2  A1 
B4  B3  B2  B1 
C4  C3  C2  C1 
D4  D3  D2  D1 
V
 
D1  D2  D3  D4 
C1  C2  C3  C4 
B1  B2  B3  B4 
A1  A2  A3  A4 
H
 
A1  B1  C1  D1 
A2  B2  C2  D2 
A3  B3  C3  D3 
A4  B4  C4  D4 
G


D4  C4  B4  A4 
D3  C3  B3  A3 
D2  C2  B2  A2 
D1  C1  B1  A1 
D
 
A4  B4  C4  D4 
A3  B3  C3  D3 
A2  B2  C2  D2 
A1  B1  C1  D1 
R1
 
D4  D3  D2  D1 
C4  C3  C2  C1 
B4  B3  B2  B1 
A4  A3  A2  A1 
R2
 
D1  C1  B1  A1 
D2  C2  B2  A2 
D3  C3  B3  A3 
D4  C4  B4  A4 
R3


A1  A3  A2  A4 
C1  C3  C2  C4 
B1  B3  B2  B4 
D1  D3  D2  D4 
IT
 
D4  D2  D3  D1 
B4  B2  B3  B1 
C4  C2  C3  C1 
A4  A2  A3  A1 
EX
 
A1  C1  B1  D1 
A3  C3  B3  D3 
A2  C2  B2  D2 
A4  C4  B4  D4 
M
 
D4  B4  C4  A4 
D2  B2  C2  A2 
D3  B3  C3  A3 
D1  B1  C1  A1 
N


A4  A2  A3  A1 
C4  C2  C3  C1 
B4  B2  B3  B1 
D4  D2  D3  D1 
X
 
D1  D3  D2  D4 
B1  B3  B2  B4 
C1  C3  C2  C4 
A1  A3  A2  A4 
Y
 
A4  C4  B4  D4 
A2  C2  B2  D2 
A3  C3  B3  D3 
A1  C1  B1  D1 
Z
 
D1  B1  C1  A1 
D3  B3  C3  A3 
D2  B2  C2  A2 
D4  B4  C4  A4 
T


B2  B1  B4  B3 
A2  A1  A4  A3 
D2  D1  D4  D3 
C2  C1  C4  C3 
A
 
B3  B4  B1  B2 
A3  A4  A1  A2 
D3  D4  D1  D2 
C3  C4  C1  C2 
V*A
 
C2  C1  C4  C3 
D2  D1  D4  D3 
A2  A1  A4  A3 
B2  B1  B4  B3 
H*A
 
B2  A2  D2  C2 
B1  A1  D1  C1 
B4  A4  D4  C4 
B3  A3  D3  C3 
G*A


C3  D3  A3  B3 
C4  D4  A4  B4 
C1  D1  A1  B1 
C2  D2  A2  B2 
D*A
 
B3  A3  D3  C3 
B4  A4  D4  C4 
B1  A1  D1  C1 
B2  A2  D2  C2 
R1*A
 
C3  C4  C1  C2 
D3  D4  D1  D2 
A3  A4  A1  A2 
B3  B4  B1  B2 
R2*A
 
C2  D2  A2  B2 
C1  D1  A1  B1 
C4  D4  A4  B4 
C3  D3  A3  B3 
R3*A


C3  C1  C4  C2 
A3  A1  A4  A2 
D3  D1  D4  D2 
B3  B1  B4  B2 
IT*A
 
B2  B4  B1  B3 
D2  D4  D1  D3 
A2  A4  A1  A3 
C2  C4  C1  C3 
EX*A
 
C3  A3  D3  B3 
C1  A1  D1  B1 
C4  A4  D4  B4 
C2  A2  D2  B2 
M*A
 
B2  D2  A2  C2 
B4  D4  A4  C4 
B1  D1  A1  C1 
B3  D3  A3  C3 
N*A


C2  C4  C1  C3 
A2  A4  A1  A3 
D2  D4  D1  D3 
B2  B4  B1  B3 
X*A
 
B3  B1  B4  B2 
D3  D1  D4  D2 
A3  A1  A4  A2 
C3  C1  C4  C2 
Y*A
 
C2  A2  D2  B2 
C4  A4  D4  B4 
C1  A1  D1  B1 
C3  A3  D3  B3 
Z*A
 
B3  D3  A3  C3 
B1  D1  A1  C1 
B4  D4  A4  C4 
B2  D2  A2  C2 
T*A

These notations are personal notations.
You have a set of 32 elements:
1 element of period 1: I,
19 elements of period 2: V, H, G, D, R2, IT, EX, M, N, X, Y, A, V*A, H*A, G*A, D*A, R2*A, Z*A, T*A,
12 elements of period 4 (all with a common square R2): R1, R3, Z, T, R1*A, R3*A, IT*A, EX*A, M*A,
N*A, X*A, Y*A.
It is possible to build a multiplication table and to demonstrate in that way you have a group.
The relations of definition
("abstract definition") allows a smarter demonstration. You have nevertheless to proceed by trial and error to find the good relations of definition among the 51 possible groups of order 32! The literature does not give exhaustive catalogues, only partial catalogues and for low orders.
When you take
s_{1} = R1, s_{2} = X*A, s_{3} = IT, s_{4} = G,
you find that all the relations of definition of group G 32 n° 7 of BURNS are verified [Reference: Josephine E. BURNS,
The Abstract Definitions of Groups of Degree 8, Am. J. Math. 37 (1915), p. 207]:
s_{1}^{4} = s_{2}^{4} = s_{3}^{2} = s_{4}^{2} = s_{1}^{3} s_{2 }s_{1 }s_{2} = s_{3 }s_{1 }s_{3 }s_{1}^{3} = s_{2}^{3} s_{1 }s_{2 }s_{1 }
= (s_{4 }s_{1})^{2} = (s_{4 }s_{2})^{2} = (s_{4 }s_{3})^{2} = (s_{3 }s_{2})^{2} = 1.

Besides, you verify that the 4 generators R1, X*A, IT et G give 32 elements.
Therefore, the 32 transformations form a group:
the group G 32 n° 7 of BURNS.
The 8 first transformations form the octic group D4.
The 16 first transformations form the group D4*C2. It is the group n° 4 of BURNSIDE among the 9 not abelian groups of order 16 [Reference: W. BURNSIDE,
Theory of Groups of finite order, Dover 1911, p.146]:
P^{4} = E, Q^{2} = E, R^{2} = E, R^{1}PR = P^{3}, P^{1}QP = Q, R^{1}QR = Q
You can verify it easily with P = R1, Q = IT, R = G
It is the direct product of {Q} by {P,R}.
We can calculate the order of the group G for magic squares of order n.
For magic squares of order 4, we have seen that group G is of order 32. It is the product of 8 (order of group octic) by 4 (number of remaining transpositions).
In magic square language, we call
transpositions permutations of rows and columns symmetrically located from the centre of the square (it is not the usual mathematical language for permutations). If the rows are called 1234, we see there are 8 transpositions:
We search permutations abcd such as a+d=b+c=5
It is sufficient to fix a and b: a can take 4 values, and b=2 values when a is given.
Then we have 4*2=8 solutions
Among these transpositions, the transposition 4321 is already in the octic group: it is the transformation R2 (symmetry about the centre).
Therefore, 32 is the product:
8*(number of transpositions /2) = 4*(number of transpositions)
This formula can be generalized for all the orders.
For order 6, the number of transpositions abcdef is 6*4*2=48
a can take 6 values, b can take 4 values when a is given, c can take 2 values when a and b are given, d=7c, e=7b, f=7a.
Then the order of group G is 4*48 = 192.
For order 2p, the number of transpositions is 2p*(2p2)*...*2=2
^{p }(p!)
and the order of group G is 4*2
^{p }(p!) = 2
^{p+2} (p!).
We have the fundamental result:
for magic squares of order 2p, the order of group G is 2^{p+2 }(p!)
We see easily that
for magic squares of order 2p +1, the order of group G is the same: 2^{p+2 }(p!).
As a matter of fact, the number of transpositions is the same as for order 2p (the middle number doesn't move).
These two groups are isomorphic; they have the same abstract definition.
Therefore,
for magic squares of order n=2p or n=2p+1, the order of group G is 2^{p+2} (p!).
This last formula seems original. The first values are:
Order of squares  Order of group G 
3  8 
4  32 
5  32 
6  192 
7  192 
8  1 536 
Note:
for semimagic of order n, the order of group G is 2 (n!)^{2}.
The demonstration is here obvious: for a given square, you have n! possible permutations of rows, idem for columns, therefore (n!)^{2 }possible transformations and you double this number with symmetry about the first diagonal.
This set is defined for example by A1+A2+B1+B2=34 (it is not the only definition: you can have also A2+B1+C4+D3=34 or A1+A3+C1+C3=34, etc.).
3.1 Group G
The group G (transformations) is of order 1 152; you have 3 elementary squares.
I identified this group: it is the only group of order 1 152 of degree 8 (reference: J. BURNS's article already quoted). Its abstract definition is:
s_{1}^{12} = s_{2}^{2} = (s_{1}^{4} s_{2})^{2}(s_{1}^{8} s_{2})^{ 2} = (s_{1}^{3} s_{2})^{2}(s_{1}^{9} s_{2})^{ 2} = (s_{1}^{3} s_{2 }s_{1}^{4} s_{2})^{ 2} = 1

with for example
s_{1} = 
D1  A4  D4  A1 
C2  B3  C3  B2 
D2  A3  D3  A2 
C1  B4  C4  B1 


s_{2} = IT = 
A1  A3  A2  A4 
C1  C3  C2  C4 
B1  B3  B2  B4 
D1  D3  D2  D4 

You can verify that all the relations of definition are true with these two generators (and you need to verify also that you have 1 152 generated elements: these relations aren't abstract definition of a subgroup).
This group has a large number of subgroups (non exhaustive list, notations of J. BURNS's article if any else indication)  it is a mine! :
G 576 n° 1, 2, 3
G 288 n° 3
G 192 n° 1, 2
G 144 n° 1
G 128
G 96 n° 3, 4, 6
G 72 degree 6
G 64 n° 2, 3, 5
G 48 n° 1
G 32 n° 6, 7, 8, 9
G 24 degree 4
G 16 n° 1 (?), 9
G 8 (octic group) and its subgroups.
For each subgroup, the abstract definition is in the literature and I can give example of generators.
I think G 1 152 has not subgroup of order 384.
3.2 Group G'
The group G' (permutations) is of order 384; you have 9 elementary squares.
This group is G 384 degree 8.
Subgroups (non exhaustive list):
G 192 n° 4, 5
G 128 and its subgroups (the same as above).
G 1 152 and G 384 have G 128 as a common subgroup (with also all its own subgroups).
3.3 Decomposition
The set of 3 456 semipandiagonal squares is made with 9 types of squares:





Type I


Type II


Type III






Type IV


Type V


Type VI






Type IV'


Type V'


Type VI'

Note: I put a double frame in a type diagram and a simple frame in a transformation diagram.
Each type has 384 squares and a group structure.
There are the 6 DUDENEY's types I, II, III, IV, V, VI (IV and IV' are the same type for DUDENEY; idem for V and V'; idem also for VI and VI').
It is funny to identify the transformations between the 9 types of squares (internal transformations for some types, but external transformations for other types). With these external transformations, you see that the 9 sets of 384 squares are isomorphic.
Note you have 2 supplementary isomorphic sets with 384 squares: they are of DUDENEY's type VI. These 768 squares aren't semipandiagonal.
They can be defined  for example  as squares of DUDENEY's type VI with condition A2+B4+C1+D3=34 or A3+B1+C4+D2=34.
You have 1 152 squares with condition A2+B4+C1+D3=34 or A3+B1+C4+D2=34: the 384 symmetrical squares and the 2 supplementary sets above.
3.4 Intermediate Square
All the 3 456 semipandiagonal squares can be built with only one intermediate square. For example, you can take (cf. BENSON & JACOBY, page 111):
A+a  B+b  C+c  D+d 
D+c  C+d  B+a  A+b 
B+d  A+c  D+b  C+a 
C+b  D+a  A+d  B+c 
This intermediate square is a GrecoLatin square: it is magic, diagonal (see definition in
The intermediate square method, § 4) and semipandiagonal.
You have no magic conditions. You have 4! =24 solutions for capital letters, 4! solutions too for lowercase letters, k=6 basis; therefore you have a set of 6*24*24 = 3 456 squares.
This demonstration is simplier than BENSON & JACOBY'.
For magic squares of order 4, you find the following graph of filiation:
K is the number of magic series (relations a+b+c+d = 34). Cf.
The transformation method, § 5.1.
The origin set of 27 648 squares is the set of semimagic squares of DUDENEY's type VI, i.e. defined by the 3 conditions:
A1+D1=A2+D2=A3+D3=17 




Type VI 
The set n° 2 of 1 216 squares is the set of 27 648 semimagic squares with condition:
A1+B2+C3+D4=34 (and D1+C2+B3+A4=34 which is a consequence of DUDENEY's type VI)
The set n° 1 of 384 squares is the set n° 2 of 1 216 squares with condition of semipandiagonal squares:
A1+A2+B1+B2=34 for example (and other induced relations).
This set n° 1 is the set of semipandiagonal squares of DUDENEY's type VI.
The set of 3 072 squares is the set of semimagic squares of DUDENEY's type IX.

Type IX 
It is also the set of 27 648 semimagic squares with condition:
D1+B2+B3+D4=34 (and other induced relations) given by transformation
A1  A2  A3  A4 
D1  B2  B3  D4 
C1  C2  C3  C4 
B1  D2  D3  B4 
The set n° 3 of 224 squares is the set of 3 072 squares with magic conditions for the diagonals (and with other induced conditions).
For order 4, it is possible to show exactly how the 7 040 magic squares are generated from capital and lowercase letter squares. For that, it is not necessary to enumerate all the possible Gallic squares, you need only to identify the elementary lowercase letter Gallic squares (and the induced squares by group G of 32 and by permutations). You search after the orthogonal squares to these elementary Gallic squares.
You have to search only Gallic squares with semiregular rows and columns because only diagonals can be irregular (you verify this property a posteriori on all the 7 040 squares and all the basis; maybe it is possible to demonstrate directly that rows and columns cannot be irregular). If you classify the Gallic squares according to the value of the parameter S1, sum of the first diagonal, you verify on the 7 040 squares that S1 can take only some values (then the sum of the second diagonal is S2=2σS1). And among these Gallic squares, you have to take only those which have orthogonal squares.
I give hereafter the result of my research.
With the first basis (1,2,3,4;0,4,8,12), I found:

there are 9 elementary squares S1=10 in a normalized position (the numbers hereafter are those of my program of enumeration), and each one has orthogonal squares:
1) 
 2) 
 3) 
 4) 
 5) 


8) 
 9) 
 12) 
 13) 

(squares # 1, 2 and 12 are Latin and magic, square # 8 is Latin diagonal)

there are 7 elementary squares S1=6, but only 6 with orthogonal squares:
(squares # 1, 2 and 12 are Latin and magic, square # 8 is Latin diagonal)

there are 9 elementary squares S1=9, but only 7 with orthogonal squares:
(square # 11 is Latin but not magic).
For each elementary square, I searched the number of induced squares and the number of orthogonal squares:
S1=10  # of elementary square  Number of induced squares (group 32 * permutations)  Number of orthogonal squares for each  Total number of orthogonal squares 
 1  4*2= 8  16  8*16= 128 
2  4*2= 8  32  8*32= 256 
3  32*2=64  12  64*12= 768 
4  16*2=32  40  32*40=1 280 
5  32*1=32  8  32*8= 256 
8  16*3=48  24  48*24=1 152 
9  32*1=32  8  32*8= 256 
12  16*1=16  32  16*32= 512 
13  16*1=16  40  16*40= 640 
256  5 248 
(There are 80 Latin magic squares, coming from the elementary squares # 1, 2, 8 and 12.
There are 48 Latin diagonal squares coming from the elementary square # 8.
You see that Latin magic squares haven't the maximum number of orthogonal squares.
The total number of orthogonal squares to Latin magic squares is 2 048  with 1 152 orthogonal to Latin diagonal squares  ).
S1=6  # of elementary square  Number of induced squares (group 32 * permutations)  Number of orthogonal squares for each  Total number of orthogonal squares 
 1  16*1=16  16  16*16=256 
2  16*1=16  16  16*16=256 
4  32*1=32  12  32*12=384 
5  32*1=32  14  32*14=448 
8  32*1=32  6  32* 6=192 
18  16*1=16  16  16*16=256 
144  1 792 
S1=10  # of elementary square  Number of induced squares (group 32 * permutations)  Number of orthogonal squares for each  Total number of orthogonal squares 
 1  32*2=64  4  64*4=256 
3  32*2=64  4  64*4=256 
6  32*2=64  4  64*4=256 
7  32*1=32  6  32*6=192 
11  32*3=96  4  96*4=384 
13  32*2=64  3  64*3=192 
23  32*1=32  8  32*8=256 
416  1 792 
Then, I built two figures: one for the 5 248 semiregular squares and another for the 1 792 irregular squares:
The 7 040 magic squares are distributed in blocks of squares. In each block, all the squares come from the same couple of two elementary Gallic squares, each one in a normalized position.
In fact, for the capital letters values (0,4,8,12), I used the same program as for lowercase letters values (1,2,3,4) with:
0→1 4→2 8→3 12→4
In that way, the Gallic squares are written in a normalized form.
The figures give the number of magic squares in each block, but naturally I identified each magic square in the list of 7 040.
All the 7 040 magic squares can be built from 9+6+7=22 elementary Gallic squares (with the first basis).
I also drove the decomposition of the 7 040 squares with the second and the third basis (with the 3 other basis, the diagrams are obtained by permutation of capital and lowercase letters, i.e. by symmetry about the first bisecting line).
I give the results:
 Basis # 1  Basis # 2  Basis #3 
Number of semiregular squares  5 248 (with 1 152 regular)  4 224 (with 1 152 regular)  4 736 (with 1 152 regular) 
Number of irregular squares  1 792  2 816  2 304 
Number of elementary squares  22  32  45 
Number of blocks of squares  36  41  41 
The set of 1 152 regular squares is different in the 3 basis. The 3 sets of 1 152 constitute the set of 3 456 semipandiagonal squares.
There are 4 224 squares which are always semiregular in all the basis
and 1 280 squares which are always irregular in all the basis.
There are 5 760 semiregular squares with at least one basis.
Note: for orders superior to 4, it is naturally impossible to do the same task up to the very end.
It is possible to have other decompositions (of the 7 040 squares) than decomposition in types. I show here 2 decompositions in sets of intermediate squares.
6.1 Decomposition with a given basis
You can make a review of the 7 040 squares with a given basis, the first one for example.
Each square generates an intermediate square, and this intermediate square generates a set of squares by tabulation with the first basis. The 7 040 squares can then be distributed in sets of disjointed intermediate squares. For each set, you can define an elementary square and search isomorphic sets (by transformations and permutations).
You find 9 elementary squares with the first basis (the numbers hereafter are those of my classical program of enumeration, which is a reduced/32 program):
1) 
1  8  12  13 
7  14  2  11 
10  3  15  6 
16  9  5  4 
 7) 
1  7  14  12 
8  13  2  11 
9  4  15  6 
16  10  3  5 
 13) 
1  14  8  11 
7  12  2  13 
10  5  15  4 
16  3  9  6 


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


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

# of square  Conditions on letters  Number of squares (with the intermediate square)  Number of isomorphic sets  Total number of squares  Remark 
13  no conditions  576  2  1 152  regular and semiregular 
1  A+D=B+C a+d=b+c  64  48  3 072  semiregular 
25  A+D=B+C a+d=b+c 2b=a+c  16  64  1 024  semiregular 
35  a+d=b+c A+2B+D+2c+2d=S  24  8  192  irregular 
39  a+d=b+c A+2B+D+b+3d=S  12  16  192  irregular 
33  A+D=B+C a+d=b+c A+2B+D+2c+2d=S  8  48  384  irregular 
17  A+D=B+C a+d=b+c A+B+2D+2a+2b=S  8  48  384  irregular 
7  A+D=B+C a+d=b+c A+B+2D+3a+c=S  4  80  320  irregular 
54  A+D=B+C a+d=b+c A+2B+D+b+3d=S  4  80  320  irregular 
 7 040 
You find that the numbers of squares in the 9 elementary sets are 576 and its divisors (64, 24, 16, 12, 8 and 4). You haven't all the divisors. With the second basis, you find also sets of 32 squares.
The 7 040 squares of order 4 can be decomposed in elementary sets of q squares, q being a divisor
of (n!)^{2}.
It is possible to have other decompositions with the 5 other basis.
6.2 Decomposition with all the basis at the same time
You can do the same at the beginning but you apply after the 6 basis to each intermediate square. Then you have to choose each elementary square in order to have disjointed sets with the 6 basis.
You find 6 elementary squares among the 9 previous elementary squares:
Sub set  # of square  Conditions on letters  Number of squares (with the intermediate square)  Total number of squares (6 basis)  Number of iso morphic sets  Total number of squares 
Basis #1  Basis #2  Basis #3  Basis #4  Basis #5  Basis #6 
I  13  no conditions  576  576  576  576  576  576  3 456  1  3 456 
II  1  A+D=B+C a+d=b+c  64  64  64  64  64  64  384  2  768 
III  25  A+D=B+C a+d=b+c 2b=a+c  16  16  0  0  0  16  48  32  1 536 
IV  33  A+D=B+C a+d=b+c A+2B+D+2c+2d=S  8  0  0  0  8  0  16  48  768 
V  54  A+D=B+C a+d=b+c A+2B+D+b+3d=S  4  0  4  0  4  0  12  32  384 
VI  7  A+D=B+C a+d=b+c A+B+2D+3a+c=S  4  0  4  4  0  4  16  8  128 
 7 040 
This decomposition is more "symmetrical".
You have the same property:
the 7 040 squares of order 4 can be decomposed in elementary sets of q squares,
q being a divisor of (n!)^{2}.
Finally, you have different ways to generate and to decompose the 7 040 magic squares of order 4:
7.1 Generations
7.2 Decompositions
Maybe, it should be possible to find other generations and decompositions.