Lanczos 계수를 이용하여 Gamma 함수값 구하기

 

 

아래의 소스는

Coefficients for the Lanczos Approximation to the Gamma Function

의 것을 Microsoft Visual C/C++ 명령줄 컴파일러 cl 로 컴파일되도록 수정한 것이다.

 

// Filename: test_gamma_lanczos.c
//
// Compile: cl test_gamma_lanczos.c /EHsc /utf-8
// Execute: test_gamma_lanczos

#include <stdio.h>
	
#define _USE_MATH_DEFINES // for C
#include <math.h>
	
  #define LG_g 5.0      // Lanczos parameter "g"
  #define LG_N 6        // Range of coefficients i=[0..N]
  const double lct[LG_N+1] = {
     1.000000000190015,
    76.18009172947146,
   -86.50532032941677,
    24.01409824083091,
    -1.231739572450155,
     0.1208650973866179e-2,
    -0.5395239384953e-5
  };

  const double ln_sqrt_2_pi = 0.91893853320467274178;
  const double         g_pi = 3.14159265358979323846;

  double lanczos_ln_gamma(double z)
  {
    double sum;
    double base;
    double rv;
    int i;
    
    if (z < 0.5) {
      return log(g_pi / sin(g_pi * z)) - lanczos_ln_gamma(1.0 - z);
    }
    z = z - 1.0;
    base = z + LG_g + 0.5;  // Base of the Lanczos exponential
    sum = 0;

    for(i=LG_N; i>=1; i--) {
      sum += lct[i] / (z + ((double) i));
    }
    sum += lct[0];
    return ((ln_sqrt_2_pi + log(sum)) - base) + log(base)*(z+0.5);
  }

  double lanczos_gamma(double z)
  {
    return(exp(lanczos_ln_gamma(z)));
  }
  
  int main()
  {
      printf("                        Gamma(0.5)               Log(Gamma(0.5))\n");
      printf("                 ---------------------------------------------------\n");
      printf("  Lanczos Gamma: %22.18f     %22.18f\n", lanczos_gamma(0.5), lanczos_ln_gamma(0.5));
      printf("       Built-in: %22.18f     %22.18f\n", exp(lgamma(0.5)), lgamma(0.5));
      printf("   Use sqrt(PI): %22.18f     %22.18f\n", sqrt(M_PI), log(sqrt(M_PI)));

      return 0;	  
  }

 

 

빌트인 Gamma 함수와 비교하여 12자리 까지만 맞는데,

이는 허용오차가 약 sqrt(PI)*N 이기 때문이다.

여기서 N은 Lanczos 계수의 개수이다, (소스에서는 LG_N+1)

Python 인터프러터를 실행하여 잠깐 그 허용오차를 알아본다.

Python 3.7.5 
>>> import math
>>> math.gamma(0.5)
1.7724538509055159
>>> import math
>>> math.sqrt(math.pi)*7
12.407176956338612

N=7 인 경우, 유효수자의 개수는 약 12개이다.