Конкурс Topswops

Имеется n карт, каждая из которых имеет свой номер a[i] от 1 до n:

a[1], a[2], ... , a[n]
Поворотом p(k) назовем следующее преобразование этой последовательности:
a[1], a[2], ... , a[n] -> a[k], a[k-1], ... , a[1], a[k+1], ..., a[n]
Для каждого начального набора {a} карт можно построить цепочку наборов {ai}, т.ч. {a0}={a} и {ai+1}=p(ai[1])({ai}). Последним набором цепочки считается такой набор, для которого ak[1]=1. Число k, равное числу поворотов начального набора, назовем длиной цепочки. Легко доказывается, что для любого начального набора длина порождаемой цепочки конечна. Задача заключается в нахождении таких наборов, которые порождают цепочки максимальной длины.

Пример максимальной цепочки для n=5:
   3 1 4 5 2
   4 1 3 5 2
   5 3 1 4 2
   2 4 1 3 5
   4 2 1 3 5
   3 1 2 4 5
   2 1 3 4 5
   1 2 3 4 5
Длина этой цепочки равна 7 (число поворотов). В конкурсе Topswops предлагалось найти наборы для всех простых n от 2 до 97. Доказательства максимальности не требовалось.
В своё время Д.Кнут опубликовал два алгоритма для поиска максимальных наборов. Приведу первый алгоритм на языке pascal:
const maxl=1000;

type perm=array[0..16] of integer;

var a,b:array[0..maxl] of perm;
    n:integer;
    x:array[0..maxl] of integer;

procedure main;
var m,l,j,k:integer;
label try;
begin
  a[0,0]:=1;m:=0;l:=0;
  try:
  x[l]:=x[l]+1;k:=x[l];
  if k0) then goto try else
    else
      if a[l,k]<>k+1 then goto try;
    a[l+1]:=a[l];
    for j:=1 to k do a[l+1,j]:=a[l,k-j];
    b[l+1]:=b[l];
    a[l+1,0]:=k+1;b[l+1,k+1]:=1;
    if l>=m then begin
      m:=l;
      write(m+1,': ');
      for j:=0 to n-1 do write(a[l+1,j],' ');writeln;
    end;
    l:=l+1;x[l]:=0;
    goto try;
  end;
  l:=l-1;
  if l>=0 then goto try;
end;

begin
  writeln('Knuth 1');
  write('N= ');readln(n);
  main;
end.
Здесь мы имеем множество наборов a[l,i] для уровней l=0,... Для каждого уровня имеется число x[l], равное максимальному индексу, для которого строится дерево, а также массив меток b[i] с обозначением уже использованных элементов (b[i]=1). Делается попытка перехода от уровня l к уровню l+1 - если a[l,k]=k+1 (заметим, что индексация карт идет от 0). При переходе к l=l+1 элемент x[l] обнуляется. Попытка перехода заключается в добавлении элемента k+1 на 0 место - это может быть впервые, если a[l,k]=0 и элемент k+1 еще не использовался, и повторно, если a[l,k]=k+1. В обоих случаях идет инверсия k предыдущего уровня. Если уже все элементы были использованы, то идет обычный обход дерева, в противном случае тестируются новые элементы. А т.к. добавление новых элементов происходит на разных ветках дерева, то происходит полный перебор ситуаций.

Обратим внимание на то, как происходит переход от n к n+1. После достижения некоторого узла, в котором находится набор из n элементов:
a[1], a[2], ..., a[n]
и для которого a[k]<>k для всех k, т.е. для концевого узла ветки дерева, производится "достройка" этого узла с помощью добавления нового числа n+1. Этот узел содержит n+1 элемент и выглядит так:
n+1, a[n], a[n-1], ..., a[1]
Ясно, что длина цепочки для этого набора на 1 больше длины цепочки набора {a} и есть шанс нарастить от этого узла новую ветку. Можно было бы взять и узел
a[1], a[2], ..., a[n], n+1
но для одного добавляемого числа получится эквивалентная ситуация.
Возможен еще один подход к переходу от n к n+1. Пусть мы имеем цепочку
   a0[1], a0[2], ..., a0[n]
   a1[1], a1[2], ..., a1[n]
   ...
   ak[1], ak[2], ..., ak[n]    ak[1]=1
