Для вычисления Гамма-функции используется аппроксимация ее логарифма. Сама Гамма, а также Бета - вычисляются через него.
Для аппроксимации гамма-функции на интервале x>0 используется следующая формула (для комплексных z):
Для получения (действительной) гамма-функции на интервале x>0 используется рекуррентная формула Gam(z+1)=z*Gam(z) и вышеприведенная аппроксимация Gam(z+1). Кроме того, можно заметить, что удобнее аппроксимировать логарифм гамма-функции, чем ее саму. Во-первых, при этом потребуется вызов только одной математической функции -- логарифма, а не двух -- экспоненты и степени (последняя все равно использует вызов логарифма), во-вторых, гамма-функция -- быстро растущая для больших x, и аппроксимация ее логарифмом снимает вопросы переполнения.
Для аппроксимации LnGam() -- логарифма гамма-функции -- получается формула:
Сама гамма-функция получается из ее логарифма взятием экспоненты.
Бета-функция тождественно равна B(z,w)=Gam(z)*Gam(w)/Gam(z+w). Ее вычисление сводится к вычислению трех логарифмов гамма-функции и экспоненте от их комбинации.
Ниже приведена программа, реализующая вычисление логарифма гамма-функции, самой гамма-функции и бета-функции.
/* Logarithm of Gamma function, Gamma itself and Beta -- for the value
of argument(s) >0. Double precision.
Copyright (c) Nikitin V.F. 2000
Calculate logarithm of gamma-function for x>0
double GammLn(double x);
Calculate gamma-function for x>0
double Gamma(double x);
Calculate beta-function for x>0, y>0
double Beta(double x,double y);
Note: no domain check is implemented.
*/
/* math.h included for log() and exp() */
#include <math.h>
/* the coefficients table */
#define CN 8
static double cof[CN]={
2.5066282746310005,
1.0000000000190015,
76.18009172947146,
-86.50532032941677,
24.01409824083091,
-1.231739572450155,
0.1208650973866179e-2,
-0.5395239384953e-5,
};
/* logarithm of gamma-function */
double GammLn(double x) {
double y,ser,*co;
int j;
/* calculate the series */
ser=cof[1]; y=x; co=cof+2;
for(j=2;j<CN;j++) {
y+=1.; ser+=(*co)/y; co++;
}
/* and the other parts of the function */
y=x+5.5;
y-=(x+0.5)*log(y);
return(-y+log(cof[0]*ser/x));
}
/* the gamma-function itself */
double Gamma(double x) {
return(exp(GammLn(x)));
}
/* the beta function */
double Beta(double x,double y) {
return(exp(GammLn(x)+GammLn(y)-GammLn(x+y)));
}