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

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.

Answer 2

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 3

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 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

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.