Если теперь во всех наборах заменить 1 на n+1, а на n+1 место записать 1, то получим цепочку
   a0[1], a0[2], ..., a0[n], 1
   a1[1], a1[2], ..., a1[n], 1
   ...
   ak[1], ak[2], ..., ak[n], 1  ak[1]=n+1
которую, естественно, можно продолжить далее.
Как для первого варианта, так и для второго можно рассмотреть добавление не только одной карты, а нескольких. Но выращивать дерево по первому варианту в вычислительном отношении сложнее, чем продолжение цепочки по второму варианту для различных перестановок добавляемых элементов. К сожалению, для обоих вариантов остается проблема выбора набора {a}, т.к. практически всегда (скорее всего - всегда) этот набор не является максимальным.

П.п.п.

Далее мы рассмотрим представление наборов в виде так называемых первых первых последовательностей (п.п.п.). Это представление я узнал от Павловского (форум dxdy.ru). Рассмотрение последовательности первых элементов цепочек идет от статьи Morales L. and Sudborough, I.H. "A Quadratic Lower Bound for Reverse Card Shuffle." Presented at 26th S.E. Conf. Comb. Graph Th. Computing. Preprint 1995. Например для вышеприведенной цепочки n=5 эта последовательность выглядит так:
3, 4, 5, 2, 4, 3, 2, 1
Если в этой последовательности оставить только те числа, которые встречаются первый раз, то мы получим п.п. последовательность:
3, 4, 5, 2, 1
Самое приятное в п.п.п. это то, что они однозначно задают цепочку и начальный набор. Конечно, для полной строгости этого утверждения следовало бы потребовать, чтобы в п.п.п. присутствовали бы все используемые числа, хотя даже для случая неполного использования чисел цепочка восстанавливается однозначно, но с "дырками".

Два варианта перехода от n к n+1 на языке п.п.п. смотрятся совсем просто. Пусть
p[1], p[2], ..., p[n-1], 1
исходный набор из n элементов, тогда переход к набору из n+1 элементов производится
- по первому варианту
n+1, p[1], p[2], ..., p[n-1], 1
- по второму варианту
p[1], p[2], ..., p[n-1], n+1, 1
Второй вариант легко обобщается на k добавляемых элементов. Для любой перестановки q[1], ...,q[k] k элементов n+1, n+2, ..., n+k можно рассматривать обобщенный второй вариант
p[1], p[2], ..., p[n-1], q[1], ..., q[k], 1
в котором исходный набор полностью сохраняет свою цепочку. По первому же варианту для сохранения исходной цепочки возможно только следующее наращивание
n+k, ..., n+2, n+1, p[1], p[2], ..., p[n-1], 1
Понятно, что можно рассматривать и сочетания обоих вариантов.

На языке п.п.п. можно также рассматривить и цепочки. Приведем пример цепочки для n=5, которая рассматривалась ранее:
   3 4 5 2 1
   4 5 2 3 1
   5 2 4 3 1
   2 4 3 1
   4 3 2 1
   3 2 1
   2 1 
   1 
Переход к следующей п.п.п. происходит путем перестановки первого элемента на новое место при сохранении позиций остальных элементов. К концу цепочки размер п.п.п. начинает уменьшаться, т.к. часть элементов исключаются из процесса. Но эти элементы можно было бы и оставить:
   3 4 5 2 1
   4 5 2 3 1
   5 2 4 3 1
   2 4 3 1 5
   4 3 2 1 5
   3 2 1 4 5
   2 1 3 4 5
   1 2 3 4 5
Эта картинка напоминает приведенный алгоритм Д.Кнута - справа от 1 лежат числа, которые добавляются на первое место при построении дерева.
Приведу процедуру преобразования п.п.п. при переходе к следующему элементу цепочки
procedure next(var a:ppp);
var i,k:integer;aa:nabor;
    b:byte;
