What is the optimal way of generating all possible 3X3 magic squares?

Question

Magic square: sum of any row, column, or diagonal of length is always equal to the same number. All 9 numbers are distinct positive integers.

I am doing it this way in JavaScript, but what is the optimal way of generating all of them?

function getMagicSquare() {

let myArray = [
    [4, 9, 2],
    [3, 5, 7],
    [8, 1, 5]
];

for (let index1 = 1; index1 < 10; index1++) {
    for (let index2 = 1; index2 < 10; index2++) {
        for (let index3 = 1; index3 < 10; index3++) {
            for (let index4 = 1; index4 < 10; index4++) {
                for (let index5 = 1; index5 < 10; index5++) {
                    for (let index6 = 1; index6 < 10; index6++) {
                        for (let index7 = 1; index7 < 10; index7++) {
                            for (let index8 = 1; index8 < 10; index8++) {
                                for (let index9 = 1; index9 < 10; index9++)
                                // if numbers are not distinct for each loop, I can break the loop and make it a bit faster
                                {
                                    const mySet = new Set();
                                    mySet.add(index1).add(index2).add(index3).add(index4).add(index5).add(index6).add(index7).add(index8).add(index9)
                                    if ((mySet.size === 9))
                                        if (
                                            (index1 + index2 + index3 === index4 + index5 + index6) &&
                                            (index4 + index5 + index6 === index7 + index8 + index9) &&
                                            (index7 + index8 + index9 === index1 + index4 + index7) &&
                                            (index1 + index4 + index7 === index2 + index5 + index8) &&
                                            (index2 + index5 + index8 === index3 + index6 + index9) &&
                                            (index3 + index6 + index9 === index1 + index5 + index9) &&
                                            (index1 + index5 + index9 === index3 + index5 + index7)
                                        ) {
                                            myArray[0][0] = index1;
                                            myArray[0][1] = index2;
                                            myArray[0][2] = index3;
                                            myArray[1][0] = index4;
                                            myArray[1][1] = index5;
                                            myArray[1][2] = index6;
                                            myArray[2][0] = index7;
                                            myArray[2][1] = index8;
                                            myArray[2][2] = index9;

                                            console.log(myArray);

                                        }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

}

Second question: What if I want to generate NxN magic squares?

Answer 1

First, you have to buy the fact that the center value, i.e., yourArray[1][1], must be 5 for all magic squares, otherwise turn to this answer for the reasoning.

With the belief that the center value, let's denote it as c, must be 5, it's obvious albeit may need some pondering that, we can construct an initial matrix using the allocation of numbers below, where all elements must be larger than 0, also a and b cannot be 2 times the other, otherwise you get identical numbers in the matrix, so a=1 and b=3 is the only valid combinations.

c+a      c-a-b      c+b
c-a+b      c        c+a-b
c-b      c+a+b      c-a

If you find it's hard to figure out the correct layout, you need to give some pondering on the concept of Moment of inertia and Center of mass. Shortly speaking, you have to keep the element with heavier weight, i.e., c+a+b, c-a-b, etc, as close to center as possible, and keep the element with lighter weight away from the center.

Once we have the initial matrix to start with, remaining steps are much easier, given the property of magic square, rotating the matrix by multiple of 90 degree, mirroring the matrix, rotating the mirror by multiple of 90 degree, as well as the original matrix, give us all valid magic squares:

m = [[4, 9, 2],
     [3, 5, 7],
     [8, 1, 6]]

def rotateBy45Degree( A ):
    # col3 -> row1
    row1 = [A[i][2] for i in range( 3 )]
    # col2 -> row2
    row2 = [A[i][1] for i in range( 3 )]
    # col1 -> row3
    row3 = [A[i][0] for i in range( 3 )]
    return [row1, row2, row3]

print( m )

next = rotateBy45Degree( m )
while next != m:
    print( next )
    next = rotateBy45Degree( next )

# mirror
m[0][0], m[1][0], m[2][0], m[0][2], m[1][2], m[2][2] = m[0][2], m[1][2], m[2][2], m[0][0], m[1][0], m[2][0]
print( m )

next = rotateBy45Degree( m )
while next != m:
    print( next )
    next = rotateBy45Degree( next )

Answer 2

Here is my solution. It's very fast but has a significant memory footprint. This implementation would not be good for squares of order 4+ as then you would be attempting to load into memory 16 factorial permutations.

Let me know what you all think.

import numpy as np
import itertools as it

a = it.permutations((range(1, 10)))
b = np.fromiter(it.chain(*a), dtype=np.uint8).reshape(-1,3,3)

ma = np.array([15,15,15])

rs = np.sum(b.reshape(-1,3), axis=1).reshape(-1,3) #row sums
cs = np.sum(b, axis=1) #col sums
ds = np.trace(b, axis1=2, axis2=1) #diag sums
dr = np.trace(np.flip(b,axis=1), axis1=2) #diag flipped sums

i = np.all(rs == ma, axis=1) & np.all(cs == ma, axis=1) & (ds == 15) & (dr == 15)
r = b[i]

print(r)

Answer 3

1. There is only one 3×3 magic square (up to rotations and mirroring). The optimal way of generating it is to hardcode it (or the 8 rotations and mirror images of it) in the program.
2. Enumeration of N×N magic squares is an open problem. It is not even known how many 6×6 magic squares exist (it is estimated that the number is about 1.77e19) wikipedia.

Answer 4

Not much optimised but somewhat cleaner and easy to understand approach.

bool is_magic(int n, vector<int> vec, int m_sum) {
    vector<vector<int>> values(3, vector<int>(3, -1));

    for(int i=0; i<9; i++) values[i/3][i%3] = vec[i];

    for (int i=0; i<n; i++) {
        int r_sum = 0;
        int c_sum = 0;
        int ld_sum = 0;
        int rd_sum = 0;
        
        for (int j=0; j<n; j++) {
            r_sum+=values[i][j];
            c_sum+=values[j][i];
            ld_sum+=values[j][j];
            rd_sum+=values[j][n-1-j];
        }
        if (r_sum!=m_sum) return false;
        if (c_sum!=m_sum) return false;
        if (ld_sum!=m_sum) return false;
        if (rd_sum!=m_sum) return false;
    }
    return true;
}

void magicSquare(int n) {
    vector<int> values = {1,2,3,4,5,6,7,8,9};
    int m_sum = accumulate(values.begin(), values.end(), 0)/3;
    
    vector<vector<int> > combs;
    do {
        if(is_magic(n, values, m_sum)) combs.push_back(values);
    }  while(next_permutation(values.begin(), values.end())); 
}

output=>

2 7 6 
9 5 1 
4 3 8 

2 9 4 
7 5 3 
6 1 8 

4 3 8 
9 5 1 
2 7 6 

4 9 2 
3 5 7 
8 1 6 

6 1 8 
7 5 3 
2 9 4 

6 7 2 
1 5 9 
8 3 4 

8 1 6 
3 5 7 
4 9 2 

8 3 4 
1 5 9 
6 7 2 

[Done] exited with code=0 in 1.65 seconds

Answer 5

Below is the Javascript implementation for generating all 8 3x3 magic squares.

Algorithm used:

1. Generate one 3x3 magic square (geeksforgeeks article).
2. Derive the remaining magic squares by reflections and rotations (based on Presh Talwalkar's blog). There is another method where you can generate the first 4 sets of 3x3 magic square and then derive the remaining 4 by subtracting 10 from each element.

function generateMagic3x3(n) {
  let i, j;
  i = Math.floor(n / 2);
  j = n - 1;
  let baseMatrix = [
    [],
    [],
    []
  ];
  baseMatrix[i][j] = 1;

  for (let k = 2; k <= n * n; k++) {
    i -= 1;
    j += 1;

    if (i < 0 && j === n) {
      i = 0;
      j = n - 2;
    } else if (i < 0) {
      i = n - 1;
    } else if (j === n) {
      j = 0;
    }

    if (typeof baseMatrix[i][j] === 'number') {
      i += 1;
      j -= 2;
    }

    baseMatrix[i][j] = k;
  }
  const baseMatrix2 = reflectDiag(baseMatrix);
  renderMatrix(baseMatrix)
  renderMatrix(reflectRows(baseMatrix));
  renderMatrix(reflectColumns(baseMatrix));
  renderMatrix(reflectColumns(reflectRows(baseMatrix)));
  renderMatrix(baseMatrix2);
  renderMatrix(reflectRows(baseMatrix2));
  renderMatrix(reflectColumns(baseMatrix2));
  renderMatrix(reflectColumns(reflectRows(baseMatrix2)));

};


function reflectColumns(matrix) {
  var newMatrix = matrix.map(function(arr) {
    return arr.slice();
  });
  for (let row = 0; row < matrix.length; row++) {
    newMatrix[row][0] = matrix[row][2];
    newMatrix[row][2] = matrix[row][0];
  }
  return newMatrix;
}

function reflectRows(matrix) {
  var newMatrix = matrix.map(function(arr) {
    return arr.slice();
  });
  for (let column = 0; column < matrix.length; column++) {
    newMatrix[0][column] = matrix[2][column];
    newMatrix[2][column] = matrix[0][column];
  }
  return newMatrix;
}

function reflectDiag(matrix) {
  var newMatrix = matrix.map(function(arr) {
    return arr.slice();
  });
  for (let row = 0; row < matrix.length; row++) {
    for (let column = 0; column < matrix.length; column++) {
      if (row !== column) {
        newMatrix[row][column] = matrix[column][row];
      }
    }
  }
  return newMatrix;
}

function renderMatrix(matrix) {
  const table = document.createElement('table');
  let resBox = document.getElementById('res')
  for (let row = 0; row < matrix.length; row++) {
    const tr = table.insertRow(row);
    for (let column = 0; column < matrix.length; column++) {
      const cell = tr.insertCell(column);
      cell.innerHTML = matrix[row][column];
    }
  }
  resBox.appendChild(table)
}


generateMagic3x3(3);

table {
  border-collapse: collapse;
  display:inline-block;
  margin: 10px;
}

td {
  border: 1px solid #000;
  text-align: center;
  width: 50px;
  height: 50px;
}

<div id='res'></div>

Answer 6

Here's a pretty naive implementation using state-space search with basic pruning to generate all possible magic squares of given dimension n: https://ideone.com/0aewnJ

from collections import defaultdict, deque
from copy import copy, deepcopy
import time

def magicSum(n):
    return int((n*n * (n*n + 1)) / 6)

def validate(sumDict, n):
    for k, v in sumDict.items():
        if v > magicSum(n):
            return False
    return True

def check(sumDict, n):
    for k, v in sumDict.items():
        if v != magicSum(n):
            return False
    return True

def isValid(m, n):
    rowSum = defaultdict(int)
    colSum = defaultdict(int)
    diagSum = defaultdict(int)

    isLeft = False

    for i in range(n):
        for j in range(n):
            if m[i][j] == 0: isLeft = True
            rowSum[i] += m[i][j]
            colSum[j] += m[i][j]
            if i == j: diagSum[0] += m[i][j]
            if i + j == n - 1: diagSum[n - 1] += m[i][j]

    if isLeft:
        return (validate(rowSum, n) and validate(colSum, n) and validate(diagSum, n))       
    return (check(rowSum, n) and check(colSum, n) and check(diagSum, n))        

def next(cur, m, n):
    possible = []
    for i in range(n): 
        for j in range(n):
            if m[i][j] == 0:
                nextM = deepcopy(m)
                nextM[i][j] = cur
                if isValid(nextM, n):
                    possible.append(nextM)
    return possible

def printM(m):
    for i in range(len(m)):
            print(m[i])
    print("\n")

def gen(n):
    startM = [[0 for x in range(n)] for y in range(n)]
    magic = []
    Q = deque([(1, startM)])
    while len(Q):
        state = Q.popleft()
        cur = state[0]
        m = state[1]
        if cur == n * n + 1:
            magic.append(m)
            printM(m)
            continue
        for w in next(cur, m, n):
            Q.append((cur + 1, w))
    return magic

start_time = time.time()
magic = gen(3)
elapsed_time = time.time() - start_time
print("Elapsed time: ", elapsed_time)

Output:

[6, 1, 8]
[7, 5, 3]
[2, 9, 4]


[8, 1, 6]
[3, 5, 7]
[4, 9, 2]


[6, 7, 2]
[1, 5, 9]
[8, 3, 4]


[8, 3, 4]
[1, 5, 9]
[6, 7, 2]


[2, 7, 6]
[9, 5, 1]
[4, 3, 8]


[4, 3, 8]
[9, 5, 1]
[2, 7, 6]


[2, 9, 4]
[7, 5, 3]
[6, 1, 8]


[4, 9, 2]
[3, 5, 7]
[8, 1, 6]


Elapsed time:  13.479725122451782

Though must I say that it performs a bit poorly in terms of run time than expected but I guess it's still good for a start. Will try to optimize it further in a while.