F1 Source Code: MagicSquares5x5-All.f1

///////////////////////////////////////////////////////////////////////////////
//
// 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
This page was created by F1toHTML