Разбиения ко-лексикографически без рекурсии (PARI/GP)

[kthxbye]

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

Ниже представлены несколько разбиений в ко-лексикографическом порядке:

$\begin{bmatrix} 1 \end{bmatrix}$

$\begin{bmatrix} 1 & 1 \\ 2 \end{bmatrix}$

$\begin{bmatrix} 1 & 1 & 1\\ 2 & 1 \\ 3 \end{bmatrix}$

$\begin{bmatrix} 1 & 1 & 1 & 1\\ 2 & 1 & 1 \\ 3 & 1 \\ 2 & 2 \\ 4 \end{bmatrix}$

$\begin{bmatrix} 1 & 1 & 1 & 1 & 1 \\ 2 & 1 & 1 & 1 \\ 3 & 1 & 1 \\ 2 & 2 & 1 \\ 4 & 1 \\ 3 & 2 \\ 5 \end{bmatrix}$

Здесь прослеживается интересный паттерн. Мы можем взять лишь

$$\begin{bmatrix} 1 \\ 2 \\ 3 \\ 2 & 2 \\ 4 \\ 3 & 2 \\ 5 \end{bmatrix}$$

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

Идея не нова и была представлена в последовательности A210941 еще в далеком 2012.

Мне стало интересно, как генерировать члены последовательности этих разбиений без рекурсии. Вот что у меня получилось:

\\ Основная функция это list(vec).
\\ Входные данные это массив индексов k.
\\ Выходные данные это массив соответствующих разбиений.
\\ Здесь b4(k) вычисляется без рекурсии.
\\ Так что, эта программа хороша, когда вам нужно вычислить небольшое число разбиений.
\\ В особенности для значительно больших k.

