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

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

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

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

pathfinder (10.10.2010 19:05:40)

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

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

const86 (10.10.2010 19:06:53)

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

mclaudt (10.10.2010 19:18:43)

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

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

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

annoynimous (10.10.2010 19:21:47)

Да-да, именно так.

annoynimous (10.10.2010 19:22:54)

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

metar (10.10.2010 19:28:45)

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

#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 (10.10.2010 19:29:18)

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

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

metar (10.10.2010 19:35:01)

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

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

pathfinder (10.10.2010 19:36:39)

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

annoynimous (10.10.2010 19:38:05)

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

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

mclaudt (10.10.2010 19:43:31)

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

annoynimous (10.10.2010 19:49:46)

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

metar (10.10.2010 19:49:52)

читал, читал. Нормальная идея. Но... трудозатратно :)

annoynimous (10.10.2010 19:54:13)

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 (10.10.2010 19:58:43)

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

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

mclaudt (10.10.2010 20:03:26)

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

#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 (10.10.2010 20:05:21)

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

Yareg (10.10.2010 20:20:05)

Ага, т.е. Вы реализовали восходящую рекурсию. Понятно. Что ж, неплохо, думаю, пойдет

annoynimous (10.10.2010 20:42:42)

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

annoynimous (10.10.2010 20:43:37)

Можно сразу в цикл переделать...

Yareg (10.10.2010 21:02:00)

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

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

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

vertexua (10.10.2010 21:10:23)

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

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

annoynimous (10.10.2010 21:13:21)

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

annoynimous (10.10.2010 21:13:55)

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

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

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

vertexua (10.10.2010 21:17:07)

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

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

разложим 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 (10.10.2010 22:22:53)

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

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

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

www_linux_org_ru (10.10.2010 22:25:16)

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

vertexua (10.10.2010 22:26:17)

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

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

annoynimous (10.10.2010 22:35:41)

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

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

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

mclaudt (10.10.2010 22:39:43)

> Ну хоть 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 (10.10.2010 22:40:40)

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

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

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

www_linux_org_ru (10.10.2010 22:43:02)

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

vertexua (10.10.2010 22:45:15)

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

www_linux_org_ru (10.10.2010 22:45:57)

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

annoynimous (10.10.2010 22:47:31)

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

vertexua (10.10.2010 22:48:05)

кстати, что за книжка Бейкера?

www_linux_org_ru (10.10.2010 22:48:25)

Эта. Был перевод ещё. На этой доказательство на восьмой.

mclaudt (10.10.2010 22:57:09)

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

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

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

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

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

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

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

mclaudt (11.10.2010 0:59:54)

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

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

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

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

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

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

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

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

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

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

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

www_linux_org_ru (11.10.2010 1:44:44)

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

www_linux_org_ru (11.10.2010 1:49:30)

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

mclaudt (11.10.2010 1:51:12)

Я реализовал Ваш алгоритм на схеме, в которой есть поддержка 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))

annoynimous (11.10.2010 4:01:35)

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

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

mclaudt (11.10.2010 4:07:15)

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

annoynimous (11.10.2010 4:27:58)

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

mclaudt (11.10.2010 4:53:51)