Перестановки. Примеры генераторов магических квадратов 4x4 и 5x5

Практически любая простая задача по составлению алгоритма усложняется, если потребовать достижения максимальной скорости. В первой части статьи мы рассмотрели генератор магических квадратов, который для порядка n=4 работал около минуты. Сейчас я покажу аналогичный генератор, но делающий ту же задачу за 0.01 секунды. Приведу полный текст программы (Free Pascal):

uses crt,dos;

var p:array[1..100] of integer;
    n,l,k,i:integer;
    c:longint;
    ti:longint;
    sum,nn:integer;
    stop:boolean;
    f:text;

function key:char;
begin
  if keypressed then key:=readkey else key:=#0;
end;

procedure time(x:byte);
var h,m,s,s100:word;
    r:real;
begin
  gettime(h,m,s,s100);
  if x=0 then ti:=s100+(s+m*60+h*longint(3600))*100 else begin
    r:=(s100+(s+m*60+h*longint(3600))*100-ti)/100;
    writeln('Time: ',r:0:2,' sec');
  end;
end;

procedure act;
begin
  inc(c);writeln(f,c,':');
  writeln(f, p[1]:l, p[2]:l, p[3]:l, p[4]:l);
  writeln(f,p[13]:l, p[9]:l, p[7]:l,p[14]:l);
  writeln(f,p[15]:l, p[8]:l,p[10]:l,p[16]:l);
  writeln(f, p[5]:l,p[11]:l,p[12]:l, p[6]:l);
end;

{ a1  a2  a3  a4
  a5  a6  a7  a8
  a9 a10 a11 a12
 a13 a14 a15 a16
a1 a2 a3 (a4) a13 (a16) a7 (a10) a6 (a11) (a14) (a15) a5 (a8) (a9) (a12)
 1  2  3   4    5    6   7    8   9   10    11    12  13  14   15    16  }

function good(x:integer):integer;
var s:integer;
begin
  case x of
    4:begin s:=sum-(p[1]+p[2]+p[3]);
      if s < p[1] then good:=0 else good:=s end;
    6:begin s:=sum-(p[1]+p[4]+p[5]);
      if (s < p[1])or(p[5] < p[1])or(p[5] < p[4]) then good:=0 else good:=s end;
    8:good:=sum-(p[4]+p[5]+p[7]);
   10:good:=sum-(p[1]+p[9]+p[6]);
   11:good:=sum-(p[2]+p[9]+p[8]);
   12:good:=sum-(p[3]+p[7]+p[10]);
   14:good:=sum-(p[13]+p[9]+p[7]);
   15:good:=sum-(p[1]+p[13]+p[5]);
   else good:=sum;
  end;
end;

procedure perm(x:integer);
var i,d,q:integer;
begin
  if stop then exit;
  if (x=nn) then begin act;if key=#27 then stop:=true;exit end;
  q:=good(x);if q=0 then exit;
  d:=p[x];
  if q=sum then
  for i:=x to nn do begin
    p[x]:=p[i];p[i]:=d;
    perm(x+1);
    p[i]:=p[x];p[x]:=d;
  end else
  for i:=x to nn do if p[i]=q then begin
    p[x]:=p[i];p[i]:=d;
    perm(x+1);
    p[i]:=p[x];p[x]:=d;
    break;
  end
end;

begin
  writeln('Magic 4x4');
  n:=4;nn:=n*n;c:=0;
  for i:=1 to nn do p[i]:=i;
  if paramstr(1)<>'' then begin
    assign(f,paramstr(1));reset(f);
    if ioresult <> 0 then exit;
    for i:=1 to nn do read(f,p[i]);
    close(f);
  end;
  sum:=0;k:=0;
  for i:=1 to nn do begin
    sum:=sum+p[i];
    if p[i]>k then k:=p[i]
  end;
  l:=2;while k>=10 do begin inc(l);k:=k div 10 end;
  if (sum mod n)<>0 then exit;
  sum:=sum div n;
  writeln('Summa=',sum);
  assign(f,'mag4x4.txt');rewrite(f);writeln(f,'Summa=',sum);
  stop:=false;
  time(0);
  perm(1);
  time(1);
  close(f);
  writeln('Count=',c);
