Рассчитать биномиальный коэффициент очень надежно
Каков наилучший метод для вычисления биномиального коэффициента в 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.
Решение
это
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));
}
Другие решения
Вы можете использовать асимптотически более эффективную рекуррентную формулу:
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).
detector