Разбиения с ограничениями: кто быстрее?

[wrest]

Общее число разбиений $p(n)$ целого числа $n$ на слагаемые вычисляют (для больших $n$) по формуле Рамануджана-Харди (быстро сходящийся ряд).

Однако, если вводить ограничения на разбиения, например "слагемые не превышают k" или "число слагаемых не больше чем k", то формулы в виде быстро сходящегося ряда нет.

Например, в pari/gp есть встроенная функция, которая легко (на моем компе за примерно 12 секунд) вычисляет $p(10^9) \approx 1,6 \cdot 10^{35218}$

Однако, попытки вычислять количество разбиений $n$ на не более чем $k$ слагаемых уже при $n=10^5;k=10^4$ требуют намного бОльшего времени, а также огромных ресурсов по памяти. Для некоторых $k$ специального вида быстрые формулы существуют, но не в общем случае.

Изучение Википедии и смежных источников свелось к двум алгоритмам.
1. "Прямой" рекурсивный алгоритм.
$p(0,0)=1,p(n>1,0)=0$
$n>k:p(n,k)=p(n-1,k)+p(n-k,k)$
$n \le k:p(n,k)=p(n,n)$
Для исключения повторных вычислений одного и того же, этот алгоритм требует создания и наполнения массива размером примерно $nk$ элементов, так что при больших $n,k$ просто не хватает памяти. Ну и медленно работает.

p(n,k) =
{
    if(n==0 && k==0,return(0));
    local(NP=vector(k,i,numbpart(i)),M=vector(n,i,vector(k)));
    my(memo_p = (x,y)-> my(res);
       if(y<2,return(y));
       if(res=M[x][y],return(res));
       res=if(y>=x,NP[x],self()(x,y-1)+self()(x-y,y));
       M[x][y]=res; 
       return(res));
    memo_p(n,k);
}

2. "Хитрый" крученый алгоритм, сути которого я не понимаю, но который использует для вычисления массив размером около $k^2$ элементов, работает существенно быстрее. На pari/gp это выглядит так:

p_1(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[kx]);
}

Внутренний цикл работает примерно $nk$ раз.
"Хитрый" алгоритм быстрее "Прямого" примерно в $1,6$ раз:

? p(2000,1000)
*** vector: Warning: increasing stack size to 16000000.
*** vector: Warning: increasing stack size to 32000000.
time = 1,857 ms.
%1 = 4720819175618825183434073956853110605486740202
? > p_1(2000,1000)
*** vector: Warning: increasing stack size to 16000000.
time = 1,170 ms.
%2 = 4720819175618825183434073956853110605486740202
Для сравнения, вычисление всех $p(n,n)$ и их суммы от $n=1$ до $n=1000$
? sum(i=1,10000,numbpart(i));
time = 141 ms.

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

[wrest]

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

1. "Прямой" рекурсивный алгоритм.
$p(0,0)=1,p(n>1,0)=0$
$n > k:p(n,k)=p(n-1,k)+p(n-k,k)$
$n \le k:p(n,k)=p(n,n)$

Следует читать:
$\text{if } n\ge k \text{ then } p(n,k)=p(n,k-1)+p(n-k,k)$
$\text{if } n < k \text{ then } p(n,k)=p(n,n)$
С начальными условиями:
$p(0,0)=1,p(n>0,0)=0$

[Dmitriy40]

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

этот алгоритм требует создания и наполнения массива размером примерно $nk$ элементов, так что при больших $n,k$ просто не хватает памяти. Ну и медленно работает.

Для уменьшения памяти можно сохранять не все варианты комбинаций $n$ и $k$, а только лишь некоторые, например с условием $nk<10^8$ (так легко регулировать потребление памяти), остальные же продолжая вычислять рекурсивно. А ещё лучше - предрасчитать лишь небольшие значения $n$ и $k$, для прочих же класть в хеш-таблицу по ключу $nk$ (а лучше конечно по двум отдельно), тут ускорение оценить сложно, оно зависит от порядка вычислений и размера таблицы.

[wrest]

Dmitriy40 в сообщении #1405503 писал(а):

Для уменьшения памяти можно сохранять не все варианты комбинаций $n$ и $k$, а только лишь некоторые, например с условием $nk<10^8$ (так легко регулировать потребление памяти), остальные же продолжая вычислять рекурсивно.

