Home More Samples
```///////////////////////////////////////////////////////////////////////////////
//
// This program uses brute force to find all solutions to 5x5 magic squares.
// It is very important to determine number combinations where there are no
// solutions possible as early as possible. In order to do that, we specify
// the order in which the numbers are placed on the square. After each number
// is placed, we perform a position dependent sanity test. We do not want to
// test more nor less then it is neccessary.
// The placement order was chosen in such a way as to allow most efficient
// sanity tests. The tests are basically these:
//
//   1. Any three numbers in a line (row,column or diagonals) add up to less then 65
//   2. Any four numbers in a line (row, column or diagonals) add up to less than 65
//   3. Any five numbers in a line (row, column or diagonals) add up to exactly 65.
//
// In addition to these tests, there are a few tests to eliminate solutions that are
// symmetrical. Having a notation as follows:
//
//      a1 a2 a3 a4 a5
//      b1 b2 b3 b4 b5
//      c1 c2 c3 c4 c5
//      d1 d2 d3 d4 d5
//      e1 e2 e3 e4 e5
//
// we eliminate symmetrical solutions by requesting
//
//      a1 > a5 and a1 > e1 and a1 > e5 and a5 > e1
//
// Other requirements (see http://www.hbmeyer.de/backtrack/mq5/mag5the.htm)
// which can speed up calculations are possible, such as
//
//      26 <= a1 + a5 + e1 + e5 <= 78
//      52 <= a1 + a5 + e1 + e5 + 2*c3 <= 104
//
///////////////////////////////////////////////////////////////////////////////
//
// To generate all solutions, run query:
//
//      all MagicSquares5x5All()
//
// You may want to generate only a subset of all magic squares, for example
// to generate only 100000 solutions, use the following query:
//
//      all MagicSquares5x5All() & RtlTrimSolutions(100000)
//
///////////////////////////////////////////////////////////////////////////////
//
// Using Core2 @ 2.4GHz generates upon average about 800 solutions per second.
// At this rate it still takes about 95 hours to generate all 275305224
// solutions:
//
// Number of solutions: 275305224   Number of backtracks: 10806920438401
// Elapsed time: 94:54:16
//
///////////////////////////////////////////////////////////////////////////////

local Test = None | A1gtA5 |A1gtE1_A5gtE1 | A1gtE5 |
RowFull | ColFull|DiagsIncomplete|
RowIncomplete | ColIncomplete | RowColIncomplete |
ColFullRowIncomplete | ColFullRowFull | RowFullColIncomplete |
DiagIncomplete | DiagFull | InvDiagIncomplete | InvDiagFull

local Square = [0..4]->[0..4]->I[0..25]
local AvailableNumbers = [0..24]->>I[0..25]
local AllPosInfo = [0..24]->>(row:I[0..4],col:I[0..4],test:Test)

///////////////////////////////////////////////////////////////////////////////
//
pred MagicSquares5x5All() iff
// Initialize the "square", i.e. the 5x5 rectangle which we want to fill with
// the numbers 1..25. We intialize all position as 0 (no numbers placed yet).
square:.Square & Dupl(5,Dupl(5,0),square) &
CalculateMagicSquare5x5(square) & PrintFormattedSolution(square,0)

///////////////////////////////////////////////////////////////////////////////
//
// This program uses the following order for placing the numbers on the square:
//
//       1 10 12 11  2
//      20  6 18  8 21
//      22 14  5 17 23
//      24  9 19  7 25
//       3 13 16 15  4
//
//  a1 > a5 and a1 > e1 and a1 > e5 and a5 > e1
//
// For example, this will happen after the first few numbers are placed:
//
//
// Position 1 at (0,0): place 'a1' number: No test neccessary
//
//  a1  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//
// Position 2 at (0,4): place 'a1' number: Test for a1 > a5.
//
//  a1  .  .  . a5
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//
// Position 3 at (4,0): place 'e1' number: Test for a1 > e1 & a5 > e1
//
//  a1  .  .  . a5
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//  e1  .  .  .  .
//
// Position 4 at (4,4): place 'e5' number: Test for a1 > e5 & 26 <= a1+a5+e1+e5 <= 78
//
//  a1  .  .  . a5
//   .  .  .  .  .
//   .  .  .  .  .
//   .  .  .  .  .
//  e1  .  .  . e5
//
//
// Position 5 at (2,2): place 'c3' number: Test for a1+c3+e5 < 65 and e1+c3+e5 < 65 &
//                                         52 <= a1+a5+e1+e5+2*c3 <= 104
//  a1  .  .  . a5
//   .  .  .  .  .
//   .  . c3  .  .
//   .  .  .  .  .
//  e1  .  .  . e5
//
//
// Position 6 at (1,1): place 'b2' number: Test for a1+b2+c3+e5 < 65
//
//  a1  .  .  . a5
//   . b2  .  .  .
//   .  . c3  .  .
//   .  .  .  .  .
//  e1  .  .  . e5
//
//
// Position 7 at (3,3): place 'd4' number: Test for a1+b2+c3+d4+e5 = 65
//
//  a1  .  .  . a5
//   . b2  .  .  .
//   .  . c3  .  .
//   .  .  . d4  .
//  e1  .  .  . e5
//
//
// Position 8 at (1,3): place 'b4' number: Test for e1+c3+b4+a5 < 65
//
//  a1  .  .  . a5
//   . b2  . b4  .
//   .  . c3  .  .
//   .  .  . d4  .
//  e1  .  .  . e5
//
//
// Position 9 at (3,1): place 'd2' number: Test for e1+d2+c3+b4+a5 = 65
//
//  a1  .  .  . a5
//   . b2  . b4  .
//   .  . c3  .  .
//   . d2  . d4  .
//  e1  .  .  . e5
//
//
// Position 10 at (0,1): place 'a2' number: Test for a1+a2+a5 < 65 & a2+b2+d2 < 65
//
//  a1 a2  .  . a5
//   . b2  . b4  .
//   .  . c3  .  .
//   . d2  . d4  .
//  e1  .  .  . e5
//
// Position 11 at (0,3): place 'a4' number: Test for a1+a2+a4+a5 < 65 & a4+b4+d4 < 65
//
//  a1 a2  . a4 a5
//   . b2  . b4  .
//   .  . c3  .  .
//   . d2  . d4  .
//  e1  .  .  . e5
//
// Position 12 at (0,2): place 'a3' number: Test for a1+a2+a3+a4+a5 < 65
//
//  a1 a2 a3 a4 a5
//   . b2  . b4  .
//   .  . c3  .  .
//   . d2  . d4  .
//  e1  .  .  . e5
//
//  ....
//  ....
//  ....
//  ....
//
//
// Position 24 at (3,0): place 'd1' number: Test for a1+b1+c1+d1+e1 = 65 & d1+d2+d3+d4 < 65
//
//  a1 a2 a3 a4 a5
//  b1 b2 b3 b4 b5
//  c1 c2 c3 c4 c5
//  d1 d2 d3 d4  .
//  e1 e2 e3 e4 e5
//
// Position 25 at (3,4): place 'd5' number: Test for a1+b1+c1+d1+e1 = 65 & d1+d2+d3+d4+d5 = 65
//
//  a1 a2 a3 a4 a5
//  b1 b2 b3 b4 b5
//  c1 c2 c3 c4 c5
//  d1 d2 d3 d4 d5
//  e1 e2 e3 e4 e5
//
local pred CalculateMagicSquare5x5(square:.Square) iff
allpos :> AllPosInfo & allpos  =  [
(0,0,None),                     //1
(0,4,A1gtA5),                   //2
(4,0,A1gtE1_A5gtE1),            //3
(4,4,A1gtE5),                   //4
(2,2,DiagsIncomplete),          //5
(1,1,DiagIncomplete),           //6
(3,3,DiagFull),                 //7
(1,3,InvDiagIncomplete),        //8
(3,1,InvDiagFull),              //9
(0,1,RowColIncomplete),         //10
(0,3,RowColIncomplete),         //11
(0,2,RowFull),                  //12
(4,1,RowColIncomplete),         //13
(2,1,ColFull),                  //14
(4,3,RowColIncomplete),         //15
(4,2,RowFullColIncomplete),     //16
(2,3,ColFullRowIncomplete),     //17
(1,2,RowColIncomplete),         //18
(3,2,ColFullRowIncomplete),     //19
(1,0,RowColIncomplete),         //20
(1,4,RowFullColIncomplete),     //21
(2,0,RowColIncomplete),         //22
(2,4,RowFullColIncomplete),     //23
(3,0,ColFullRowIncomplete),     //24
(3,4,ColFullRowFull)            //25
] &

// We use a set of numbers 1..25 to fill the square. The sequential order of
// the numbers was intentionally chosen to correspond to a a valid solution.
// This is done in order to allow us to verify our filling sequence step-by-step
// easily...
numbers :. AvailableNumbers &
numbers := [25,24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1] &

// Start the actual computation
CoverAllSquares(square,numbers,allpos,0)

///////////////////////////////////////////////////////////////////////////////
//
// Cover all 25 cells
//
///////////////////////////////////////////////////////////////////////////////
local pred
CoverAllSquares(square:.Square,numbers:.AvailableNumbers,apos:<AllPosInfo,cnt:<I) iff
if cnt < 25 then
row = apos(cnt).row & col = apos(cnt).col &
PlaceOnePiece(square,numbers,row,col,apos(cnt).test) &

// Uncomment the following line to see the filling sequence in slow motion...
//      PrintFormattedSolution(square,0)& Pause(2) &
CoverAllSquares(square,numbers,apos,cnt+1)
end

///////////////////////////////////////////////////////////////////////////////
// The position at row,col is unused.
// Place there one number and perform a position dependent test.
///////////////////////////////////////////////////////////////////////////////
local pred
PlaceOnePiece(square:.Square,numbers:.AvailableNumbers, row:<I[0..4],col:<I[0..4],test:<Test) iff
p:>I[0..24] &                       // get an unused number index
number = numbers(p) & number <> 0 & // get an unused available number
square(row,col):= number &

// Perform the position dependent test.
case test of
None => true;
A1gtA5 => number < square(0,0);
A1gtE1_A5gtE1 => number < square(0,0) & number < square(0,4);
A1gtE5 => number < square(0,0) &
num = square(0,0) + square(0,4) + square(4,0) + number &
num <= 78 & num >= 26;
ColFullRowFull => CheckRowColFull(square,row,col);
RowFull => CheckRowFull(square,row);
ColFull => CheckColFull(square,col);
RowIncomplete => CheckRowIncomplete(square,row);
ColIncomplete => CheckColIncomplete(square,col);
RowColIncomplete => CheckRowColIncomplete(square,row,col);
DiagsIncomplete => CheckDiagIncomplete(square) & CheckInvDiagIncomplete(square) &
num = square(0,0) + square(0,4) + square(4,0) + square(4,4) + number + number &
num <= 104 & num >= 52;
ColFullRowIncomplete => CheckColFull(square,col) & CheckRowIncomplete(square,row);
RowFullColIncomplete => CheckRowFull(square,row) & CheckColIncomplete(square,col);
DiagFull => CheckDiagFull(square);
InvDiagFull => CheckInvDiagFull(square);
DiagIncomplete => CheckDiagIncomplete(square);
InvDiagIncomplete => CheckInvDiagIncomplete(square);
end &

// If we came here, we passed the test.
// mark the number we just placed as used so it will not be used again in this
// magic square.
numbers(p) := 0

///////////////////////////////////////////////////////////////////////////////
//
local pred CheckRowColIncomplete(square:.Square, row:<I[0..4], col:<I[0..4])  iff
square(row,0) + square(row,1) + square(row,2) + square(row,3) + square(row,4) < 65 &
square(0,col) + square(1,col) + square(2,col) + square(3,col) + square(4,col) < 65

local pred CheckRowColFull(square:.Square, row:<I[0..4], col:<I[0..4])  iff
square(row,0) + square(row,1) + square(row,2) + square(row,3) + square(row,4) = 65 &
square(0,col) + square(1,col) + square(2,col) + square(3,col) + square(4,col) = 65

local pred CheckRowFull(square:.Square, row:<I[0..4])  iff
square(row,0) + square(row,1) + square(row,2) + square(row,3) + square(row,4) = 65

local pred CheckColFull(square:.Square, col:<I[0..4])  iff
square(0,col) + square(1,col) + square(2,col) + square(3,col) + square(4,col) = 65

local pred CheckRowIncomplete(square:.Square, row:<I[0..4])  iff
square(row,0) + square(row,1) + square(row,2) + square(row,3) + square(row,4) < 65

local pred CheckColIncomplete(square:.Square, col:<I[0..4])  iff
square(0,col) + square(1,col) + square(2,col) + square(3,col) + square(4,col) < 65

local pred CheckDiagFull(square:.Square) iff
square(0,0) + square(1,1) + square(2,2) + square(3,3) + square(4,4) = 65

local pred CheckDiagIncomplete(square:.Square ) iff
square(0,0) + square(1,1) + square(2,2) + square(3,3) + square(4,4) < 65

local pred CheckInvDiagFull(square:.Square) iff
square(0,4) + square(1,3) + square(2,2) + square(3,1) + square(4,0) = 65

local pred CheckInvDiagIncomplete(square:.Square) iff
square(0,4) + square(1,3) + square(2,2) + square(3,1) + square(4,0) < 65

///////////////////////////////////////////////////////////////////////////////
// Print the solution in a more presentable way
local proc PrintFormattedSolution(square:<Square,row:<I) iff
if row < Len(square) then
Print('\n') & PrintOneRow(square(row),0) & PrintFormattedSolution(square,row+1)
else
Print('\n')
end

local proc PrintOneRow(r:<[0..]->I,col:<I) iff
if col < Len(r) then
if r(col) < 10 then
Print(' ')
end &
Print(' ',r(col)) & PrintOneRow(r,col+1)
end

```