begin
  for i:=1 to n do aa[i]:=0;
  i:=0;b:=a[1];
  repeat
    if aa[1]=0 then begin inc(i);aa[1]:=a[i]; end;
    inv(aa[1],aa);
  until (aa[1]=b)or(aa[1]=1);
  for k:=1 to i-1 do a[k]:=a[k+1];a[i]:=b;
end;
В этой процедуре используется поворот обычного набора inv(aa[1],aa). Т.к. поворот фактически определяет быстродействие программы, то его желательно максимально оптимизировать. Я записал эту процедуру на ассемблере, который понимает компилятор Free Pascal:
procedure inv(k:longint;var a:nabor);assembler;
asm
    mov esi,a
    mov edi,esi
    add edi,k
    dec edi
@n: mov al,[esi]
    mov ah,[edi]
    mov [edi],al
    mov [esi],ah
    inc esi
    dec edi
    cmp edi,esi
    jg @n
end;
Понятно, что для различных машин возможны иные решения. Для предлагаемой задачи достаточно было использовать байты для задания элементов наборов. В наше время 8-битовые операции уже ушли в прошлое, поэтому может показаться, что можно было хотя бы использовать 16-битовые операции. Но это сопряжено с небольшими изменениями кода, которые, как это ни удивительно, только увеличивали время исполнения.

Конкурс Topswops после окончания был продолжен в виде новой задачи, которую сформулировал Tom Rokicki. Необходимо найти как можно более длинную цепь для N=10000. Для этой задачи не мудрствуя лукаво я переписал процедуру поворота в виде:
procedure inv(k:longint;var a:nabor);assembler;
asm
    mov esi,a
    mov edi,esi
    add edi,k
    add edi,k
    sub edi,2
@n: mov ax,[esi]
    mov bx,[edi]
    mov [edi],ax
    mov [esi],bx
    add esi,2
    sub edi,2
    cmp edi,esi
    jg @n
end;
т.к. уже необходимо работать с 16-битовыми элементами.

