2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 16:01 


05/09/16
11469
Общее число разбиений $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$ натуральных слагаемых? Какие еще есть алгоритмы?

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 17:35 


05/09/16
11469
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$

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 18:17 
Заслуженный участник


20/08/14
11068
Россия, Москва
wrest в сообщении #1405468 писал(а):
этот алгоритм требует создания и наполнения массива размером примерно $nk$ элементов, так что при больших $n,k$ просто не хватает памяти. Ну и медленно работает.
Для уменьшения памяти можно сохранять не все варианты комбинаций $n$ и $k$, а только лишь некоторые, например с условием $nk<10^8$ (так легко регулировать потребление памяти), остальные же продолжая вычислять рекурсивно. А ещё лучше - предрасчитать лишь небольшие значения $n$ и $k$, для прочих же класть в хеш-таблицу по ключу $nk$ (а лучше конечно по двум отдельно), тут ускорение оценить сложно, оно зависит от порядка вычислений и размера таблицы.

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 18:55 


05/09/16
11469
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}$

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 23:36 
Заслуженный участник


20/08/14
11068
Россия, Москва
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
Но вообще согласен, алгоритм медленный.

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 23:43 
Заслуженный участник


27/04/09
28128
wrest в сообщении #1405513 писал(а):
Это какая-то... шифрограмма...
Потому что J.

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение17.07.2019, 23:53 


05/09/16
11469
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 ?

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 00:06 
Заслуженный участник


20/08/14
11068
Россия, Москва
О, я понял: Вы запускаете 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()
Это не влияет, у меня в коде пробовал много разных оптимизаций, включая и эту, заметного эффекта нет.

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 00:11 


05/09/16
11469
Dmitriy40 в сообщении #1405566 писал(а):
На моей памяти это первый пример когда 64-х битная заметно медленнее.

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

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 01:10 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 01:51 
Заслуженный участник


27/04/09
28128
wrest в сообщении #1405564 писал(а):
Можете вторую расшифровать, которая pnkd ?
Я не знаю J, но в курсе его похожести на APL и того, что это довольно… экзотические языки по сравнению с распространённым сейчас алголо/паскале/си-подобным синтаксисом. Расшифровал бы если бы ориентировался, но увы…

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 01:58 
Заслуженный участник
Аватара пользователя


16/07/14
8355
Цюрих
wrest в сообщении #1405468 писал(а):
количество разбиений натурального числа $n$ на не более чем $k$ натуральных слагаемых?
Формула из википедии (и ваш код, и мой) считают число разбиений на слагаемые, не превышающие $k$, а не на не более чем $k$ слагаемых.

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 04:33 
Заслуженный участник


20/08/14
11068
Россия, Москва
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)$. Здорово!

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 07:27 


05/09/16
11469
mihaild в сообщении #1405572 писал(а):
Формула из википедии (и ваш код, и мой) считают число разбиений на слагаемые, не превышающие $k$, а не на не более чем $k$ слагаемых.
Эти числа равны 8-)

 Профиль  
                  
 
 Re: Разбиения с ограничениями: кто быстрее?
Сообщение18.07.2019, 12:11 


05/09/16
11469
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$, так и для больших.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 24 ]  На страницу 1, 2  След.

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group