end.
Процедура перебора осталась практически прежней, только слегка изменилось отсечение вариантов. И это "слегка" дало такой неплохой результат. Хочу сразу сказать, что этот вариант еще не лучший - просто хотелось сохранить относительную универсальность.

Первое, что необходимо отметить, это последовательность вычислений параметров квадрата. Мы стремимся найти квадрат
  a1  a2  a3  a4
  a5  a6  a7  a8
  a9 a10 a11 a12
 a13 a14 a15 a16
где числа a1, a2, ..., a16 берутся из заданного набора. Расположение этих чисел должно удовлетворять так называемым требованиям магичности, т.е. некоторым линейным соотношениям типа a1+a2+a3+a4=sum. Например, выбрав первые три числа a1, a2, a3, можно сразу вычислить a4 и просто проверить его наличие в исходном наборе. Среди 16 чисел квадрата только 7 являются независимыми, остальные вычисляются. В приведенной программе принята следующая последовательность определения элементов:
a1 a2 a3 (a4) a13 (a16) a7 (a10) a6 (a11) (a14) (a15) a5 (a8) (a9) (a12)
Здесь в скобках показаны вычисляемые величины. Последовательность определения элементов можно выбирать по разному и от этого выбора зависит скорость работы. Как это делать? Мы сейчас не будем заниматься этой задачей, просто отметим ее.

Следующий важный момент, не решенный в приведенном алгоритме, это проверка наличия вычисленного числа среди неиспользованных чисел набора. Сейчас применен последовательный просмотр, что, конечно, нельзя считать удачным решением.

06.04.2011
Для магических квадратов порядка n=5 попробовал применить следующую схему определения элементов (в скобках - вычисляемые):
{
                   1  2   3   4   5
                  20  6  18  10  21
                  22  7  11  15  23
                  24  8  19  13  25
                  12  9  17  16  14

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)}

function good(x:integer):integer;
var s:integer;
begin
  good:=sum;
  case x of
    5:begin s:=sum-(p[1]+p[2]+p[3]+p[4]);
      if (s < p[1])or(p[4] < p[2]) then good:=0 else good:=s end;
    7:if p[6] < p[1] then good:=0;
    9:begin s:=sum-(p[2]+p[6]+p[7]+p[8]);
      if p[8] < p[1] then good:=0 else good:=s end;
   11:if p[10] < p[1] then good:=0;
   12:begin s:=sum-(p[5]+p[8]+p[10]+p[11]);
      if (s < p[1])or(s < p[5]) then good:=0 else good:=s end;
   14:begin s:=sum-(p[1]+p[6]+p[11]+p[13]);
      if (s < p[1])or(p[13] < p[1]) then good:=0 else good:=s end;
   16:good:=sum-(p[4]+p[10]+p[15]+p[13]);
   17:good:=sum-(p[12]+p[9]+p[16]+p[14]);
   19:good:=sum-(p[3]+p[18]+p[11]+p[17]);
   21:good:=sum-(p[20]+p[6]+p[18]+p[10]);
   23:good:=sum-(p[22]+p[7]+p[11]+p[15]);
   24:good:=sum-(p[1]+p[20]+p[22]+p[12]);
  end;
end;
Оказалось, что программа по скорости обошла ранее разработанные мной программы. Решил протестировать ее на стандартном наборе 1,2,...,25. Как известно, впервые полный перебор магических квадратов порядка 5 был осуществлен Р. Шрёппелем в 1973г. на компьютере PDP-10. Прогон занял около 100 ч. машинного времени. Протестировал - около 4 часов (1.9 Ггц - один процессор). Самое приятное - результат 68.826.306 квадрата подтвердился :-)
Дополнительно протестировал набор из простых чисел:
3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 101 103
Результат - 58257 квадратов за 16.4 мин. Теперь можно пытаться улучшать эти результаты.

09.04.2011
Первая попытка улучшения связана с небольшим изменением порядка проверки элементов генерируемого квадрата:
{        1   3  (5)   4    2
        22  10 (23)   7  (21)
        18 (16) (9) (14) (19)
       (24)  8 (25)  11   20
         6  15 (17)  13  (12)   }
