Комбинаторное двенадцатизадачие: раскладываем шары по урнам

[arseniiv]

Можно кстати на общую мемоизацию не полагаться в случаях, когда отношение частичного порядка $A\prec B$ на аргументах функции, определённое как «$f(A)$ используется для вычисления $f(B)$», более-менее просто — тогда можно сразу написать код, вычисляющий нужные значения в нужном порядке и заканчивающий интересующим.

[maxal]

arseniiv, это называется динамическим программированием.

[arseniiv]

Да, знаю. Впрочем мне конечно стоило упомянуть, что это оно — не знаю почему не сделал.

[_Ivana]

Некоторые разделяют ДП и мемоизацию. А я думаю что это 2 стороны одной медали. И разница только в том, кто осуществляет патологическую топологическую сортировку аргументов - человек или пароход компьютер. Мне кажется, что проще переложить это с плеч программиста на плечи компьютера. (нескромная ссылка на мой стрим по этой теме https://www.youtube.com/watch?v=OnUXYVd-EuY)

[Geen]

wrest в сообщении #1402714 писал(а):

Geen в сообщении #1402697 писал(а):

А какое это имеет отношение к "шарам в урнах"?

Я не понял вопрос, сорри. Можете переформулировать?

Ну Вы же тему назвали "раскладывать шары по урнам", а не "12 простых комбинаторных зада на PARI/GP"...
(и, кстати, как я понимаю, в первом "решении" ошибка)

[arseniiv]

_Ivana
Я думаю, описанное мной — только частный случай ДП. Например есть динамическое программирование по профилю; я про него особо не думал и им не занимался, но на первый взгляд оно плохо должно превращаться в код с мемоизацией чего-нибудь.

[_Ivana]

Не сталкивался с ДП по профилю, но если там также итеративно преобразуется некоторая структура до определенных пор, то это точно так же можно реализовать с помощью рекурсивной функции, трактуя эту структуру как обобщенную мемоизацию (возможно, не значения функции на аргументах). Хотя конечно вы правы - чем гадать, лучше ознакомиться с явлением. Может при случае полюбопытствую.

[wrest]

Geen в сообщении #1402801 писал(а):

Ну Вы же тему назвали "раскладывать шары по урнам", а не "12 простых комбинаторных зада на PARI/GP"...

У вас претензия к названию темы?

Geen в сообщении #1402801 писал(а):

(и, кстати, как я понимаю, в первом "решении" ошибка)

Поправьте, пожалуйста.

[wrest]

В общем, я написал ещё один код, который вычисляет количество разбиений числа $n$ на не более чем $k$ частей по алгоритму, который уложен в исходных кодах pari/gp
Алгоритм описан там вот так:

(Оффтоп)

/* for loop over partitions of integer n.
* nbounds can restrict partitions to have length between nmin and nmax
* (the length is the number of non zero entries) and
* abounds restrict to integers between amin and amax.
*
* Start from central partition.
* By default, remove zero entries on the left.
*
* Algorithm:
*
* A partition of k is an increasing sequence v1,... vn with sum(vi)=n
* The starting point is the minimal k-partition of n: a,...a,a+1,.. a+1
* (a+1 is repeated r times with n = a * k + r).
*
* The procedure to obtain the next partition:
* - find the last index i<n such that v{i-1} != v{i} (that is vi is the start
* of the last constant range excluding vn).
* - decrease vi by one, and set v{i+1},... v{n} to be a minimal partition (of
* the right sum).
*
* Examples: we highlight the index i
* 1 1 2 2 3
* ^
* 1 1 1 3 3
* ^
* 1 1 1 2 4
* ^
* 1 1 1 1 5
* ^
* 0 2 2 2 3
* ^
* This is recursive in nature. Restrictions on upper bounds of the vi or on
* the length of the partitions are straightforward to take into account. */

Результат утешает не очень.
Для разбиений числа 70 на не более чем 30 частей, ответ 3 910 071 получается:
Рекуррентный алгоритм указанный в Википедии считает за 64 секунды, рекурсия делается 48 290 521 раз:

(Оффтоп)

$P(0,0)=1$
$P(i,0)=0$ для $i>0$
$P(n,k)=P(n,n)$ для $k>n$
$P(n,k)=P(n,k-1)+P(n-k,k)$
Реализованная так:

p(n,k)={
if(n==0& k==0, return (1));
if(k<1,return(0));
if(k>=n,return(numbpart(n)));
return(p(n,k-1)+p(n-k,k))}

Алгоритм из исходников pari/gp считает за 53 секунды, цикл проходится столько же раз какой и результат, т.е. 3 910 071.

(Оффтоп)

p2(n,k)={
my(v=vector(k),count=1,a=n\k,r=n%k,vi=0,n1=0);
for(i=1,k-r,v[i]=a);
for(i=k-r+1,k,v[i]=a+1);
while(#v>1,
k=#v;
vi=1;for(i=2,k-1,if(v[i-1]<>v[i],vi=i));
v[vi]--;n1=n-sum(i=1,vi,v[i]);k=k-vi;a=n1\k;r=n1%k;
for(i=vi+1,vi+k-r,v[i]=a);
for(i=vi+k-r+1,vi+k,v[i]=a+1);
if(v[1]==0,v=vecextract(v,[2..#v]));
count++);
return(count)}

Подсчет через #partitions() работает 0,3 секунды...

(Оффтоп)

? p2(70,30)
time = 52,542 ms.
%1 = 3910071
? p(70,30)
time = 1min, 3,945 ms.
%2 = 3910071
? #partitions(70,,[1,30])
*** partitions: Warning: increasing stack size to 16000000.
*** partitions: Warning: increasing stack size to 32000000.
*** partitions: Warning: increasing stack size to 64000000.
*** partitions: Warning: increasing stack size to 128000000.
*** partitions: Warning: increasing stack size to 256000000.
*** partitions: Warning: increasing stack size to 512000000.
*** partitions: Warning: increasing stack size to 1000000000.
time = 298 ms.
%3 = 3910071

В этой связи, если кто-то покажет как все-таки на pari/gp считать это быстрее, используя средства pari/gp и не слишком усложняя код, было бы здорово!

[wrest]

В общем, в pari/gp нашлась мемоизация!
С ней, подсчет количества разбиений работает мгновенно. Ускорение просто таки поражает воображение!

Однако, написать единую цельную функцию, которая бы ничего не требовала кроме двух входных чисел, считала бы количество раскладок и при этом не оставляла бы за собой никакого мусора (типа наполненного словаря), не вышло.

Надо делать глобальную переменную со словарём (ниже в коде -- M), эту переменную можно создавать и внутри функции и потом удалять, но если переменная существует до вызова функции, то она в таком случае будет удалена.

На случай если надо считать много разбиений, можно (и нужно) закрыть глаза на указанные недостатки. Но как коротенький сниппет, не канает. Ну может я ещё чего не усмотрел в доках pari/gp и можно как-то сделать функцию (которая собсно считает) внутри другой функции (где создается словарь который потом удаляется).

Код, использующий мемоизацию, считающий количество разбиений числа n на не более чем k частей:

M=Map();
memo_p(n,k)={
 my(res=0);
 if(n==0 && k==0,
   return(1));
 if(k<1,
   return(0));
 if(mapisdefined(M,[n,k],&res),
   return(res));
 if(k>=n,
   res=numbpart(n);
   mapput(M,[n,k],res);
   return(res));
 res=memo_p(n,k-1)+memo_p(n-k,k);
 mapput(M,[n,k],res);
 return(res)}

Работа:
? memo_p(1000,300)
%1 = 24060233693068843871790330541103
? ##
*** last result computed in 2,887 ms.

[maxal]

Тема перемещена из форума «Computer Science» в форум «Околонаучный софт»
Причина переноса: не указана.

[wrest]

maxal в сообщении #1402525 писал(а):

Вычислять числа Белла через сумму чисел Стирлинга тоже не самый оптимальный способ (альтернативы см. по ссылке
).

Посмотрел как это делают в OEIS, числа Белла: A000110.
Первый способ такой:
{a(n) = my(m); if( n<0, 0, m = contfracpnqn( matrix(2, n\2, i, k, if( i==1, -k*x^2, 1 - (k+1)*x))); polcoeff(1 / (1 - x + m[2, 1] / m[1, 1]) + x * O(x^n), n))}; /* Michael Somos */
Второй и третий способы очень медленные.
Четвертый способ - как раз через сумму чисел Стирлинга.
A000110(n)=sum(k=0, n, stirling(n, k, 2)) \\ M. F. Hasler, Nov 30 2018

Первый способ работает примерно в 2 раза быстрее, но начиная с $n \approx 1850$ становится недостаточно стека размером в один гигабайт.

Думаю, что оставлю сумму чисел Стирлинга, тем более что это более наглядно, как мне кажется.

[maxal]

В этой теме предлагается использовать модулярную арифметику для вычисления чисел Белла.

[wrest]

В связи со значительным прогрессом в программировании вычисления количества разбиений (partitions) целого числа $n$ на $k$ слагаемых, блок задач 10-12 переписываем так:

Цитата:

10. Сколько существует способов разложить $n$ одинаковых шаров в $k$ одинаковых урн?

N10(n,k)=#partitions(n,,[1,k])

Этот способ лаконичный, но для больших аргументов медленный или невыполнимый - т.к. создаётся массив разбиений, объем которого может быть очень велик, после чего подсчитывается количество элементов массива.
Следующие две функции подсчитывают требуемое количество намного быстрее и потребляют меньше памяти.
10.1 "Прямой" подсчет по рекуррентной формуле $P(0,0)=1;P(n>1,0)=0;P(n,k)=P(n,k-1)+P(n-k,k)$ с запоминанием промежуточных результатов. Для больших аргументов переполняется массив в котором запоминаются предыдущие вычисления.

(Оффтоп)

N10_1(n,k)=
  {
    local(M = Map());
    my(memo_p = (N,K)-> my(res);
       if(N==0 && K==0,
          return(1));
       if(K<1,
          return(0));  
       if(mapisdefined(M,[N,K],&res),
          return(res));
       if(K>=N,
          res=numbpart(N);
          mapput(M,[N,K],res);
          return(res));
       res=self()(N,K-1)+self()(N-K,K);
       mapput(M,[N,K],res);
       return(res));
    memo_p(n,k);
  }

10.2. "Хитрый" подсчет созданием матрицы размером [k+1] x [k] и последовательным вычислением, требует примерно n x k вычислений. Позаимствован в списке рассылки пользователей pari/gp у Karim Belabas ( http://pari.math.u-bordeaux.fr/archives ... 00015.html ) Промежуточная структура тоже может быть очень большой, но скорее ограничение тут по скорости работы:

(Оффтоп)

N10_2(n,k) =
  { 
my(k1 = k+1, kx = k1, a = 1, t = vector(k1, i, vector(k)));
  t[k1][1] = 1;
  for (i = 1, n-1,
    my(ky = a);
    t[a] = vector(k, j, if (!(ky--), ky = k1);
                        if (j == 1, t[ky][1]
                                      , t[kx][j-1] + t[ky][j]));
    kx = a; if (a++ > k1, a = 1));
  vecsum(t[a-1]);
}

Сравнительная скорость вычислений:

(Оффтоп)

? N10(70,30)
time = 203 ms.
%1 = 3910071
? N10_1(70,30)
time = 31 ms.
%2 = 3910071
? N10_2(70,30)
%3 = 3910071
? ##
*** last result computed in 0 ms.
Берем значения побольше. Первый метод -- не хватает памяти.
Два других:
? N10_1(1000,300)
time = 2,637 ms.
%4 = 24060233693068843871790330541103
? N10_2(1000,300)
time = 172 ms.
%5 = 24060233693068843871790330541103
? N10_2(10000,3000)
time = 28,407 ms.
%6 = 36167251325636291851786878829976856243927867494801874150751474019863050338273687745917678081066889663513003

11. Сколько существует способов разложить $n$ одинаковых шаров в $k$ одинаковых урн так, чтобы в каждой урне было не более одного шара (при условии $k\ge n$)?

N11(n,k)=if(k>=n,return(1),return(0))

12. Сколько существует способов разложить n одинаковых шаров в k одинаковых урн так, чтобы в каждой урне было не менее одного шара (при условии $n\ge k$)?

N12(n,k)=#partitions(n,,[k,k])

Замечания те же что и для задачи 10, так что только привожу "быстрый" код для решения:

(Оффтоп)

N12_2(n,k)=
  { 
my(k1 = k+1, kx = k1, a = 1, t = vector(k1, i, vector(k)));
  t[k1][1] = 1;
  for (i = 1, n-k-1,
    my(ky = a);
    t[a] = vector(k, j, if (!(ky--), ky = k1);
                        if (j == 1, t[ky][1]
                                      , t[kx][j-1] + t[ky][j]));
    kx = a; if (a++ > k1, a = 1));
  vecsum(t[a-1]);
}