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

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

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

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

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

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

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

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

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

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

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

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

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

#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;
}

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

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

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

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

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

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

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

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

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

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

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);
}

?

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

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

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

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

#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;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

разложим 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!

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

На самом деле нам известны 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-ки.

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

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

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

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

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

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

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