LINUX.ORG.RU

О вычислении биномиальных коэффициентов в целочисленной арифметике


0

5

Не буду напомнинать, что такое биномиальные коэффициенты. Суть проблемы вот в чем: допустим, я хочу вычислять биномиальные коэффициенты эффективно. Я использую такое рекурсивное определение

C(n,k+1) = (n-k)/(k+1) C(n,k) 
следующее из определения этих коэффициентов как C(n,k) = n!/k!/(n-k)!

Будучи реализовано в виде хвостовой рекурсии или итеративно, оно хорошо, работает, быстро и эффективно. Проблема вот в чем: Оно требует _рациональной_арифметики_, т.е. манипулирования дробями, хотя сами биномиальные коэффициенты — целые, что видно из их другого определения

C(n+1,k) = C(n,k) + C(n,k-1)
(«треугольника Паскаля»)

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

Возникает вопрос вынесенный в заголовок: а есть ли способ вычисления биномиальных коэффицентов, более эффективный, чем треугольник Паскаля, но так же использующий только целочисленную арифметику?

но так же использующий только целочисленную арифметику?

C(n,k+1) = ( (n-k)*C(n,k) )/(k+1)

Я тебя не совсем понял. Почему не устраивает такой вариант в использованием целочисленного деления?

pathfinder ★★★★
()

> C(n,k+1) = (n-k)/(k+1) C(n,k)

А если попробовать C(n,k+1) = (n-k) * C(n,k) / (k+1) в целых числах?

const86 ★★★★★
()
Ответ на: комментарий от pathfinder

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

mclaudt
()
Ответ на: комментарий от pathfinder

Ага, хороший вопрос.

Первая проблема то, что числа (n-k)C(n,k) переполняются раньше, чем (n-k)/k C(n,k)

Вторая — величина C(n,k) не известна в момент вычисления, если реализовывать проблему рекурсивно. Я не храню таблиц коэффициентов.

annoynimous ★★★★★
() автор топика

Прошу заметить, что в случае треугольника паскаля для вычисления C(n+1,k) вам совершенно не нужна вся матрица, а только предыдущая ее строка (n-ая). Если вам нужен конкретный биномиальный коэффициент, то обойдетесь малой кровью, храня две строки: текущую и предыдущую.

metar ★★★
()
Ответ на: комментарий от annoynimous

Вот пример кода, который я использую в данный момент:

#include "binomial.h" 

/* Binomial coefficients calculations */

double binomial_helper(unsigned n, unsigned k, double factor);

double
binomial(unsigned n, unsigned k)
{
// Check input
  if ((n < k) || ( k<0))
    return 0.0;

// Using symmetry C^k_n = C^{n-k}_n
  if (k < (n-k))
     return binomial_helper(n, k, 1.0);
  else
     return binomial_helper(n, n-k, 1.0);
}

// Tail recursion by k: C^k_n = C^{k-1}_n * (n-k+1)/k
double
binomial_helper(unsigned n, unsigned k, double factor)
{
  if(k>0)
    return binomial_helper(n, k-1, ((n-k+1.0)/k)*factor);.
  else
    return factor;
}

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от annoynimous

Пардоньте, а почему б не вычислять все значения в двойном цикле последовательно в матрицу N на K?
Если хочется рекурсии, то очень нерационально запускать два рекурсивных вызова, лучше запоминать пары (n, k) в set или, опять же, в матрицу, и проверять перед каждым запуском, не посчитали ли вы уже на самом деле результат раньше.

P.S. man мемоизация
man динамическое программирование

metar ★★★
()
Ответ на: комментарий от mclaudt

>Суть в том что раз удивительным образом выходят всегда целые числа, то можно обойтись без деления вообще. Не знаю, насколько это оправдано вычислительно, но аналитически записать выражение без знаменателя действительно интересно.

Ага, можно обойтись без деления (man «треугольник Паскаля»), только в этом нет ничего интересного.

pathfinder ★★★★
()
Ответ на: комментарий от metar

Я же написал, что это нерациональное использование памяти. Все равно промежуточные значиения нужно хранить и всего лишь ради того, чтобы вычислить 1(!) коэффициент. Скажем, мой код быстро вычисляет C(1000,500) на обычной double арифметике.

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от pathfinder

>>(man «треугольник Паскаля»), только в этом нет ничего интересного.

В треугольнике Паскаля ты можешь записать только рекурсивную формулу (это действительно неинтересно, хотя для понимания полезно), а вот в явном виде без дроби расписать существующее выражение (n)*(n-1)*...*(n-k)/k! это похоже непростая задача теории чисел. То, что результат всегда целый - это удивительный и наинеочевиднейший факт.

mclaudt
()
Ответ на: комментарий от mclaudt

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

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от annoynimous

> Я же написал, что это нерациональное использование памяти.
А вы и первое мое сообщение прочитали? Матрица 2*K - имхо, весьма щадяще.

metar ★★★
()
Ответ на: комментарий от annoynimous

Why not