Мне кажется, что оптимизировать "прямой" алгоритм особой выгоды не будет. Ну например, предвычисление всех $p(n)$ от 1 до $k$ по формуле Рамануджана дало ускорение в 6 раз (я первоначально не отправлял их на рекурсию, а вычислял через numbpart()), при том что не все они потом нужны (используется где-то 60-70%). Использование массива вместо списка дало ускорение ещё в два раза. Оказалось что помещение в список ( mapput() ) и поиск там (ismapdefined() )-- существенно затратнее по времени, массив наполняется в итоге не полностью (процентов на 80% наверное) и занимает места немного больше чем список, но выигрыш по времени существенный. Если менять массив обратно на список, будет выигрыш по памяти но, как я думаю, очень сильный проигрыш по скорости. Другое дело, что например при $n=10^4,k=3\cdot 10^3$ pari/gp для windows просто молча падает...
Тут именно в вопрос в применении какого-то другого алгоритма.
Ну например, "Хитрый" алгоритм взят примерно отсюда: https://code.jsoftware.com/wiki/Essays/ ... th_k_Parts и код на J выглядит так:

pnkv=: 4 : 0 " 0
 n=. x [ k=. y
 j=. <"1 (1+i.k) */ _1 1
 t=. ((1+k)*(-*n),1) {. ,: 0 1x
 for. i.(1<k)*0>.n-1 do. t=. (}.t), (0,j{t) + |.!.''{:t end.
 {:t
)

Это какая-то... шифрограмма... И там говорится что он быстро работает для малых $k$
А для больших $k$ (т.е. для $k$ поближе к $n$) есть другая шифрограмма

pnkd=: 4 : 0
  m=. y <. s=. x-y
  t=. >:i.s
  p=. 1,s$0x
  for_i. >:i.m do.
    for_j. (<:i)}.t do. p=. (+/(j,j-i){p)j}p end. end.
)

Но понять что она делает я не смог. Возможно её расшифровка и включение в "Хитрый" алгоритм может ускорить расчет для "больших" $k$, я бы по крайней мере проверил бы на каких $k$ алгоритмы сравниваются по скорости, и тогда в программе можно выбирать как считать в зависимости от значения $\dfrac{k}{n}$

[Dmitriy40]

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

Другое дело, что например при $n=10^4,k=3\cdot 10^3$ pari/gp для windows просто молча падает...

Это походу баг в PARI/GP, он неправильно инициализирует векторы (и матрицы и прочее) в памяти, пытается хранить только реально занятые ячейки, и иногда ошибается и идёт обращение к несуществующей ячейке. Для исправления достаточно сделать массив M глобальным (на ошибку это не влияет) и перед p(10000,3000) вызвать p(1000,10) и p(10000,10) (именно для сохранения массива между вызовами он и нужен глобальный), в таком порядке работает.
Не знаю как у Вас p(2000,1000) посчиталось за 1.8с, у меня ровно ваш код считает 10.5с, а мой 8.2с; ну а p(10000,3000) мой посчитал за три минуты:

(Код)

default(parisizemax,1G);
m=matrix(10000,5000);
{p(n,k)=my(x);
   k=min(k,n);
   if(n==0 || k==1, return(1));
   if(k==0, return(0));
   if(n<=matsize(m)[1] && k<=matsize(m)[2] && m[n,k]>0, return(m[n,k]));
   x=p(n,k-1)+p(n-k,k);
   if(n<=matsize(m)[1] && k<=matsize(m)[2], m[n,k]=x);
   return(x);
}
gettime();
printf("%u, %0.2fs\n",p(1000,10),gettime()/1000.)
printf("%u, %0.2fs\n",p(10000,10),gettime()/1000.)
printf("%u, %0.2fs\n",p(10000,3000),gettime()/1000.)

Вызов и вывод:
? \r pnk.gp
  ***   Warning: new maximum stack size = 1000000000 (953.674 Mbytes).
  *** matrix: Warning: increasing stack size to 16000000.
  *** matrix: Warning: increasing stack size to 32000000.
  *** matrix: Warning: increasing stack size to 64000000.
  *** matrix: Warning: increasing stack size to 128000000.
  *** matrix: Warning: increasing stack size to 256000000.
  *** matrix: Warning: increasing stack size to 512000000.
968356321790171, 0.08s
778400276435728381405745, 0.67s
36167251325636291851786878829976856243927867494801874150751474019863050338273687745917678081066889663513003, 172.65s

Но вообще согласен, алгоритм медленный.

[arseniiv]

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

Это какая-то... шифрограмма...

Потому что J.

[wrest]

Dmitriy40 в сообщении #1405561 писал(а):

Не знаю как у Вас p(2000,1000) посчиталось за 1.8с, у меня ровно ваш код считает 10.5с,

Ну не знаю, у меня p(2000,1000) (код -- копипаст из стартового поста) считает за 14 секунд на хилом андроид-планшете 4-летней давности:
? p(2000,1000)
*** vector: Warning: increasing stack size to 8000000.
*** vector: Warning: increasing stack size to 16000000.
time = 14,070 ms.
%131 = 4720819175618825183434073956853110605486740202
?
На ноуте (тоже хилом, примерно тех же лет) - в 10 раз быстрее.

-- 18.07.2019, 00:04 --

Dmitriy40 в сообщении #1405561 писал(а):

ну а p(10000,3000) мой посчитал за три минуты:

Видимо потому, что вы p(n,n) считаете рекурсивно, а я - встроенной функцией numbpart()
Ну и проверок многовато, плюс в целом вектор векторов работает немного быстрее чем матрица.

В общем, лишние круги нарезаете, замедляете вычисления :)

