Problem

Given an even integer n where n % 4 == 2 (i.e. n = 4k + 2), construct a normal magic square of order n containing the integers 1..n^2 exactly once so that every row, every column and the two main diagonals sum to the magic constant M = n*(n^2+1)/2.

Singly-even orders (6, 10, 14, …) are traditionally constructed by splitting the grid into four (n/2 x n/2) odd-order magic subsquares and then swapping columns and a few cells. The Strachey (4-subsquare / composite) method below is simple to implement and widely used.

Examples

Example 1

1
2
3
4
5
6
7
8
9

Input: n = 6
Output: 
35  1  6 26 19 24
3  32  7 21 23 25
31  9  2 22 27 20
8  28 33 17 10 15
30  5 34 12 14 16
4  36 29 13 18 11
Explanation: One valid output
Similar Problems

Solution

Method 1 - 4-subsquare method

Intuition

  • Let n be the order where n % 4 == 2 and set m = n/2 (so m is odd). Build four m x m odd-order magic squares using the Siamese method and place them as four quadrants:

    | A | C | | D | B |

  • Use offsets so the four subsquares cover consecutive blocks of numbers:

    • A: 1 .. m^2
    • B: m^2+1 .. 2*m^2
    • C: 2*m^2+1 .. 3*m^2
    • D: 3*m^2+1 .. 4*m^2 (which equals n^2)

  • Compute k = (n - 2) / 4. The high-level swaps are:

    1. Swap the leftmost k columns of A with the corresponding columns of D.
    2. Swap the rightmost k-1 columns of C with the corresponding columns of B.
    3. Swap two special central cells between A and D: the middle cell in the leftmost column and the central cell of the subsquare. These correct the center coupling between quadrants.

The procedure yields a valid magic square for every n = 4k + 2.

Approach

  • Build four m x m odd-order magic squares using the Siamese method and place them as quadrants:

    | A | C | | D | B |

  • Use offsets so the quadrants cover consecutive ranges: A = 1..m^2, B = m^2+1..2*m^2, C = 2*m^2+1..3*m^2, D = 3*m^2+1..4*m^2.

  • Perform these operations in sequence:

    1. Swap the leftmost k = (n-2)/4 columns between A and D.
    2. Swap the rightmost k-1 columns between C and B.
    3. Swap two special central cells between A and D (use mid = m/2).

These steps preserve the row/column sums and yield a normal magic square of order n = 4k + 2.

Pseudocode

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

function generate(n):
  if n is odd: return siamese(n)
  if n % 4 == 0: return generate_doubly_even(n)

  m = n/2
  base = siamese(m)          # m x m odd magic square
  ans = zero matrix n x n

  # place subsquares with offsets (use m2 = m*m)
  for i in 0..m-1:
    for j in 0..m-1:
      ans[i][j]       = base[i][j] + 0        # A
      ans[i][j + m]   = base[i][j] + 2*m*m    # C
      ans[i + m][j]   = base[i][j] + 3*m*m    # D
      ans[i + m][j+m] = base[i][j] + 1*m*m    # B

  k = (n - 2) / 4
  swap leftmost k columns between A and D
  swap rightmost (k-1) columns between C and B
  mid = m/2
  swap two special cells between A and D using mid

  return ans

Notes:

  • k = (n - 2) / 4 controls the number of full-column swaps.

  • mid = m/2 is used when swapping the two special central cells.

Dry run

Here n = 6, m = 3, m^2 = 9, k = (6-2)/4 = 1. We use the standard 3×3 Lo-Shu (Siamese) pattern for each subsquare and add offsets of 0, 9, 18, 27 respectively.

Lo-Shu (base m = 3):

1
2
3

8 1 6
3 5 7
4 9 2

Make the subsquares by adding offsets:

1
2
3
4

A (1..9)  = base + 0
B (10..18)= base + 9
C (19..27)= base + 18
D (28..36)= base + 27

Place them as quadrants A | C over D | B to compose the 6×6 before swaps:

1
2
3
4
5
6
7

 8  1  6 | 26 19 24
 3  5  7 | 21 23 25
 4  9  2 | 22 27 20
 ---------+--------