function good(x:integer):integer;
var s:integer;
begin
  good:=sum;
  case x of
    3:if p[2] < p[1] then good:=0;
    5:if p[4] < p[3] then good:=0 else good:=sum-(p[1]+p[2]+p[3]+p[4]);
    7:if (p[6] < p[1])or(p[6] < p[2]) then good:=0;
    8:if p[7] < p[1] then good:=0;
    9:if p[8] < p[1] then good:=0 else good:=sum-(p[2]+p[6]+p[7]+p[8]);
   11:if p[10] < p[1] then good:=0;
   12:begin s:=sum-(p[1]+p[9]+p[10]+p[11]);
      if (s < p[1])or(p[11] < p[1]) then good:=0 else good:=s end;
   14:good:=sum-(p[4]+p[7]+p[11]+p[13]);
   16:good:=sum-(p[3]+p[10]+p[8]+p[15]);
   17:good:=sum-(p[6]+p[15]+p[12]+p[13]);
   19:good:=sum-(p[18]+p[16]+p[9]+p[14]);
   21:good:=sum-(p[20]+p[12]+p[19]+p[2]);
   23:good:=sum-(p[22]+p[7]+p[10]+p[21]);
   24:good:=sum-(p[1]+p[18]+p[22]+p[6]);
  end;
end;
Теперь из набора простых чисел, который рассматривался ранее, генерация всех 58257 квадратов прошла за 11.8 мин. Исходные тексты генераторов квадратов 5x5 и 4x4 и исполняемые программы лежат здесь.

10.04.2011
Очередная попытка ускорения работы программы - 8.4 мин. при генерации квадратов из простых чисел.
         1   3  (5)   4    2
        22  10 (23)   7  (21)
        18 (16) (9) (14) (19)
       (24)  8 (25)  11   20
         6  15 (17)  13  (12)  

procedure perm(x:byte);
var i,o:byte;q,d:longint;
begin
  if stop then exit;
  o:=2;q:=1;
  case x of
    2,7,8,10,11:begin q:=p[1];o:=1 end;
    4:begin q:=p[3];o:=1 end;
    5:q:=sum-(p[1]+p[2]+p[3]+p[4]);
    6:begin q:=p[2];o:=1 end;
    9:q:=sum-(p[2]+p[6]+p[7]+p[8]);
   12:begin q:=sum-(p[1]+p[9]+p[10]+p[11]);if (q < p[1]) then exit end;
   14:q:=sum-(p[4]+p[7]+p[11]+p[13]);
   16:q:=sum-(p[3]+p[10]+p[8]+p[15]);
   17:q:=sum-(p[6]+p[15]+p[12]+p[13]);
   19:q:=sum-(p[18]+p[16]+p[9]+p[14]);
   21:q:=sum-(p[20]+p[12]+p[19]+p[2]);
   23:q:=sum-(p[22]+p[7]+p[10]+p[21]);
   24:q:=sum-(p[1]+p[18]+p[22]+p[6]);
   25:begin act;if key=#27 then stop:=true;exit end;
   else o:=0;
  end;
  if q<=0 then exit;
  d:=p[x];i:=x;
  case o of
    0:repeat
        p[x]:=p[i];p[i]:=d;
        perm(x+1);
        p[i]:=p[x];p[x]:=d;
        inc(i)
      until i>nn;
    1:repeat
        if p[i]>q then begin
          p[x]:=p[i];p[i]:=d;
          perm(x+1);
          p[i]:=p[x];p[x]:=d;
        end;
        inc(i)
      until i>nn;
    2:repeat
        if p[i]=q then begin
          p[x]:=p[i];p[i]:=d;
          perm(x+1);
          p[i]:=p[x];p[x]:=d;
          exit;
        end;
        inc(i)
      until i>nn;
  end;
end;
Исходные тексты генераторов квадратов 5x5 и исполняемые программы лежат [здесь]. Получение полного набора стандартных квадратов на моей машине потребовало 1.31 часа, что уже приемлемо.



 
X