b1(n) = {my(A = 1, B, v1); v1 = [1];
until(v1[A]>n, v1 = concat(v1, 0); B = 1;
while((binomial(B+1, 2) - binomial(B\2+1, 2)) <= A,
v1[A+1] += (-1)^((B-1)\2)*v1[A - binomial(B+1, 2) + binomial(B\2+1, 2) + 1]; B++); A++);
A-1}
list(vec) = {my(b2(n) = my(v1);
v1 = vector(n, i, vector(i\2+1, j, i==1));
for(i = 2, n, v1[i][i\2+1] = 1; v1[i][1] = vecsum(v1[i-1]);
for(j = 2, i\2, v1[i][j] = v1[i-1][j-1] - if((i-j) >= 2*(j-1), v1[i-j][j-1])));
for(i = 2, n, forstep(j = i\2, 1, -1, v1[i][j] += v1[i][j+1]));
v1,
vv1 = b2(b1(vecmax(vec))),
b3(n) = my(A = 0); until(vv1[A][1] > n, A++); A-1,
b4(n) = my(A = b3(n),
n = if(n == vv1[A][1], n, vv1[A][1] + vv1[A+1][1] - n),
B = n - vv1[A][1], C, v1);
v1 = []; while(!(B == 0), C = 1;
until(vv1[A+1][C]<=B, C++); v1 = concat(v1, C-1);
n -= vv1[A][1] - vv1[A - C + 1][1] + vv1[A+1][C];
A = b3(n); B = n - vv1[A][1]); v1 = Vecrev(v1);
v1 = concat(A + (#v1 > 0), v1));
vec = vector(#vec, i, b4(vec[i]))}

Мне бы хотелось как-то оценить скорость работы этого алгоритма. Для этого, его желательно с чем-нибудь сравнить. Интересуют алгоритмы без рекурсии.

Если интересует история разработки алгоритма, то буду рад поделиться деталями.

[wrest]

kthxbye в сообщении #1634684 писал(а):

Мне бы хотелось как-то оценить скорость работы этого алгоритма.

Много текста, а покажите как вы его запускаете и что получаете.

[kthxbye]

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

kthxbye в сообщении #1634684 писал(а):

Мне бы хотелось как-то оценить скорость работы этого алгоритма.

Много текста, а покажите как вы его запускаете и что получаете.

Ваша рекомендация кажется мне несколько странной, ну да ладно.

Запускаю например так:

v1 = list(vector(20, i, i))
for(i=1, 20, print(v1[i]))

Получаю следующее:

[1]
[2]
[3]
[2, 2]
[4]
[3, 2]
[5]
[2, 2, 2]
[4, 2]
[3, 3]
[6]
[3, 2, 2]
[5, 2]
[4, 3]
[7]
[2, 2, 2, 2]
[4, 2, 2]
[3, 3, 2]
[6, 2]
[5, 3]

Что же мы получили? Мы получили строки треугольника A210941 с индексами от $1$ до $20$ по возрастанию.

Но плюс алгоритма в том, что можно в качестве входных данных задавать значения не по-порядку, а абсолютно любые. Например:

v1 = list(vector(20, i, 2^i))
for(i=1, 20, print(v1[i]))

Соответственно он выдает следующее:

[2]
[2, 2]
[2, 2, 2]
[2, 2, 2, 2]
[4, 2, 2, 2]
[4, 3, 3, 2]
[11, 3]
[3, 3, 3, 3, 3, 2]
[5, 3, 3, 3, 2, 2, 2]
[8, 5, 2, 2, 2, 2, 2]
[4, 4, 3, 3, 3, 3, 2, 2, 2]
[17, 8, 2, 2]
[10, 8, 6, 4, 4]
[21, 11, 2, 2]
[8, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2]
[13, 11, 3, 3, 3, 3, 2, 2, 2, 2]
[12, 10, 5, 3, 3, 3, 3, 3, 2, 2, 2]
[14, 10, 9, 9, 6, 2, 2]
[18, 8, 8, 7, 5, 5, 5]
[17, 14, 10, 8, 4, 4, 2, 2]

Это строки треугольника A210941, индексы которых это степени двойки от $1$ до $20$ по возрастанию.

Другой пример:

v1 = list([100, 27, 99, 91, 63, 176, 2301, 453, 921, 6817])
for(i=1, 10, print(v1[i]))

Соответствующие строки треугольника A210941:

[7, 6]
[3, 3, 3]
[8, 5]
[11, 2]
[8, 2, 2]
[15]
[10, 7, 3, 3, 3]
[7, 3, 3, 3, 3]
[15, 5, 2]
[9, 8, 7, 7]

[Null]

Если я не ошибся можно из $a_1,\dots,a_k$(включая единицы) получить следующее так:
Если $a_k=a_{k-1}=\dots=a_{k-l+1}$($l$ штук), то
если $l>1$ то следующее $a_1,\dots,a_{k-l+1}+1$ и $a_k(l-1)-1$ единиц на конце.
если $l=1$ и $a_{k-1}=\dots=a_{k-m}$($m$ штук)
то следующее $a_1,\dots,a_{k-m}+1$ и $a_k+a_{k-1}(m-1)-1$ единиц на конце.

[Null]

Вы хотите получить строку по номеру?

[kthxbye]

Null в сообщении #1634748 писал(а):

Вы хотите получить строку по номеру?

Именно так. Мой алгоритм с этим успешно справляется, но я не уверен что он работает достаточно быстро.

[Null]

Идея такая(тут с учетом единичных слагемых):
В начале создаем таблицу $a_{i,j}$ - количество способов представить $i$ в виде суммы натуральных невозрастающих слагемых $<j$ - сложность $O(n^2)$
$a_{i,j}=a_{i,j-1}+[i\ge j-1]a_{i-j+1,j}$
$a_{0,j}=1$
$a_{i,1}=0$
$a_{0,1}=1$
Тогда $m$-тый элемент последовательности считается в цикле(Сложность O(n) на самом деле):
Если $a_{n,l+1}\ge m >a_{n,l}$, то первый элемент последовательности - $l$,
и далее переходим к $m:=m-a_{n,l};n=n-l$

У вас написано что-то похожее.
Ускорить можно если научиться считать $a_{i,j}$ быстро.

[wrest]

Null в сообщении #1634979 писал(а):

Ускорить можно если научиться считать $a_{i,j}$ быстро.

Тут у меня припасена функция, вычисляющая количество разбиений post1405311.html#p1405311
называется N10_2(n,k)

[kthxbye]

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

Тут у меня припасена функция, вычисляющая количество разбиений post1405311.html#p1405311
называется N10_2(n,k)

Замечательная функция! А можно ли ее как-нибудь немного переделать? В моем алгоритме используются значения

$$T(n, k) = \sum\limits_{i=k}^{n}R(n,i)$$

где $R(n, k)$ это A026794. Там очень похожая рекурсия.

[wrest]

kthxbye в сообщении #1635018 писал(а):

А можно ли ее как-нибудь немного переделать?

Конечно, переделывайте, тем более что я сам взял её из Иртернета.

Я переделать не могу, поскольку традиционно не понимаю чего вы хотите

[wrest]

kthxbye в сообщении #1635018 писал(а):

где $R(n, k)$ это A026794
.

Это я могу написать,
R(n,k)=if(n==k,1,if(k*2>n,0,#partitions(n-k,[k,n-k])))

kthxbye в сообщении #1635018 писал(а):

В моем алгоритме используются значения

$$T(n, k) = \sum\limits_{i=k}^{n}R(n,i)$$

А, ну и это:
T(n,k)=sum(x=k,n,if(n==k,1,if(k*2>n,0,#partitions(n-k,[k,n-k]))))