35 28 33 | 17 10 15
30 32 34 | 12 14 16
31 36 29 | 13 18 11

Step 1 — swap leftmost k = 1 column of A and D (column index 0):

Before swap (A and D shown):

1
2
3
4
5
6
7
8
9

A:
8 1 6
3 5 7
4 9 2

D:
26 19 24
21 23 25
22 27 20

After swapping column 0 of A with column 0 of D:

1
2
3
4
5
6
7
8
9

A becomes:
35 1 6
30 5 7
31 9 2

D becomes:
8 28 33
3 32 34
4 36 29

Step 2 — swap the two special central cells between A and D (using mid = 1):

Swap A[mid][0] <-> D[mid][0] and A[mid][mid] <-> D[mid][mid].

Final 6×6 after all swaps:

1
2
3
4
5
6

35  1  6 26 19 24
3  32  7 21 23 25
31  9  2 22 27 20
8  28 33 17 10 15
30  5 34 12 14 16
4  36 29 13 18 11

This matches the classical 6×6 Strachey example and every row, column and main diagonal sums to 111.

Code

1
2


##### C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

using namespace std;

class Solution {
 public:
  vector<vector<int>> generate(int n) {
    if (n % 2 == 1) return generateOdd(n);
    if (n % 4 == 0) return generateDoublyEven(n);

    int m = n / 2;
    int m2 = m * m;
    // build base m x m odd magic square (Siamese method)
    auto base = siamese(m);

    vector<vector<int>> ans(n, vector<int>(n));
    // offsets: A:0, B:m2, C:2*m2, D:3*m2
    for (int i = 0; i < m; ++i)
      for (int j = 0; j < m; ++j) {
        ans[i][j] = base[i][j] + 0;            // A
        ans[i][j + m] = base[i][j] + 2 * m2;  // C
        ans[i + m][j] = base[i][j] + 3 * m2;  // D
        ans[i + m][j + m] = base[i][j] + m2;  // B
      }

    int k = (n - 2) / 4;
    // swap leftmost k columns of A and D
    for (int col = 0; col < k; ++col)
      for (int row = 0; row < m; ++row)
        swap(ans[row][col], ans[row + m][col]);

    // swap rightmost (k-1) columns of C and B
    for (int col = m - (k - 1); col < m; ++col)
      for (int row = 0; row < m; ++row)
        swap(ans[row][col + m], ans[row + m][col + m]);

    // final special swaps (middle of leftmost column and central cell)
    int mid = m / 2;
    swap(ans[mid][0], ans[mid + m][0]);
    swap(ans[mid][mid], ans[mid + m][mid]);

    return ans;
  }

 private:
  vector<vector<int>> siamese(int m) {
    vector<vector<int>> a(m, vector<int>(m, 0));
    int num = 1;
    int i = 0, j = m / 2;
    while (num <= m * m) {
      a[i][j] = num++;
      int ni = (i - 1 + m) % m;
      int nj = (j + 1) % m;
      if (a[ni][nj]) i = (i + 1) % m;
      else { i = ni; j = nj; }
    }
    return a;
  }

  // Minimal helpers for parity cases; straightforward implementations.
  vector<vector<int>> generateOdd(int n) { return siamese(n); }
  vector<vector<int>> generateDoublyEven(int n) {
    vector<vector<int>> ans(n, vector<int>(n));
    int cnt = 1;
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        ans[i][j] = cnt++;
    int total = n * n + 1;
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j) {
        int r = i % 4, c = j % 4;
        if (r == c || r + c == 3) ans[i][j] = total - ans[i][j];
      }
    return ans;
  }
};
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

