감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 C 언어 소스

Lanczos 알고리즘은 Stirling 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp, sin, pow 함수들의 구현에 따라 오차가 더 있을 수 있다. 소스를 보면 (조금이라도 유효수자의 개수가 더 많은 계산을 하기 위해) double 타입이 아니라 long double 타입으로 처리하고 있음을 알 수 있다.

            long double frac(long double x)

            long double gamma2(long double y)

Unix/Linux의 C 컴파일러 gcc 를 사용하면 math 라이브러리에 tgamma 라는 함수가 이미 있다. 이 함수와 (자체 구현된) gamma2 함수에 의한 계산 결과를 비교하기 바란다.

// Filename: testLanczos-02.c
//
// Compile: gcc -o testLanczos-02 -lm testLanczos-02.c
// Execute: ./testLanczos-02
// Output:
//      Gamma(10) = (9)! = 362880
//     tgamma(10) = (9)! = 362880
//
// See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function
// See: http://en.wikipedia.org/wiki/Gamma_function
// See: http://http://jany.tistory.com/53
// See: http://support.realsoftware.com/listarchives/realbasic-nug/2003-08/msg03123.html

#include <stdio.h>
#include <math.h>

long double frac(long double x)
{
  long double result;

  if (x >= 0.0)
      result = fabs(x) - floor(x);
  else
      result = -(fabs(x) - floor(x));
  return result;
}

long double gamma2(long double y)
{
        long double pi = 3.1415926535897932385;
        long double sq2pi = 2.50662827463100050241577;  // sqrt(2Pi)
        long double g = 607./128;   // best resu'ts when 4<=g<=5
        long double t, w, gam, x = y - 1.0;
        int i, cg = 14;
        long double c[16];
        
        // Lanczos approximation for the complex plane
        // calculated using vpa digits(256)
        // the best set of coeffs was selected from
        // a solution space of g=0 to 32 with 1 to 32 terms
        // these coeffs really give superb performance
        // of 15 sig. digits for 0<=real(z)<=171
        // coeffs should sum to about g*g/2+23/24
        
        // http://www.numericana.com/answer/info/godfrey.htm
        
          c[ 1] =        0.99999999999999709182;   // thiBasic arrays start at 1 ?
          c[ 2] =       57.156235665862923517;
          c[ 3] =      -59.597960355475491248;
          c[ 4] =       14.136097974741747174;
          c[ 5] =       -0.49191381609762019978;
          c[ 6] =        0.33994649984811888699/10000;
          c[ 7] =        0.46523628927048575665/10000;
          c[ 8] =       -0.98374475304879564677/10000;
          c[ 9] =        0.15808870322491248884/1000;
          c[10] =       -0.21026444172410488319/1000;
          c[11] =        0.21743961811521264320/1000;
          c[12] =       -0.16431810653676389022/1000;
          c[13] =        0.84418223983852743293/10000;
          c[14] =       -0.26190838401581408670/10000;
          c[15] =        0.36899182659531622704/100000;
          
        // printf(" y = %LG\n", y);
       //  printf(" x = %LG\n", x);

        if ( x < 0.0 ) {
            x = -x;
            if (frac(x) == 0.0) {
                gam = pow(10.0, 4932);
            }
            else {
                t = c[1];
                for (i = 1; i <= cg; i++) {
                    t = t + c[i+1]/(x+i);
                }
                w = x + g + 0.5;
                gam = pow(w, x+0.5)*exp(-w)*sq2pi*t;
                gam = pi*x/(gam*sin(pi*x));
            }
        }
        else {
            t = c[1];
            for (i = 1; i <= cg; i++) {
                t = t + c[i+1]/(x+i);
            }
            w = x + g + 0.5;
            // printf("   w = %LG\n", w);
            gam = pow(w, x+0.5)*exp(-w)*sq2pi*t;
        }
        return gam;
 }

int main()
{
    long double y, x = 1.0;
    while (x > 0) {
        printf("Input x(0 to exit): ");
        scanf("%LF", &x);
        // printf("x = %LG\n", x);
        if (x == 0.0) {
            break;
        }
        y = gamma2(x);
        printf(" Gamma(%LG) = (%LG)! = %30LF\n", x, x - 1, y);
        printf("tgamma(%LG) = (%LG)! = %30lf\n", x, x - 1, tgamma(x));
    }
    // printf("All done. Press any key to finish\n");
    return 0;
}

/*
100! =  93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Execute: ./testLanczos-02
Result:
Input x(0 to exit): 101
 Gamma(101) = (100)! = 93326215443944119652178484825626229006662095017918287040779978859652504703806235528941436474245838774007726779438736471616107541007146050712395913562333118464.000000
tgamma(101) = (100)! = 93326215443947553183793338240612302366006877769275129529769934567734921397913148008797092443225880267680360266584467988186500321788890593731926090238149001216.000000
Input x(0 to exit): 0
*/