Рассчитать биномиальный коэффициент очень надежно

Question

Рассчитать биномиальный коэффициент очень надежно

Каков наилучший метод для вычисления биномиального коэффициента в C ++? Я видел некоторые фрагменты кода, но мне показалось, что это всегда выполнимо только в каком-то определенном регионе. Мне нужен очень, очень, очень надежный расчет. Я попробовал это с гамма-функцией:

unsigned n=N;
unsigned k=2;
number = tgammal(n + 1) / (tgammal(k + 1) * tgammal(n - k + 1));

но он отличается уже при n = 8, k = 2 из 1 (а при n = 30, k = 2 происходит сбой) ..
Мне «только» нужно вычислить примерно как минимум n = 3000 с k = 2.

5

binomial-coefficients c++gamma-function math

Решение

Другие решения

Вы можете использовать асимптотически более эффективную рекуррентную формулу:

constexpr inline size_t binom(size_t n, size_t k) noexcept
{
return
(        k> n  )? 0 :          // out of range
(k==0 || k==n  )? 1 :          // edge
(k==1 || k==n-1)? n :          // first
binom(n - 1, k - 1) * n / k;   // recursive
}

Это будет использовать только Хорошо) операции для расчета C (n, k).

4

Источник

Accepted Answer

это

constexpr inline size_t binom(size_t n, size_t k) noexcept
{
return
(        k> n  )? 0 :          // out of range
(k==0 || k==n  )? 1 :          // edge
(k==1 || k==n-1)? n :          // first
(     k+k < n  )?              // recursive:
(binom(n-1,k-1) * n)/k :       //  path to k=1   is faster
(binom(n-1,k) * n)/(n-k);      //  path to k=n-1 is faster
}

требует O (тт {к, п-к}) Операции надежны и могут быть выполнены во время компиляции.

объяснение Биномиальные коэффициенты определяются как B(n,k)=k!(n-k)!/n!откуда мы видим, что B(n,k)=B(n,n-k), Мы также можем получить рекуррентные отношения n*B(n,k)=(n-k)*B(n-1,k)=k*B(n-1,k-1), Более того, результат тривиален для k=0,1,n,n-1,

За k=2, результат тоже тривиально (n*(n-1))/2,

Конечно, вы можете сделать это и с другими целочисленными типами. Если вам нужно знать биномиальный коэффициент, который превышает наибольший представимый целочисленный тип, вы должны использовать приблизительные методы: double вместо. В этом случае использование гамма-функции является предпочтительным

#include <cmath>
inline double logbinom(double n, double k) noexcept
{
return std::lgamma(n+1)-std::lgamma(n-k+1)-std::lgamma(k+1);
}
inline double binom(double n, double k) noexcept
{
return std::exp(logbinom(n,k));
}

8