(8 votes, average: 5.00 out of 5)

## Magic Squares

While working on an idea for maximizing a return on a financial portfolio using magic squares, we though we would share our Python code in case someone else may have some use for it.

```import string

def make_magic_square(n):
if n < 1 or n == 2:
return False
if n % 2 == 1:
return make_odd_magic_square(n)
else:
if n % 4 == 0:
return make_even_4n_magic_square(n)
elif (n-2) % 4 == 0:
return make_even_4np2_magic_square(n)
else:
return False

def make_panmagic_square(n):
if n < 4:
return False
if n % 2 == 0:
return make_even_panmagic_square(n)
else:
return make_odd_panmagic_square(n)

def make_odd_magic_square(n):
"""
Source: http://blogs.mathworks.com/cleve/2012/11/05/magic-squares-part-2-algorithms/
This method is called the Siam
Build I and J matrices as follows:

1 1 1 1 1		1 2 3 4 5
2 2 2 2 2		1 2 3 4 5
I=> 3 3 3 3 3	J=> 1 2 3 4 5
4 4 4 4 4		1 2 3 4 5
5 5 5 5 5		1 2 3 4 5
"""
if n < 1 or n%2 == 0: return False	#only allow odd squares 2n-1, n>0

J = [range(1, n+1)] * n
I = transpose(J)
"""
Using these indices, we generate two more matrices A and B
A = mod(I + J + (n-3)/2, n)
B = mod(I + 2*J - 2, n)

3 4 0 1 2		1 3 0 2 4
4 0 1 2 3		2 4 1 3 0
A=> 0 1 2 3 4	B=> 3 0 2 4 1
1 2 3 4 0		4 1 3 0 2
2 3 4 0 1		0 2 4 1 3
"""
A = [[(I[i][j] + J[i][j] + (n-3)/2) % n for i in range(n)] for j in range(n)]
B = [[(I[i][j] + 2*J[i][j] - 2) % n for j in range(n)] for i in range(n)]
"""
Both A and B are fledgling magic squares. They have equal row, column and
diagonal sums. But their elements are not the integers from 1 to n^2. Each
has duplicated elements between 0 and n-1. The final statement:

M = n*A + B + 1

produces a matrix whose elements are integers between 1 and n^2 and which
has equal row, column and diagonal sums. What is not obvious, is that
there are no duplicates and includes all of the integers between 1 and n^2
making it a magic square.
"""
return [[n*A[i][j] + B[i][j] + 1 for j in range(n)] for i in range(n)]

def make_even_4np2_magic_square(n):
"""
Source: http://www.math.wichita.edu/~richardson/mathematics/
magic%20squares/even-ordermagicsquares.html
The algorithm starts by partitioning the square into 4 blocks of n/2 each;
one for each quadrant.  We 	fill those blocks with a known odd sized magic
square, A,  as: A in Q II, A + (n/2)^2 	in Q III, A + 2*(n/2)^2 in Q I, and
A + 3*(n/2)^2 in Q IV.

This maintains the proper column totals, but not the row or diagonal totals.
Correcting row totals is achieved by swaping selective columns (0, 1, 9) of the upper
10x10 matrix with the lower half.  The columns affected are 0, 1, 2, ..., k
and n-1, n-2, ... n-k.  (k=(n-2)/4)
Example (zero based):
n=6	  columns 0
n=10  columns 0, 1, 9
n=14  columns 0, 1, 2, 12, 13
n=18  columns 0, 1, 2, 3, 15, 16, 17

The final step is to switch 2 values in certain columns to correct the diagonals
and make a magic square. (nx-k-1, [k-1, k]) <-> (nx+k, [k-1, k])
Example (zero based):
n=6		(1, 0) <-> (4, 0), (1, 1) <-> (4, 1)
n=10	(2, 1) <-> (7, 1), (2, 2) <-> (7, 2)
n=14	(3, 2) <-> (10, 2), (3, 3) <-> (10, 3)
n=18	(4, 3) <-> (13, 3), (4, 4) <-> (13, 4)

Fill in 4 Quads with 5x5 M.S.               Swap upper and lower column sections
17  24   1   8  15  67  74  51  58  65			92  99   1   8  15  67  74  51  58  40
23   5   7  14  16  73  55  57  64  66			98  80   7  14  16  73  55  57  64  41
4   6  13  20  22  54  56  63  70  72			79  81  13  20  22  54  56  63  70  47
10  12  19  21   3  60  62  69  71  53			85  87  19  21   3  60  62  69  71  28
11  18  25   2   9  61  68  75  52  59			86  93  25   2   9  61  68  75  52  34
92  99  76  83  90  42  49  26  33  40			17  24  76  83  90  42  49  26  33  65
98  80  82  89  91  48  30  32  39  41			23   5  82  89  91  48  30  32  39  66
79  81  88  95  97  29  31  38  45  47			 4   6  88  95  97  29  31  38  45  72
85  87  94  96  78  35  37  44  46  28			10  12  94  96  78  35  37  44  46  53
86  93 100  77  84  36  43  50  27  34			11  18 100  77  84  36  43  50  27  59
You will see our typical 5 x 5 magic
square in Q II (top left corner)

The finished M.S after the 2 switches
92  99   1   8  15  67  74  51  58  40
98  80   7  14  16  73  55  57  64  41
79   6  88  20  22  54  56  63  70  47
85  87  19  21   3  60  62  69  71  28
86  93  25   2   9  61  68  75  52  34
17  24  76  83  90  42  49  26  33  65
23   5  82  89  91  48  30  32  39  66
4  81  13  95  97  29  31  38  45  72
10  12  94  96  78  35  37  44  46  53
11  18 100  77  84  36  43  50  27  59
"""

nx, k = n/2, (n-2)/4
if n < 6 or nx % 2 == 0: return False	#only allow even squares 4n+2, n>1

#make an odd nx x nx magic square
A = make_odd_magic_square(nx)

#fill in each quadrant with an augmentation of A according to algorithm
I = A + sc_add(A, 3*nx*nx)

#create initial square by concatenating I and J - column sums are "magic"
B = [I[i]+J[i] for i in range(n)]

#swap upper and lower halves of specific columns to make row sums "magic"
for j in range(k) + range(n-1, n-k, -1):
for i in range(nx):
B[i][j], B[i+nx][j] =  B[i+nx][j], B[i][j]

#switch middle values for 2 columns to make diagonals (and square) magic
B[nx-k-1][k-1], B[nx+k][k-1] = B[nx+k][k-1], B[nx-k-1][k-1]
B[nx-k-1][k], B[nx+k][k] = B[nx+k][k], B[nx-k-1][k]

return B

def make_even_4n_magic_square(n):
"""
Source: http://mathworld.wolfram.com/MagicSquare.html
An elegant method for constructing magic squares of doubly even order n=4m
is to draw Xs through each 4 x 4 subsquare and fill all squares in sequence.
Then replace each entry A[i][j] on a crossed-off diagonal by (n^2+1) - A[i][j]
or, equivalently, reverse the order of the crossed-out entries.

Here is a mask of the crossed out matrix:
1 0 0 1 1 0 0 1
0 1 1 0 0 1 1 0
0 1 1 0 0 1 1 0
1 0 0 1 1 0 0 1
1 0 0 1 1 0 0 1
0 1 1 0 0 1 1 0
0 1 1 0 0 1 1 0
1 0 0 1 1 0 0 1
"""
if n < 4 or n%4: return False	#only allow even squares 4n, n>0

c, cms, A = 1, n*n + 1, [[0]*n for i in range(n)]
for i in range(n):
for j in range(n):
A[i][j] = cms-c if i%4 == j%4 or (i+j)%4 == (n-1)%4 else c
c += 1

return A

def make_even_panmagic_square(n):
"""
Description and algorithm from Wikipedia
https://en.wikipedia.org/wiki/Pandiagonal_magic_square
A pandiagonal magic square remains pandiagonally magic not only under
rotation or reflection, but also if a row or column is moved from one
side of the square to the opposite side. As such, an n x n pandiagonal
magic square can be regarded as having 8n^2 orientations.
No (4n+2)x(4n+2) panmagic square exists.
"""
if n < 4 or n%4: return False	#only allow 4n even squares
nx = n/2
x = range(1, n+1)
"""
Integrating a few steps into one, compose a matrix as follows:
(example for n=2 (8x8))
[1,2,3,4,8,7,6,5]
[8,7,6,5,1,2,3,4]
...
...
repeat pattern until full, then rotate 90 degrees
"""
A = (x[:nx]+x[:nx-1:-1], x[:nx-1:-1]+x[:nx]) * nx
B = rotate_ccw(A)

"""
Build the final square by multiplying the B square by n, adding
the A square and subtract n from each cell.
Example: A + n x B - n
"""
return [[a + n*b - n for a, b in zip(r1, r2)] for r1, r2 in zip(A, B)]

def make_odd_panmagic_square(n):
#order 6n +/- 1: 5, 7, 11, 13, 17, 19, 23, ...
if n > 4 and (n%6 == 1 or n%6 == 5):
A = [[(j*2 + i) % n + 1 for j in range(n)] for i in range(n)]

#order 6n + 3: 9, 15, 21, 27, ...
elif n > 8 and n%6 == 3:
q, dir = n / 3, 1
B = [[1,2,3],
[5,6,4],
[9,7,8]]

for i in range(10, n, 3):
B.append(range(i, i+3)[::dir])
dir = -dir

A = [[0]*n for i in range(n)]
for i in range(q):
for j in range(3):
t = B[i % q][j % 3]
for k in range(n):
A[(i+k) % n][(j + 3*k) % n] = t
else:
return False

T = transpose(A)
return [[a + n*b - n for a, b in zip(r1, r2)] for r1, r2 in zip(A, T)]

def check_magic_square(matrix, seq=True, mc=0):
"""
seq: Boolean used to restrict the square to having distinct numbers from 1..n^2
mc: Integer to specify a magic constant
A (normal) magic square is a square array(1) of numbers consisting of the distinct
positive integers 1, 2, ..., n^2(2) arranged such that the sum of the
n numbers in any horizontal(3), vertical(4), or main diagonal(5) line is always the
same number (Kraitchik 1942, p. 142; Andrews 1960, p. 1; Gardner 1961, p. 130;
Madachy 1979, p. 84; Benson and Jacoby 1981, p. 3; Ball and Coxeter 1987, p. 193),
known as the magic constant(6)

A square that fails to be magic only because one or both of the main diagonal sums
do not equal the magic constant is called a semimagic square(5)

It is an unsolved problem to determine the number of magic squares of an arbitrary
order, but the number of distinct magic squares (excluding those obtained
by rotation and reflection) of order n=1, 2, 3, ... are 1, 0, 1, 880, 275305224, ...
(Sloane's A006052; Madachy 1979, p. 87).(6)
"""
global err_message

#(6) check for invalid size and shape (non-square)
n = square_size(matrix)
if (n < 1 or n == 2):
err_message = "The array must be square and can't be <1 or 2."
return False

#(1) quick check for a square array
if (not isinstance(matrix[0], (list, tuple)) or
any(len(row) != n for row in matrix)):
err_message = "Must be a square, i.e. the same number of rows and columns."
return False

#(6) calculate the magic constant and working variables
if not mc: mc = n * (n*n + 1) / 2
d1, d2, r, c, a = 0, 0, [0]*n, [0]*n, set()
for i in range(n):
d1, d2 = d1+matrix[i][i], d2+matrix[i][n-i-1]
for j in range(n):
r[i] += matrix[i][j]
c[j] += matrix[i][j]

#(2) check distinct positive integers 1, 2, ..., n^2
if seq == True and not all(x in a for x in range(1, n*n+1)):
err_message = "The numbers 1 through " + str(n*n) + \
" must appear once and only once in the array."
return False

#(3) check that all horizontal lines equal the magic constant
row_sum = set(r)
if (len(row_sum) != 1 or mc not in row_sum):
err_message = "All the rows didn't add to the magic constant, " + str(mc)
return False

#(4) check that all vertical lines equal the magic constant
col_sum = set(c)
if (len(col_sum) != 1 or mc not in col_sum):
err_message = "All the columns didn't add to the magic constant, " + str(mc)
return False

#(5) check that both diagonal lines equal the magic constant
if (d1 != mc or d2 != mc):
err_message = "One or both diagonals didn't add to the magic constant, " \
+ str(mc) + "\nThis is a semimagic square."
return False

return True

def check_panmagic_square(A):
"""
If all diagonals (including those obtained by wrapping around) of a magic
square sum to the magic constant, the square is said to be a panmagic
square (also called a diabolic square or pandiagonal square). There are no
panmagic squares of order 3 or 4n+2.
"""
n = square_size(A)
if n < 1 or n == 3 or (n-2)%4 == 0: return False

if check_magic_square(A):	#make sure A is magic before panmagic test
mc = n * (n*n + 1) / 2	#magic constant

for i in range(n):
s1 = sum(A[(i-j) % n][j] for j in range(n))
s2 = sum(A[(i+j) % n][j] for j in range(n))
if s1 != mc or s2 != mc: return False
return True

return False

def check_associative_square(A):
"""
Reference: https://en.wikipedia.org/wiki/Associative_magic_square
An associative magic square is a magic square for which every pair
of numbers symmetrically opposite to the center sum up to the same
value.  There are no associative magic squares of singly-even order.
"""
n = square_size(A)
if n < 1 or n%4 != 0 and n%2 == 0: return False

if check_magic_square(A):	#make sure A is magic before associative test
amc = n*n + 1	#associative magic constant
odd_square = n % 2
if all(A[i][j] + A[n-i-1][n-j-1] == amc
for i in range(n) for j in range(n/2 + odd_square)):
return True

return False

def check_bimagic_square(A):
"""
If replacing each number in A with its square produces another magic
square, then the square is said to be a bimagic square (or doubly
magic square).
"""

if check_magic_square(A):	#make sure A is magic before bimagic test
n = len(A)
np = n*n
mc = n * (np + 1) * (2*np + 1) / 6	#magic constant

if check_magic_square([[x*x for x in row] for row in A], False, mc):
return True

return False

"""
Basic matrix operations in native Python.  Most of the time we simulate
these operations by accessing the matrix in a different order.
"""
def square_size(A):
try:
n = len(A)
except:
return 0

if (not isinstance(A[0], (list, tuple)) or
any(len(row) != n for row in A)):
return 0
return n

def transpose(A):
"""
reflect A over its main diagonal (which runs from top-left to bottom-right)
Example:
1, 2, 3		1, 4, 7
4, 5, 6  => 2, 5, 8
7, 8, 9		3, 6, 9

"""
return [list(a) for a in zip(*A)]

def rotate_cw(A):
"""
rotatet A clockwise (right-turn, -90 degrees)
Example:
1, 2, 3		7, 4, 1
4, 5, 6  => 8, 5, 2
7, 8, 9		9, 6, 3

"""
return [list(a) for a in zip(*A[::-1])]

def rotate_ccw(A):
"""
rotate A counterclockwise (left-turn, 90 degrees)
Example:
1, 2, 3		3, 6, 9
4, 5, 6  => 2, 5, 8
7, 8, 9		1, 4, 7

"""
return [list(a) for a in zip(*A)[::-1]]

def flip(matrix):
"""
flip A over horizontally
Example:
1, 2, 3		7, 8, 9
4, 5, 6  => 4, 5, 6
7, 8, 9		1, 2, 3

"""
return matrix[::-1]

"""
Scalar add n to matrix A
"""
return [[x+n for x in row] for row in A]

def print_matrix(matrix):
"""
Source: David Robinson from
http://stackoverflow.com/questions/8747500/pythons-pretty-printing-of-matrix
"""
max_lens = [max([len(str(r[i])) for r in matrix])
for i in range(len(matrix[0]))]

print "\n".join(["".join([string.rjust(str(e), l + 2)
for e, l in zip(r, max_lens)]) for r in matrix])

def demo(A):
if check_magic_square(A):
print_matrix(A)
print "this is a magic square"
else:
print err_message
return False
if check_panmagic_square(A): print "**panmagic"
if check_associative_square(A):	print "**associative"
if check_bimagic_square(A):	print "**bimagic"
print

def swap_col(A, n, m):
l = len(A)
for i in range(l):
A[i][n], A[i][m] = A[i][m], A[i][n]
return A

def swap_row(A, n, m):
l = len(A)
for i in range(l):
A[n][i], A[m][i] = A[m][i], A[n][i]
return A

"""
n = 4*1
print "4n (doubly even) pandiagonal or panmagic square, n=",n/4
M = make_even_panmagic_square(n)
for x in range(1):
demo(M)
N = rotate_cw(M)
demo(rotate_cw(M))
demo(rotate_cw(rotate_cw(M)))
demo(rotate_ccw(M))
demo(flip(M))
demo(flip(rotate_ccw(M)))
demo(flip(rotate_cw(M)))
demo(flip(rotate_cw(rotate_cw(M))))
M=swap_row(N, 0,1)
M=swap_row(M, 2,3)

"""

n = 7
print "Odd magic square, n=",n
M = make_magic_square(n)
demo(M)

print "Odd panmagic square, 6n+1, 6n+5"
n = 7
M=make_panmagic_square(n)
demo(M)

n = 11
M=make_panmagic_square(n)
demo(M)

print "Odd panmagic square, 6n+3"
n = 9
M=make_panmagic_square(n)
demo(M)

n = 4*2
print "4n (doubly even),  magic square, n=",n/4
M = make_magic_square(n)
demo(M)

n = 4*2
print "4n (doubly even) panmagic square, n=",n/4
M = make_panmagic_square(n)
demo(M)

print "8 x 8 bimagic square test check"
T =[[56, 34,  8, 57, 18, 47,  9, 31],
[33, 20, 54, 48,  7, 29, 59, 10],
[26, 43, 13, 23, 64, 38,  4, 49],
[19,  5, 35, 30, 53, 12, 46, 60],
[15, 25, 63,  2, 41, 24, 50, 40],
[ 6, 55, 17, 11, 36, 58, 32, 45],
[61, 16, 42, 52, 27,  1, 39, 22],
[44, 62, 28, 37, 14, 51, 21,  3]]
demo(T)

n=4*2 + 2
print "Even, 4n + 2, magic square, n=",(n-2)/4
M = make_magic_square(n)
demo(M)

print "Even, 4n, square associative test"
T = [[2,13,16,3],
[11,8,5,10],
[7,12,9,6],
[14,1,4,15]]
demo(T)

print "Odd square associative test"
T = [[25,22,5,10,3],
[7,12,11,17,18],
[2,6,13,20,24],
[8,9,15,14,19],
[23,16,21,4,1]]
demo(T)
```

## Discussion

### One Response to “Panmagic (Pandiagonal) and magic square generation, testing”

1. mind blowing! Thank you

Posted by Eduardo Antonio Fernandez Diaz | November 5, 2016, 8:05 PM