Вычисление гамма-функции и бета-функции
Страница 2. Программа вычисления логарифма гамма-функции, гамма-функции и бета-функции

Программа вычисления логарифма гамма-функции, гамма-функции и бета-функции

Нижеприведенная программа реализует вычисление логарифма гамма-функции, самой гамма-функции и бета-функции.

/* Логарифм гамма-фукнции, сама гамма и бета для положительных аргументов двойной точности.

 Вычисление логарифма гамма-функции для x>0
 double GammLn(double x);
 
 Вычисление гамма-функции для x>0
 double Gamma(double x);

 Вычисление бета-функции для x>0, y>0
 double Beta(double x,double y);
*/

/* math.h included for log() and exp() */
#include <math.h>

/* таблица коэффициентов */
#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,
};

/* логарифм гамма-фукнции */
double GammLn(double x) {
 double y,ser,*co;
 int j;
 /* вычисоение последовательностей */
 ser=cof[1]; y=x; co=cof+2;
 for(j=2;j<CN;j++) {
 y+=1.; ser+=(*co)/y; co++;
 }
 /* и других частей функции */
 y=x+5.5;
 y-=(x+0.5)*log(y);
 return(-y+log(cof[0]*ser/x));
}

/* сама гамма-функция */
double Gamma(double x) {
 return(exp(GammLn(x)));
}

/* бета-функция */
double Beta(double x,double y) {
 return(exp(GammLn(x)+GammLn(y)-GammLn(x+y)));
}