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

Written on .

ОГЛАВЛЕНИЕ

Для вычисления гамма-функции используется аппроксимация логарифма гамма-функции. Сама же гамма, а также бета - вычисляются через него.

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

Gam(z+1)=(z+g+0.5)z+0.5exp(-(z+g+0.5))*Sqrt(2*pi)*[a0+a1/(z+1)+a2/(z+2)+...+an/(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(C0(C1+C2/(x+1)+C3/(x+2)+...+C7/(x+8))/x)

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

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

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


 

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

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

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

   Вычисление логарифма гамма-функции для 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)));
}
форум firemonkey