Magic Squares 4x4


Targets

Magic squares were made famous by German artist Albrecht Dürer. There are a lot of possibilities to arrange the numbers 1 to 16 in order to give the sum of 34 added horizontally and vertically and even in both diagonals. Mathematicians like to call it an 8-dimensional vector space because you can choose 8 values freely.

As said in final remark 2, for his copper engraving Dürer selected a very special one of all possible solutions. The sketch shown below will give you all the solutions.

The magic square of Dürer, highlighted in the yellow box, is shown enlarged in the picture below.

How it was done

normal sequence:
1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16
modified sequence:
1 2 3 4
5 6 7 8
9 11 13 14
10 12 15 16
In this approach a recursive way was implemented. Actually, if no tricks were used the Arduino needs some hours to find all solutions. The speed for finding them all can be improved considerably by eliminating unwanted selections and symmetrical solutions at an early stage.

Also, if possible you should give priority to calculating over brute-force. That is why cell 13 is calculated before iterating cell 10.

Mind you, the procedure Try is called 441,354,018 times.
(Of course, there is also an iterative algorithm.)

As arrays in C always start with index 0 all indexes have to be decreased by 1.


#define FILENAME "magic-squares-4x4, recursive"
// gives all 880 solutions
const byte cols = 4;
const byte max = cols * cols;
const byte sum = (1 + max) * max / 2 / cols;
byte p[max];
const byte seq[] = {0,  1,  2,  3, 
                    4,  5,  6,  7, 
                    8, 10, 12, 13, 
                    9, 11, 14, 15};
int num = 0;

void setup() {
  Serial.begin(115200); // set terminal baud rate accordingly 
  Serial.println(FILENAME);
  long t1 = millis();
  place(0, 0xFFFF);  // all 16 numbers are available 
  Serial.print(millis() - t1);
  Serial.println(" milliseconds. The End.");
}

void loop() {} // nothing more to do

void place(byte pos, int avail) {
  byte i;
  switch (pos) {
    case  0: case 1: case 2: case 4: case 5: case 6: case 8: case 10:  
             for (i = max; i >= 1; i--) { Try(pos, i, avail); } break;
    case  3: i = rem3(0,  1,  2); 
                    if (p[0] > i) Try(pos, i, avail); break;   
             // row 1, mirrored at vertical axis
    case  7: i = rem3(4,  5,  6); Try(pos, i, avail); break;   // row 2
    case  9: i = rem3(0,  4,  8); 
    if ((p[0] > i) && (p[3] > i)) Try(pos, i, avail); break; // column 1
    case 11: i = rem3(1,  5, 10); Try(pos, i, avail); break;   // column 2
    case 12: i = rem3(5,  6, 10); Try(pos, i, avail); break;   // central block 
    case 13: i = rem3(8, 10, 12); Try(pos, i, avail); break;   // row 3
    case 14: i = rem3(2,  6, 12); Try(pos, i, avail); break;   // column 3
    case 15: i = rem3(3,  7, 13); 
                    if (p[0] > i) Try(pos, i, avail); break;   
             // column 4, mirrored at main diagonal
    case 16: if (rem3(0, 5, 12) == p[15]) printAll(); // main diagonal 
  }
} 
Serial terminal window
  
...
365
 +----|----|----|----+
 | 16 |  2 |  9 |  7 |
 +----|----|----|----+
 |  5 | 11 |  4 | 14 |
 +----|----|----|----+
 |  3 | 13 |  6 | 12 |
 +----|----|----|----+
 | 10 |  8 | 15 |  1 |
 +----|----|----|----+

366
 +----|----|----|----+
 | 16 |  3 |  2 | 13 |
 +----|----|----|----+
 |  5 | 10 | 11 |  8 |
 +----|----|----|----+
 |  9 |  6 |  7 | 12 |
 +----|----|----|----+
 |  4 | 15 | 14 |  1 |
 +----|----|----|----+

367
 +----|----|----|----+
 | 16 |  3 |  2 | 13 |
 +----|----|----|----+
 |  9 |  6 |  7 | 12 |
 +----|----|----|----+
 |  5 | 10 | 11 |  8 |
 +----|----|----|----+
 |  4 | 15 | 14 |  1 |
 +----|----|----|----+
...
void Try(byte pos, byte value, int avail) {   if ((value < 1) || (value > max)) return; // bad number   int mask = 1 << value-1;   if (avail & mask) { // test if still available     p[pos] = value;     place(pos + 1, avail & ~mask);   // remove from the available list   } } byte rem3(byte a, byte b, byte c) { // return the remainder   return sum - p[a] - p[b] - p[c]; } void printAll() {   String s = "\n +----+----+----+----+";   Serial.print(++num);   Serial.println(s);   for (byte i = 0; i < max; i++) {     printIt(p[seq[i]]);     if ((i+1) % cols == 0) Serial.println(" |"+s);   }   Serial.println(); } void printIt(byte x) {   if (x < 10) Serial.print(" | "); else Serial.print(" | ");   Serial.print(x); }

Final remark 1:

The Serial terminal output was chosen for simplicity. You might prefer to display the solutions on a 4x20 char LCD or a TFT. You also can show the results using 16 LED displays each showing 1½ digits; controlling of LEDs will be much faster (but you will need a lot of external hardware).

Have a look at the video showing a short sequence of the recursive backtracking algorithm on a 1.8" TFT (actually, most of the CPU time is wasted for displaying the data):


Final remark 2:

Note that Dürer choose a solution that showed the year of the completion of his work in the middle of the bottom cells. He also benefited from the fact that the numbers "4" and "1" also found in the bottom row can be interpreted as his initials "D" and "A" (in the olden days, the ASCII code had not been invented yet).

Final remark 3:

If the time the Arduino takes to calculate it all appears too long for you, just take an Arduino-Due or an Intel Galileo which will perform the algorithm much faster.

Final remark 4:

If you ain't got an Arduino just go to H.B. Meyer: Math'-pages and let your browser do the job. (373 lines of Javascript plus 63 lines HTML will be working for you.)




contact: nji(at)gmx.de