class Solution {
    public int[][] generate(int n) {
        if (n % 2 == 1) return generateOdd(n);
        if (n % 4 == 0) return generateDoublyEven(n);

        int m = n / 2;
        int m2 = m * m;
        int[][] base = siamese(m);
        int[][] ans = new int[n][n];

        // place A, C, D, B with offsets 0, 2*m2, 3*m2, m2 respectively
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < m; ++j) {
                ans[i][j] = base[i][j];                // A
                ans[i][j + m] = base[i][j] + 2 * m2;   // C
                ans[i + m][j] = base[i][j] + 3 * m2;   // D
                ans[i + m][j + m] = base[i][j] + m2;   // B
            }
        }

        int k = (n - 2) / 4;
        // swap leftmost k columns of A and D
        for (int col = 0; col < k; ++col)
            for (int row = 0; row < m; ++row)
                swap(ans, row, col, row + m, col);

        // swap rightmost (k-1) columns of C and B
        for (int col = m - (k - 1); col < m; ++col)
            for (int row = 0; row < m; ++row)
                swap(ans, row, col + m, row + m, col + m);

        int mid = m / 2;
        swap(ans, mid, 0, mid + m, 0);
        swap(ans, mid, mid, mid + m, mid);

        return ans;
    }

    private void swap(int[][] a, int r1, int c1, int r2, int c2) {
        int t = a[r1][c1]; a[r1][c1] = a[r2][c2]; a[r2][c2] = t;
    }

    private int[][] siamese(int m) {
        int[][] a = new int[m][m];
        int num = 1;
        int i = 0, j = m / 2;
        while (num <= m * m) {
            a[i][j] = num++;
            int ni = (i - 1 + m) % m;
            int nj = (j + 1) % m;
            if (a[ni][nj] != 0) i = (i + 1) % m;
            else { i = ni; j = nj; }
        }
        return a;
    }

    private int[][] generateOdd(int n) { return siamese(n); }

    private int[][] generateDoublyEven(int n) {
        int[][] ans = new int[n][n];
        int cnt = 1;
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                ans[i][j] = cnt++;
        int total = n * n + 1;
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                int r = i % 4, c = j % 4;
                if (r == c || r + c == 3) ans[i][j] = total - ans[i][j];
            }
        return ans;
    }
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

class Solution:
    def generate(self, n: int) -> list[list[int]]:
        if n % 2 == 1:
            return self.siamese(n)
        if n % 4 == 0:
            return self.generate_doubly_even(n)

        m = n // 2
        m2 = m * m
        base = self.siamese(m)
        ans = [[0] * n for _ in range(n)]

        # place A, C, D, B with offsets 0, 2*m2, 3*m2, m2
        for i in range(m):
            for j in range(m):
                ans[i][j] = base[i][j]
                ans[i][j + m] = base[i][j] + 2 * m2
                ans[i + m][j] = base[i][j] + 3 * m2
                ans[i + m][j + m] = base[i][j] + m2

        k = (n - 2) // 4
        # swap leftmost k columns of A and D
        for col in range(k):
            for row in range(m):
                ans[row][col], ans[row + m][col] = ans[row + m][col], ans[row][col]

        # swap rightmost (k-1) columns of C and B
        for col in range(m - (k - 1), m):
            for row in range(m):
                ans[row][col + m], ans[row + m][col + m] = ans[row + m][col + m], ans[row][col + m]

        mid = m // 2
        ans[mid][0], ans[mid + m][0] = ans[mid + m][0], ans[mid][0]
        ans[mid][mid], ans[mid + m][mid] = ans[mid + m][mid], ans[mid][mid]

        return ans

    def siamese(self, m: int) -> list[list[int]]:
        a = [[0] * m for _ in range(m)]
        num = 1
        i, j = 0, m // 2
        while num <= m * m:
            a[i][j] = num
            num += 1
            ni, nj = (i - 1) % m, (j + 1) % m
            if a[ni][nj]:
                i = (i + 1) % m
            else:
                i, j = ni, nj
        return a

    def generate_doubly_even(self, n: int) -> list[list[int]]:
        ans = [[0] * n for _ in range(n)]
        cnt = 1
        for i in range(n):
            for j in range(n):
                ans[i][j] = cnt
                cnt += 1
        total = n * n + 1
        for i in range(n):
            for j in range(n):
                r, c = i % 4, j % 4
                if r == c or r + c == 3:
                    ans[i][j] = total - ans[i][j]
        return ans

Complexity

  • ⏰ Time complexity: O(n^2) — building the subsquares and performing swaps visits O(n^2)
  • 🧺 Space complexity: O(n^2) for the output matrix.