int binomial_iter(int n,int k,int i,int prev)
{
  return i==k?prev:binomial_iter(n,k,i+1,(n-i)*prev/(i+1));
}

int binomial(int n,int k)
{
  return k<n-k?binomial_iter(n,k,0,1):binomial_iter(n,n-k,0,1);
}

?

Деление есть, но оно целочисленное...

Yareg ★★★
()
Ответ на: комментарий от annoynimous

>>делимость произведениz k подряд идущих чисел на k! это довольно простая теорема теории чисел.

Не путаешь? На k - простая, а вот на k! - я бы с удовольствием заценил его простоту.

mclaudt
()
Ответ на: комментарий от annoynimous

Насчет трудозатрат. Быдлокод за три минуты, но думаю, разберетесь.

#include <cstdio>
#include <algorithm>

const int MAX_K = 100;

int main()	{
	int n, k;
	int cur = 1, prev = 0;
	scanf("%d%d", &n, &k);
	long long c[2][MAX_K+1] = {{0}};
	c[prev][0] = 1;
	for (int i = 1; i <= n; i++)	{
		c[cur][0] = 1;
		for (int j = 1; j <= k; j++)
			c[cur][j] = c[prev][j-1]+c[prev][j];
		std::swap(cur, prev);
	}
	printf("%lld\n", c[prev][k]);
	return 0;
}
metar ★★★
()
Ответ на: комментарий от Yareg

long long слишком быстро переполняется, но там, где не переполняется, работает мгновенно. Рекурсия хвостовая, т.е. выоптимизировывается в цикл.

Yareg ★★★
()

Думаю тему обсудили достаточно, потому не буду советовать динамическое программирование которое явно напрашивается.

Позволю небольшой вопросик на примере этой задачи. Например Haskell - pure functional. Функции без побочных эффектов позволяют вычислять значение функции с определенными значениями параметров только один раз.

Будет ли Хаскелль кешировать вычисленные значения C(n,k)? Если будет, то даже C(n+1,k) = C(n,k) + C(n,k-1) становится эффективным и прекрасным.

vertexua ★★★★★
()
Ответ на: комментарий от vertexua

> Функции без побочных эффектов позволяют вычислять значение функции с определенными значениями параметров только один раз.

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

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от Yareg

Не, в функциональном стиле оно смотрится естественнее.

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от annoynimous

целочисленных

Точнее любых точно сравниваемых. Это так, глобально )

А вот что на практике?

vertexua ★★★★★
()
Ответ на: комментарий от annoynimous

> Да, пока простого доказательства привести не могу.

куда уж проще то? достаточно посчитать

разложим k! на множители:

k! = П p^( [k/p] + [k/p²] + [k/p³] + ... )

очевидно, что [x+y]>=[x]+[y]

поэтому для любого а имеем [(m+n)/pª]>=[m/pª]+[n/pª]

посему (m+n)! делится не просто на m!, а на произведение m!n!

З.Ы. при этом тут может быть много нетривиального

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от vertexua

> Будет ли Хаскелль кешировать вычисленные значения C(n,k)? Если будет, то даже C(n+1,k) = C(n,k) + C(n,k-1) становится эффективным и прекрасным.

А догадается он кэшировать только то, что надо (ну там последние две строки треугольника паскаля) и не кэшировать то, что не надо (весь остальной треугольник, который может быть в 10000 раз больше?)

Что-то сомневаюсь.

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от www_linux_org_ru

[x] — это целая часть от x?

Ок, доказательство принимается.

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от www_linux_org_ru

>>куда уж проще то? достаточно посчитать

Это скатанное с 25-ой страницы Алана Бейкера доказательство ни хрена не является простым.

зы hibow поешь говна.

mclaudt
()
Ответ на: комментарий от vertexua

> Ну хоть best-effort не?

в нормальном языке я бы предложил хинты вроде

@pure @memoize(forget_all(n-2,*)) Int binomial(int n, int k) { return binomial(n,k)+binomial(n,k-1); }

но при этом приобретается зависимость времени вычисления от порядка вычислений

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от mclaudt

> Это скатанное с 25-ой страницы Алана Бейкера доказательство ни хрена не является простым.

Не читал Алана Бейкера.

А что может быть проще, чем тупо посчитать показатели простых чисел в разложении?

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от www_linux_org_ru

Это нормальный вариант, но кажется это похоже на то, что уже есть в хаскелле. Я посмотрел как там сделана мемоизация, ничего не понял, но я не вчитывался подробно. Сейчас немного спешу. Подробнее изучу потом

vertexua ★★★★★
()
Ответ на: комментарий от mclaudt

а формула разложения на множители скатана ЕМНИП у Чебышева... кстати, возможно он первый ее и выписал для какой-то практической цели

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от www_linux_org_ru

Да я знаю про нее. Сам дошел когда-то, когда решал задачу, на сколько нулей оканчивается число 100!

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от vertexua

Сдесь тупо все хранится в мап, и это можно регулировать немного поменяв реализацию. Практически все вручную, никакого автоматизма. Просто используется ленивость. И еще я не знаю что значит !!

