Аппроксимация логарифма биномиальных коэффициентов для очень больших чисел

Question

В настоящее время я борюсь с вычислением биномиальных коэффициентов для очень больших чисел, скажем, "n выбирать k" с n <10 000 000 и n <k.Это необходимо в контексте вычисления гипергеометрических распределений вероятностей. </p>

До этого момента я пробовал много подходов для обработки больших чисел, которые являются результатом этих вычислений.Однако проблема в том, что мне не нужно вычислять эти биномиальные коэффициенты один раз, а сто тысяч раз.Это означает, что обычные подходы к вычислению факториалов слишком дороги, а стандартные типы данных, такие как long long int, слишком ограничены и не могут содержать эти числа.

Я уже пробовал мультиточные типы данных из библиотеки Boost, но, как я уже говорил, выполнение вычислений много раз приводит к крайне низкой производительности.Я также пробовал многопоточность, используя OpenMP, но выигрыш в производительности все еще был слишком низким.В результате я переключился на вычисление логарифма биномиальных коэффициентов, чтобы числа были небольшими.Хотя это решило проблему больших количеств, это не ускорило процесс.Вот почему я опробовал приближение Стирлинга логарифмических биномиальных коэффициентов.Мое текущее решение выглядит так:

#include <math.h>

long double calc_hgeom(unsigned int k, unsigned int n, unsigned int K, unsigned int N)
{
    long double hprob = std::exp((log_C(K, k) + log_C(N-K, n-k)) - log_C(N, n));
    return hprob;
}

long double log_C(unsigned int u, unsigned int m)
{
    long double C = u * std::log(u) - m * std::log(m) - (u-m) * std::log(u-m)) + 0.5 * (std::log(u) - std::log(m) - std::log(u-m) - std::log(2*M_PI));
    return C;
}

Однако результаты довольно сильно отличаются от фактических значений, до 7%.Отсюда мой вопрос: существует ли эффективный способ вычисления логарифма биномиальных коэффициентов или моё приближение может быть улучшено для повышения точности?

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

Answer 1

Рассмотрим функцию lchoose R ...

> choose(10000, 5000) 
[1] Inf
> lchoose(10000, 5000)
[1] 6926.641

Хранилище исходного кода на языке R является отличным источником идей для подобных проблем.

См. https://github.com/wch/r-source/blob/trunk/src/nmath/choose.c

Хитрость здесь в том, чтобы работать с преобразованными в ln входами, чтобы избежать переполнения.

Обратите внимание, что код находится под лицензией GNU.

Answer 2

Вы должны использовать Формула аппроксимации Стерлинга для n! , который применительно к биномиальным коэффициентам дает:

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

Если этого недостаточно, т. Е. Если у вас все еще есть термины, которые почти сводят друг друга на нет, рассмотрите возможность применения расширения Тейлора к одной из переменных.