Вычисление гамма-функции и бета-функции

Для аппроксимации гамма-функции на интервале x>0 используется формула (для комплексных z) такого вида:

Gam(z+1)=(z+g+0.5)^z+0.5exp(-(z+g+0.5))*Sqrt(2*pi)*[a₀+a₁/(z+1)+a₂/(z+2)+...+a_n/(z+n)+eps]

Она похожа на аппроксимацию Стирлинга, но в ней имеется корректирующая серия. Для значений g=5 и n=6, проверено, что величина погрешности eps не превышает 2*10^-10. Кроме того, погрешность не превышает этой величины на всей правой половине комплексной плоскости: Re z > 0.

Для получения действительной гамма-функции на интервале x>0 используется рекуррентная формула Gam(z+1)=z*Gam(z) и вышеприведенная аппроксимация Gam(z+1). Также можно заметить, что удобнее аппроксимировать логарифм гамма-функции, чем ее саму. Во-первых, при этом потребуется вызов только одной математической функции -- логарифма, а не двух -- экспоненты и степени (последняя все равно использует вызов логарифма), во-вторых, гамма-функция -- быстро растущая для больших x, и аппроксимация ее логарифмом снимает вопросы переполнения.

Для аппроксимации LnGam() -- логарифма гамма-функции -- получается формула:

log(Gam(x))=(x+0.5)*log(x+5.5)-(x+5.5)+log(C₀(C₁+C₂/(x+1)+C₃/(x+2)+...+C₇/(x+8))/x)

Значения коэффициентов C_k являются табличными данными (см. в программе).

Сама гамма-функция получается из ее логарифма взятием экспоненты.

Бета-функция тождественно равна B(z,w)=Gam(z)*Gam(w)/Gam(z+w). Ее вычисление сводится к вычислению трех логарифмов гамма-функции и экспоненте от их комбинации.