vertexua ★★★★★
()

А почему, собственно, при реализации через рациональные числа, числитель и знаменатель будут сильно расти.

У вас есть доказательство? Может, не все так плохо?

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

Сорри, что сумбурно изложил, но более точное изложение потребует некоторого времени.

anonymous
()
Ответ на: комментарий от www_linux_org_ru

>>А что может быть проще, чем тупо посчитать показатели простых чисел в разложении?

Давай различать простоту и очевидность. Теорему можно доказать в пару строк, но она может оставаться неочевидной (т.е. понятной без доказательства вообще).

Так вот я именно про очевидность. Наверняка существует красивая механистическая аналогия, и нужность формализма отпадет вообще.

mclaudt
()
Ответ на: комментарий от anonymous

Вот более точное предложение.

На самом деле нам известны 2 формулы: C(n,k)=n/(n-k)*C(n-1,k) C(n,k)=(n+1-k)/k*C(n,k-1)

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

Как найти правило чередования шагов? Разумеется, мы можем некоторый анализ здесь произвести - и я даже рекомендую это сделать. Это один вариант.

Но второй вариант я вам рекомендую обязательно попробовать - чисто для интереса.

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

Соответствующий анализ метода, думаю, проведете сами.

anonymous
()
Ответ на: комментарий от anonymous

Т.е. предлагается использовать рациональную арифметику - простейший вариант из Кнута, например, думаю вы уже реализовывали.

anonymous
()
Ответ на: комментарий от mclaudt

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

А так не вижу смысла знать Особый Ключ к двери, которая вышибается обычным пинком ноги.

Если чуть подумать, то суть дела именно в [x+y]>=[x]+[y]; k! в определенном смысле худший; например, в отрезок [1,6] попадает не больше одной 5-ки, но в этот же отрезок, смещенный куда-то (например, к [10,15]) может попасть уже две 5-ки.

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от vertexua

с моей точки зрения вариант с мемоизацией более высокоуровневый, чем ручной «разворот» кода, но там обязательно нужно иметь кучу возможностей этой мемоизацией управлять и ее оптимизировать — например, чтобы пары k,n находились не по их хэшу, а в линейном массиве, или даже в кольцевом буфере

www_linux_org_ru ★★★★★
()
Ответ на: комментарий от www_linux_org_ru

Ну да, казалось бы, можно было бы по индукции начать с отрезка [1,6] и доказать что сдвиг его на любой n до [n+1,n+6] не теряет делимости произведения на 6! Но там везде какой-то рандом.

mclaudt
()
Ответ на: комментарий от Yareg

Я реализовал Ваш алгоритм на схеме, в которой есть поддержка bignum. Действительно, все работает:

(define (binomial n k)
;; Helper function to compute C(n,k) via forward recursion 
  (define (binomial-iter n k i prev)
    (if (= k i)
      prev
     (binomial-iter n k (+ i 1) (/ (* (- n i) prev) (+ i 1)))))
;; Use symmetry property C(n,k)=C(n, n-k) 
  (if (< k (-  n k))
    (binomial-iter n k 0 1)
    (binomial-iter n (- n k) 0 1)))

(display (binomial 5000 2500))

output

$ ./binomial.scm | fold

15937186853494383756910257929359190872571294341687297385114218887251837382493431 68628721321987372082216725114025558627431929452538403647953823060045455657769129 87312231391977393396458461430832244929917614864892584730480023643569324440941033 12205443181598892271025936344581075490893586567840537362475844014203104454736819 41570500042422596205820826613570804985152082496060797280773893645848289449988113 02900284896583584182227129255685546148432527106563886001286107377899921899687103 51588991986542785823713096776713060967694323141265378890596537396147632671481841 53554042435532017355139763002046444189510653672101427075975090258426372731594098 22342143017397430156968769327523500761058075924169787048410094516380576648834071 38976421971065422831749537949686497832065589926198994254652723149215649094893728 73935457621535462289487338592584093718521896760237846849202260459597759609074699 35996353111971038349697595257112068688379060221973084397646998338508204919354733 58257987062484987988967617446097612277538575862041300733777299911585208214180922 17638362672853353538377982685094940324887720944387190008030334758815325466470631 38520719622744967680682146446102287294664405588220878554098672156017625258655808 83667267829738243155298628580852338565766233367598650344799146790065148746117755 06177114579705449884133344535068585516842987546616463555923859087650313307613765 39376966719050679737673751870783288866270017279037779325935186879750716333960082 7270303879295319252652644612708041188699645673663207276383716320

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от www_linux_org_ru

>>Если у тебя такая будет, с интересом прочту.

Из тех что лежат на поверхности - это комбинаторная трактовка. Число сочетаний обязано быть целым, а делимость на факториалы видна из перестановки шариков одного цвета.

mclaudt
()
Ответ на: комментарий от mclaudt

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

annoynimous ★★★★★
() автор топика
Ответ на: комментарий от annoynimous

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

mclaudt
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.