При отсутствии удовлетворительной теории для решения задачи приходится пользоваться перестановками карт. Например, допустим мы нашли последовательность для n. Добавим к этой последовательности еще k карт. При работе с п.п.п. мы просто дописываем эти новые карты в конец последовательности и переносим 1. Затем для различных перестановок этих k карт вычисляем длины цепочек и выбираем цепочку с максимальной длиной. При работе с n=10000 стало особенно заметно, что даже подобное простое наращивание требует тщательной оптимизации алгоритма. Каждый раз заново считать длину цепочки, которая превышает 200,000,000, уже непозволительная роскошь. Следующая процедура позволила частично справиться с задачей:
procedure fperm(i,x,y:integer;aa:nabor;c:longint);
var j,k:integer;b:word;
begin
  if stop then exit;
  if i=1 then begin
    for j:=1 to n do aa[j]:=0;
    c:=0;aa[1]:=a[1];
  end;
  while (aa[1] > 1) and (i < x) do begin
    inv(aa[1],aa);inc(c);
    if aa[1]=0 then begin inc(i);aa[1]:=a[i] end;
  end;
  if ic=x then begin
    ac:=aa;cp:=c;i_old:=i
  end;
  if x=y then begin
    while (aa[1] > 1) and (i < n) do begin
      inv(aa[1],aa);inc(c);
      if aa[1]=0 then begin inc(i);aa[1]:=a[i] end;
    end;
    if c>cm then begin
      cm:=c;write(#13,cm,':',y,' ');
      am:=a;
    end;
  end
  else begin
    if y-x>4 then if key=#27 then stop:=true;
    for j:=x to y do begin
      b:=a[x];a[x]:=a[j];a[j]:=b;aa[1]:=a[x];
      fperm(i,x+1,y,aa,c);
      b:=a[x];a[x]:=a[j];a[j]:=b;
    end;
  end;
end;
Здесь строится п.п.п. "a" (глобальная переменная), по ходу дела строится вспомогательная обычная последовательность "aa", заготовка которой передается в эту процедуру в качестве параметра. Также в эту процедуру передаются параметры длины "c" и индекс "i" для возможности продолжения строительства "aa". Параметры "x" и "y" задают отрезок "x, x+1, ..., y" индексов карт п.п.п. которые перемешиваются. Для передачи информации за пределы процедуры используются глобальные переменные "ac", "cp" и "i_old". Перебор по всем перестановкам осуществляется с помощью рекурсии, кусочек:
    for j:=x to y do begin
      b:=a[x];a[x]:=a[j];a[j]:=b;aa[1]:=a[x];
      fperm(i,x+1,y,aa,c);
      b:=a[x];a[x]:=a[j];a[j]:=b;
    end;
Наверное нужно было бы выбрать другой вариант перебора, например, лексикографический. Но, на самом деле, этот кусочек требует особого анализа, т.к. от него зависит скорость работы.

Результат зависит от количества k добавляемых карт. Число перестановок k! для k=7 уже составляет 5040 - это k потребовало для моей машины (1.9 Ггц) около 50 часов для достижения результата 203,216,214. Никакие повторные коррекции получаемой последовательности практически невозможны.

Что делать? Выход напрашивается сам - сократить количество просматриваемых перестановок. К сожалению вразумительных критериев подходящего набора перестановок пока нет.

На форуме AlZimmermannsProgrammingContests Dmitry Kamenetsky привел длины цепочек для n=9, 17, 25, ..., 3985, т.е. для k=8. Видимо из-за длительности процесса он не продолжил далее. Зарегистрированная им последовательность имеет длину 203,245,822, что чуть превышает полученный мной вариант при k=7. Это, кстати, не очень понятно, т.к., если бы он продолжил цепочку 3985 с шагом 7, то, как мне кажется, его результат был бы лучше. На первом месте пока Markus Sigg с результатом 254,151,933.

Как использовать процедуру fperm показано ниже:
procedure FPermutation;
var x,y,i,k,t,nn:integer;
begin
  window(1,5,80,15);clrscr;
  writeln('FPermutation');
  write('Number card=');readln(t);if t<2 then exit;
  for i:=1 to n do if a[i]=1 then begin k:=i;break end;
  y:=n-1;nn:=n;
  if t>=k then t:=k-1;
  if y>k-1 then y:=k-1;if y<2 then exit;
  x:=y+1-t;if x<1 then x:=1;
  am:=a;cm:=0;
  i_old:=1;cp:=0;
  write('From card(1..',x,')=');readln(i);
  stop:=false;
  time(0);
  repeat
    if i>x then i:=x;
    gotoxy(20,wherey);write(y+i-x,' ');
    ic:=i;n:=y+i-x+1;
    fperm(i_old,i,y+i-x,ac,cp);
    a:=am;ac[1]:=a[i];n:=nn;
    i:=i+t;
  until i=x+t;
  writeln;
  time(1);
  fout;
end;
Понятно, что эта часть программы не критична. Единственное требование для использования fperm заключается в том, чтобы параметр "i" ("i_old") возрастал при последовательных вызовах fperm. В этом случае возможно сокращение вычислений путем передачи вспомогательных параметров "i_old", "ac" и "cp". Для ограничения вычислений "хвостов", которые могут быть очень большими, применяется временная смена длины "n". При малых "n" я использовал многократные прогоны и не менял длину последовательности "n". Таким путем иногда удавалось сильно увеличивать итоговую длину цепочки, но объем вычислений "хвоста" (до достижения aa[]=1) был очень значительным, а эти вычисления необходимо проводить для каждой рассматриваемой перестановки. Если мы временно уменьшаем длину "n", то неявно предполагаем, что последующие элементы п.п.п. добавляют одно и то же число в результирующую длину цепочки для всех перестановок. Действительно, это справедиво, например, для следующей п.п.п.:
a[1], ..., a[x], a[x+1], ..., a[y], y+2, y+3, ..., n, 1

Если все a[i] в п.п.п. уменьшить на величину "x", то перемешиваемые участки
a[x], a[x+1], ..., a[y]
будут выглядеть как некоторые перестановки элементов
1, 2, ..., k
где k=y-x+1. Особый теоретический интерес представляет вопрос об использовании некоторой фиксированной перестановки в процессе наращивания п.п.п. Именно для таких п.п.п. была обнаружена квадратичная зависимость длины цепочки от "n" (см. упоминавшуюся выше статью Morales L. и Sudborough, I.H.).

В простейшем процессе наращивания по второму варианту итоговая последовательность выглядит как "сумма" перестановок:
[p[1,1],p[1,2],...,p[1,k1]] + [p[2,1],p[2,2],...,p[2,k2]] + ...
где каждое слагаемое [p[i,1],p[i,2],...,p[i,ki]] является некоторой перестановкой чисел [1,2,3,...,ki]. Каждое последующее слагаемое можно выбирать исходя из требования достижения максимума длины цепочки для получаемых частных последовательностей. Эксперименты показали, что такой способ дает неплохие результаты. Количество добавляемых карт ki на каждом шаге ограничено вычислительными ресурсами используемого компьютера.
Можно усовершенствовать этот способ, если строить частные последовательности, отталкиваясь от набора уже построенных последовательностей. При построении последовательности из "n" карт просто пытаемся добавлять 1,2,...,ki карт к предыдущим последовательностям из n-1,n-2,...,n-ki карт, оставляя для "n" вариант с максимальной длиной цепочки. Для каждого "n" можно оставлять не одну, а несколько последовательностей для использования в дальнейших построениях. Но при этом увеличивается объем работы.

Для n=10000 я попробовал использовать фиксированную перестановку:
[1]+[2 7 1 8 5 3 4 9 6]+[2 7 1 8 5 3 4 9 6]+...
Получившийся результат C=108.568.388 не так уж и плох, учитывая быстроту получения результата. Перестановка [2 7 1 8 5 3 4 9 6] была найдена из соображения максимума длины цепочки для короткой последовательности типа
[1]+20*[2 7 1 8 5 3 4 9 6]
Использование максимальных на каждом шаге перестановок из 8 карт должно давать результат близкий к 250.000.000. С ростом "n" скорость перебора убывает, поэтому приходится сокращать количество добавляемых карт либо использовать перебор по некоторому набору перестановок. Так до n=1000 я использовал перестановки из 9 карт, затем до n=6800 использовал максимальные перестановки из 8 карт. Дальше просто использовал уже найденный набор перестановок из 8 карт. Итог: C=228.725.337. Для анализа получемых длин цепочек удобно использовать величину
p(n)=length/n2

Эта величина для процесса добавления максимальных перестановок фиксированной длины слабо возрастает. Например при 9 картах p(1000)=2.5359. Переход на 8 карт привел к начальному снижению этой величины - p(6800)=2.4820. Дальше можно было предполагать рост этой величины, т.к. процесс снижения уже закончился. Переключение на перебор имеющихся перестановок из 8 карт привел к плавному уменьшению до p(10000)=2.2873.

Дополнение

Вычисление длины обыкновенной последовательности:
function cc(x:nabor):longint;
var i:longint;
begin
  i:=0;
  while x[1]>1 do begin inc(i);inv(x[1],x) end;
  cc:=i
end;
Вычисление длины п.п. последовательности:
function fcc(a:nabor):longint;
var i,c:longint;aa:nabor;
begin
  for i:=1 to n do aa[i]:=0;
  i:=1;c:=0;aa[1]:=a[1];
  if a[1]>1 then
  repeat
    inv(aa[1],aa);inc(c);
    if aa[1]=0 then begin inc(i);aa[1]:=a[i] end;
  until a[i]<=1;
  fcc:=c;
end;
Преобразование п.п.п. в обыкновенную последовательность:
procedure conv2nabor(var a:nabor;var c:longint);
var aa,b:nabor;i:integer;
begin
  for i:=1 to n do begin aa[i]:=0;b[i]:=i; end;
  i:=0;c:=0;
  repeat
    if aa[1]=0 then begin inc(i);aa[1]:=a[i] end;
    if aa[1]>1 then begin inv(aa[1],b);inv(aa[1],aa);inc(c) end;
  until aa[1]=1;
  for i:=1 to n do a[b[i]]:=aa[i];
end;
Преобразование обыкновенной последовательности в п.п.п.:
procedure conv2front(var a:nabor;var c:longint);
var aa,b:nabor;i:integer;
begin
  for i:=1 to nmax do begin aa[i]:=0;b[i]:=0; end;
  i:=0;c:=0;
  repeat
    if aa[a[1]]=0 then begin inc(i);b[i]:=a[1];aa[a[1]]:=1 end;
    if a[1]>1 then begin inv(a[1],a);inc(c) end;
  until (b[i]=1);
  a:=b;
end;

Результаты конкурса 10K

Tom Rokicki подвел итоги конкурса 16.03.2011:
Here's the scoreboard as of midnight tonight, end of contest.

276,790,959 Dmitry Kamenetsky
254,151,933 Markus Sigg
228,725,337 Serg Belyaev
217,111,538 Vincent Quesnoit
170,288,670 Mark Beyleveld
156,671,621 Kevin Burfitt
153,847,641 Herbert Kociemba
152,301,367 Jean Floch
150,177,788 Roy van Rijn
61,616,118 Jaroslaw Wroblewski

Looks like Dmitry wins! Congratulations to all who entered; I wasn't even
sure I was going to be able to score these things.

Thanks to everyone who entered!

-tom
На форуме AlZimmermannsProgrammingContests 17.03.2011 появилось письмо победителя конкурса Dmitry Kamenetsky. Даю свой вольный перевод основной части этого письма:

В моем подходе не очень много интересного и довольно много грубой силы. На самом деле он идентичен описанному Сергеем Беляевым. Я использовал "uniqueHeads" представление для перестановок - это список карт, которые появляются впервые на первом месте. Например, если оригинальный колода 3,1,4,2,5 то мы получаем следующую последовательность шагов:
3,1,4,5,2
4,1,3,5,2
5,3,1,4,2
2,4,1,3,5
4,2,1,3,5
3,1,2,4,5
2,1,3,4,5
1,2,3,4,5
UniqueHeads представлением будет 3,4,5,2,1. Это представление имеет то же самое число элементов, что и исходная колода и всегда заканчивается 1. Каждой колоде соответствует единственная uniqueHeads и наоборот. Не слишком трудно конвертировать одно представления в другое.

В своем подходе я итеративно строил uniqueHeads перестановки порциями длины k. На каждом этапе m (от 0), я перебирал все k! перестановок [2,3,4 ,...,k+1]+mk и оставлял перестановку, которая давала мне наибольшее количество очков. Ключом оптимизации было повторное использования перестановки (и оценки) от предыдущего (m-1)-го этапа для текущего этапа m. Я также использовал быстрый метод подсчета оценки, который был описан на этих форумах, но это дало мне только 3-х кратное ускорение и только для больших N. Мое последнее решение для 276M использовало k=9 для первые 8000 элементов и k=8 для остальных 2000 элементов. Это заняло около 2 недель вычислений.

Я заметил, что оценки растут линейно с ростом k, хотя время работы увеличивается в геометрической прогрессии. Так меня интересовало при каких больших k можно получить р-оценку 4+. Вот мои результаты для N=10K и различных k:
k=9: не закончен, но стремится к примерно 291M
k=8: 248M
k=7: 203M
k=6: 156M
k=5: 126M
При таких темпах я ожидаю, что k=12 даcт p-оценку 4+. Так что, возможно, стоит для этого стартовать продолжение проекта?

Я испытывал несколько других идей, которые не работают так хорошо, но, возможно, они будут работать у других:

* То же самое, но только фиксировать первое число в каждой стадии
* Вместо того, чтобы проверять все k! перестановок попробовать различные стратегии для сокращение их числа. Например k!/2 или k!/3 случайных перестановок. Это может быть сделано как на этом сайте. Я также пробовал фиксировать первые (скажем 3) цифры, а затем проверять (k-3)! перестановок
* Уменьшить k, но просматривать больший набор возможных значений. Например, для k=3 я в настоящее время перемешиваю [2,3,4]. Вместо этого вы могли бы перемешивать [2,3,4], [2,3,5], [2,4,5], [3,4,5]. Если мы выбираем из р значений, то мы должны рассмотреть Choose(p,k)*k! перестановок. Так хорошим компромиссом будет что-то подобное k=6 и p=11, аналогичное по сложности k=9. Во всяком случае я не проверял эту идею, но думаю, что она имеет некоторый потенциал.


1 2
 
X