Met het onderstaande stukje code (C++) reken je magische vierkanten uit met dezelfde eigenschappen als het gevonden. Na een uurtje of 10 draaien heeft hij zijn hele zoekruimte doorlopen, vind een half miljoen oplossingen en maakt een file met oplossingen die ongeveer een half GB is. De code moet nog opgeschoond worden, maar werkt dus wel

.
De oplossosing met 'SeqNr: 42177' is die de scholieren hebben gevonden. In de telegraaf stond zaterdag een reactie van een van de teleurgestelde vinders, op wat negatieve publiciteit, dat je dit echt niet 'even' met de computer doet, helaas doe je dit dus wel.
P.S. de indentatie in de code is hier een beetje weg.
#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;
template<int size> class Square
{
private:
int _square[size][size];
bool _free[size*size];
int _size2;
int _quadSum;
int _rowColSum;
int _printedCount;
public:
Square()
{
_printedCount = 0;
_size2 = size * size;
_quadSum = (_size2+1)*2;
_rowColSum = (_size2+1)*(size/2);
for (int y = 0; y < size; y++)
{
for (int x = 0; x < size; x++)
{
_square[x][y] = 0;
}
}
for (int n = 0; n < _size2; n++)
{
_free[n] = true;
}
}
Square(const Square<size>& ref)
{
_printedCount = ref._printedCount;
_size2 = ref._size2;
_quadSum = ref._quadSum;
_rowColSum = ref._rowColSum;
memcpy(_square, ref._square, sizeof(_square));
memcpy(_free, ref._free, sizeof(_free));
}
Square<size>& operator =(const Square<size>& ref)
{
memcpy(_square, ref._square, sizeof(_square));
memcpy(_free, ref._free, sizeof(_free));
return *this;
}
bool Assign(int x, int y, int n)
{
bool result = false;
if (0 < n && n <= _size2 && _free[n-1])
{
bool free = true;
// Check if range is already in use
for (int i = 0; i < y && free; i++)
{
int div = (_square[0][i]-1)/size;
if ((div*size) < n && n <= ((div+1)*size))
{
free = false;
}
div = size-div;
if (((div-1)*size) < n && n <= (div*size))
{
free = false;
}
}
if (free)
{
_square[x][y] = n;
_free[n-1] = false;
int nMirror = (_size2+1)-n;
_square[x][(size-1)-y] = nMirror;
_free[nMirror-1] = false;
result = true;
}
}
return result;
}
void Print(ostream& os)
{
_printedCount++;
os << "SeqNr: " << _printedCount << endl << endl;
int xSum[size];
int xSum1[size];
int xSum2[size];
int quadSumError = 0;
int c8SumError = 0;
int c12SumError = 0;
int diagSumError = 0;
for (int y = 0; y < size; y++)
{
int sum = 0;
int sum1 = 0;
int sum2 = 0;
for (int x = 0; x < size; x++)
{
os << setw(4) << _square[x][y];
int quadSum = 0;
quadSum += _square[x][y];
quadSum += _square[(x+1)%size][y];
quadSum += _square[x][(y+1)%size];
quadSum += _square[(x+1)%size][(y+1)%size];
if (quadSum != _quadSum)
{
quadSumError++;
}
int c8Sum = 0;
c8Sum += _square[(x+1)%size][y];
c8Sum += _square[(x+2)%size][y];
c8Sum += _square[(x+3)%size][(y+1)%size];
c8Sum += _square[(x+3)%size][(y+2)%size];
c8Sum += _square[(x+2)%size][(y+3)%size];
c8Sum += _square[(x+1)%size][(y+3)%size];
c8Sum += _square[x][(y+2)%size];
c8Sum += _square[x][(y+1)%size];
if (c8Sum != _quadSum*2)
{
c8SumError++;
}
int c12Sum = 0;
c12Sum += _square[(x+2)%size][y];
c12Sum += _square[(x+3)%size][y];
c12Sum += _square[(x+4)%size][(y+1)%size];
c12Sum += _square[(x+5)%size][(y+2)%size];
c12Sum += _square[(x+5)%size][(y+3)%size];
c12Sum += _square[(x+4)%size][(y+4)%size];
c12Sum += _square[(x+3)%size][(y+5)%size];
c12Sum += _square[(x+2)%size][(y+5)%size];
c12Sum += _square[(x+1)%size][(y+4)%size];
c12Sum += _square[x][(y+3)%size];
c12Sum += _square[x][(y+2)%size];
c12Sum += _square[(x+1)%size][(y+1)%size];
if (c12Sum != _quadSum*3)
{
c12SumError++;
}
sum += _square[x][y];
if (x < (size/2))
{
sum1 += _square[x][y];
}
else
{
sum2 += _square[x][y];
}
if

{
xSum[x] += _square[x][y];
}
else
{
xSum[x] = _square[x][y];
}
if (y < (size/2))
{
if

{
xSum1[x] += _square[x][y];
}
else
{
xSum1[x] = _square[x][y];
}
}
else
{
if (y == (size/2))
{
xSum2[x] = _square[x][y];
}
else
{
xSum2[x] += _square[x][y];
}
}
}
int diagSum = 0;
if (y < size/2)
{
diagSum += _square[0][y];
diagSum += _square[1][y+1];
diagSum += _square[2][y+2];
diagSum += _square[3][y+3];
diagSum += _square[4][y+4];
diagSum += _square[5][y+5];
diagSum += _square[6][y+5];
diagSum += _square[7][y+4];
diagSum += _square[8][y+3];
diagSum += _square[9][y+2];
diagSum += _square[10][y+1];
diagSum += _square[11][y];
}
else
{
diagSum += _square[0][y];
diagSum += _square[1][y-1];
diagSum += _square[2][y-2];
diagSum += _square[3][y-3];
diagSum += _square[4][y-4];
diagSum += _square[5][y-5];
diagSum += _square[6][y-5];
diagSum += _square[7][y-4];
diagSum += _square[8][y-3];
diagSum += _square[9][y-2];
diagSum += _square[10][y-1];
diagSum += _square[11][y];
}
if (diagSum != _rowColSum)
{
diagSumError++;
}
os << setw(8 ) << sum
<< setw(4) << sum1
<< setw(4) << sum2
<< setw(4) << diagSum << endl << endl;
}
os << endl << endl;
for (int x = 0; x < size; x++)
{
os << setw(4) << xSum[x];
}
os << setw(8 ) << quadSumError
<< setw(4) << c8SumError
<< setw(4) << c12SumError
<< setw(4) << diagSumError << endl << endl;
for (x = 0; x < size; x++)
{
os << setw(4) << xSum1[x];
}
os << endl << endl;
for (x = 0; x < size; x++)
{
os << setw(4) << xSum2[x];
}
os << endl << endl << endl << endl;
if (diagSumError == 0)
{
cout << "SeqNr: " << _printedCount << endl;
}
}
bool Propagate(int x1, int y1)
{
if (x1 == 0 && y1%4 == 2 && _square[x1][y1+1] == 0)
{
int n = _quadSum - (_square[x1][y1-2] + _square[x1][y1-1] + _square[x1][y1]);
if (size < n && n <= (_size2-size) && Assign(x1, y1+1, n))
{
return Propagate(x1, y1+1);
}
else
{
return false;
}
}
else
{
// Propagate result of the assignment at x, y
// and check if no conflict arrises
int x0 = (x1-1+size)%size;
int y2 = (y1+1)%size;
if (_square[x1][y2] == 0 && _square[x0][y1] && _square[x0][y2])
{
int n = _quadSum - (_square[x0][y1] + _square[x0][y2] + _square[x1][y1]);
if (Assign(x1, y2, n))
{
return Propagate(x1, y2);
}
else
{
return false;
}
}
}
return true;
}
void InitLoop(int x, int y, int& start, int& stop, int& dir, int& xn, int& yn)
{
dir = 1;
if (x)
{
xn = x+1;
yn = 0;
}
else
{
int hs = size/2;
if (y%4 == 2 && _square[x][y+1] == 0)
{
xn = (y+2)/hs;
yn = (y+2)%hs;
}
else
{
xn = (y+1)/hs;
yn = (y+1)%hs;
}
}
if

{
start = size+1;
stop = (_size2-size)+1;
if (y == 2)
{
// TODO Some optimization is possible here,
// the start stop window can be made smaller
}
}
else
{
// See if it is a forced value
if (x%4 == 3)
{
start = _quadSum - (_square[x-3][0] + _square[x-2][0] + _square[x-1][0]);
if (size < start && start <= (_size2-size))
{
// Wrong range so stop loop
stop = start;
}
else
{
stop = start + 1;
}
}
else
{
// Check half row value, it must be near half row value
if ((x+1) == (size/2))
{
start = _rowColSum/2-1;
for (int i = 0; i < (size/2); i++)
{
start -= _square[i][0];
}
stop = start+3;
}
else
{
// We only take the extreme values at the top row
if (x%2)
{
stop = _size2-size;
start = _size2;
dir = -1;
}
else
{
start = 1;
stop = size+1;
}
}
}
}
}
void Fill(ostream& os, int x = 0, int y = 0)
{
if (x < size)
{
// Remember current status in order to restore
Square<size> status(*this);
int start, stop, dir, xn, yn;
InitLoop(x, y, start, stop, dir, xn, yn);
for (int n = start; n != stop; n += dir)
{
// Check if number is available
if (Assign(x, y, n))
{
// Propagate results and check if no conflict has pop up
if (Propagate(x, y))
{
// There are no conflicts, assign next cell
Fill(os, xn, yn);
}
// Restore status
*this = status;
}
}
}
else
{
bool valid = true;
for (int i = 0; i < size && valid; i++)
{
int diagSum = 0;
if (i < size/2)
{
diagSum += _square[0][i];
diagSum += _square[1][i+1];
diagSum += _square[2][i+2];
diagSum += _square[3][i+3];
diagSum += _square[4][i+4];
diagSum += _square[5][i+5];
diagSum += _square[6][i+5];
diagSum += _square[7][i+4];
diagSum += _square[8][i+3];
diagSum += _square[9][i+2];
diagSum += _square[10][i+1];
diagSum += _square[11][i];
}
else
{
diagSum += _square[0][i];
diagSum += _square[1][i-1];
diagSum += _square[2][i-2];
diagSum += _square[3][i-3];
diagSum += _square[4][i-4];
diagSum += _square[5][i-5];
diagSum += _square[6][i-5];
diagSum += _square[7][i-4];
diagSum += _square[8][i-3];
diagSum += _square[9][i-2];
diagSum += _square[10][i-1];
diagSum += _square[11][i];
}
if (diagSum != _rowColSum)
{
valid = false;
}
}
if (valid)
{
Print(os);
}
}
}
};
int _tmain(int argc, _TCHAR* argv[])
{
Square<12> square;
ofstream os("out.txt");
cout << "Print zero square to see if error counters on the sum of 2x2 squares and the circles of 8 & 12 and the diagonals are OK" << endl;
square.Print(os);
square.Fill(os);
return 0;
}