-- 18.07.2019, 00:06 --

arseniiv в сообщении #1405563 писал(а):

Потому что J.

Можете вторую расшифровать, которая pnkd ?

[Dmitriy40]

О, я понял: Вы запускаете 32-х битную версию gp.exe, а я 64-х битную. 32-х битная у меня тоже оказывается быстрее, 2.6с на вашем коде, вместо 10.5с на 64-х битной. Мой код для p(2000,1000) показывает 2.4с/8.5с для 32/64 версии.
На моей памяти это первый пример когда 64-х битная заметно медленнее.

-- 18.07.2019, 00:08 --

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

Видимо потому, что вы p(n,n) считаете рекурсивно, а я - встроенной функцией numbpart()

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

[wrest]

Dmitriy40 в сообщении #1405566 писал(а):

На моей памяти это первый пример когда 64-х битная заметно медленнее.

О, вот это неожиданный такой сюрприз

[mihaild]

На моей машине первый ваш вариант работает ~1.7 секунд, второй ~1.3. Наивная реализация на плюсах - 0.2, чуть более интеллектуальная - 0.011 0.007.
https://github.com/mihaild/dxdy/tree/ma ... _partition

[arseniiv]

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

Можете вторую расшифровать, которая pnkd ?

Я не знаю J, но в курсе его похожести на APL и того, что это довольно… экзотические языки по сравнению с распространённым сейчас алголо/паскале/си-подобным синтаксисом. Расшифровал бы если бы ориентировался, но увы…

[mihaild]

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

количество разбиений натурального числа $n$ на не более чем $k$ натуральных слагаемых?

Формула из википедии (и ваш код, и мой) считают число разбиений на слагаемые, не превышающие $k$, а не на не более чем $k$ слагаемых.

[Dmitriy40]

mihaild в сообщении #1405570 писал(а):

На моей машине первый ваш вариант работает ~1.7 секунд, второй ~1.3. Наивная реализация на плюсах - 0.2, чуть более интеллектуальная - 0.011 0.07.
https://github.com/mihaild/dxdy/tree/ma ... _partition

У меня Ваш вариант single_line переписанный обратно на PARI/GP даёт времена 0.56с/2.15с для p(2000,1000) и 13с/46с для p(10000,3000) для 32/64 бит версии PARI/GP. Памяти очевидно всего $O(n)$. Здорово!

[wrest]

mihaild в сообщении #1405572 писал(а):

Формула из википедии (и ваш код, и мой) считают число разбиений на слагаемые, не превышающие $k$, а не на не более чем $k$ слагаемых.

Эти числа равны

[wrest]

mihaild в сообщении #1405570 писал(а):

На моей машине первый ваш вариант работает ~1.7 секунд, второй ~1.3. Наивная реализация на плюсах - 0.2, чуть более интеллектуальная - 0.011 0.07.

Очень круто! Я подозревал что ускорение может быть.
Давайте глянем что вычисляется.
Вот как заполняется cache для $p(20,5)$
[1, 1, 2, 3, 5, 7, 10, 13, 18, 23, 30, 37, 47, 57, 70, 84, 64, 33, 10, 1, 192]
В "single_line" алгоритме массив cache наполнен так: $cache(i)=p(i,k)$ для $ i \le n-k$, а для $n+1>i>n-k, cache(i)=p(i,n-i)$
Ну, первые $k+1$ чисел (в примере выше [1, 1, 2, 3, 5, 7]) это просто количества разбиений (без ограничений на количество слагаемых).
Мне кажется, что тут ещё как-то можно оптимизировать...

-- 18.07.2019, 12:34 --

Да, кстати, мой перевод программы mihaild на pari/gp:

p_mihaild(n,k)=
{
my(M=vector(n),res=0);
M[1]=1;
for(i=1,k,
  for(j=1,n-2*i+1,
     M[i+j]+=M[j]);
  res+=M[n-i+1]);
return(res);
}

Там ещё можно явно дописать, что внутренний цикл срабатывает только на первых $min(k,n/2)$ проходах, дописать за этим ещё один оставшийся цикл по суммированию, но это будут слёзы с точки зрения скорости.
Кстати алгоритм очень хорош тем, что он эффективен как для малых $k$, так и для больших.