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

Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp, sin, pow 함수들의 구현에 따라 오차가 더 있을 수 있다. C 언어에는 long double 타입이 있어서 좀 더 정밀한 계산을 할 수 있지만. Java 언어에는 그렇지 못하여 그냥 double 타입을 사용하였다. 비교를 위해 아파치의 commons의 Math 라이브러라에서 지원하는 Gamma.logGamma 함수를 사용한 결과도 함께 출력하도록 하였다. (commons 라이브러리는 http://commons.apache.org/math/download_math.cgi 에서 다운로드할 수 있다.

commons 2.2와 common 3.0은 차이가 있으며 서로 호환되지 않는다.

우선 import 구문이 서로 다르다.

        import org.apache.commons.math.special.*;     // commons 2.2의 경우

        import org.apache.commons.math3.special.*;     // commons 3.0의 경우

common 2.2의 경우에는 org.apache.commons.math.special.Gamma 클래스에 gamma 라는함수가 있어서 이를 직접 호출하면 되지만, commons 3.0의 경우에는 org.apache.commons.math3.special.Gamma 클래스에 gamma 라는 함수가 없으므로, logGamma 함수와 exp 함수의 합성함수로 계산하였다.

  exp(Gamma.logGamma(double)) 함수와 (자체 구현된) gamma2(double) 함수에 의한 계산 결과를 비교하기 바란다.

/**
 * Filename: TestLanczos_02.java
 *
 * Compile: javac -d . -cp .:commons-math3-3.0.jar TestLanczos_02.java
 * Execute: java -classpath .:commons-math3-3.0.jar TestLanczos_02
 *
 * Note: You can download commons-math3-3.0-bin.tar.gz or commons-math3-3.0-bin.zip
 *         at http://commons.apache.org/math/download_math.cgi
 *
 * ----------------------------------------------------------------------------
 * 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
 * ----------------------------------------------------------------------------
 * Date: 2012. 12. 12.
 * Copyright (C) Pil Ho, Kim  (pkim (AT) scripts.pe.kr)
 */

import java.io.*;
import org.apache.commons.math3.special.Gamma;

public class TestLanczos_02 {

     private static double frac(double x) {
         double result;
         if (x >= 0.0)
             result = Math.abs(x) - Math.floor(Math.abs(x));
         else
             result = -(Math.abs(x) - Math.floor(Math.abs(x)));
        return result;
    }

    public static double gamma2(double y)  {
        double pi = 3.1415926535897932385;
        double sq2pi = 2.50662827463100050241577;  // sqrt(2Pi)
        double g = 607./128;   // best resu'ts when 4<=g<=5
        double t, w, gam, x = y - 1.0;
        int i, cg = 14;
        double[] c = new double[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;
          
        if ( x < 0.0 ) {
            x = -x;
            if (frac(x) == 0.0) {
                gam = Math.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 = Math.pow(w, x+0.5)*Math.exp(-w)*sq2pi*t;
                gam = pi*x/(gam*Math.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;
            gam = Math.pow(w, x+0.5)*Math.exp(-w)*sq2pi*t;
        }
        return gam;
    }

    public static void main(String[] args) {
        double y, x = 1.0;

        while (x > 0.0) {
            BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
            String line = null;

            try {
                System.out.printf("Input x(0 to exit): ");
                line = br.readLine();
                x = Double.parseDouble((line));
            }
            catch (IOException ioe) {
                System.err.println("IO error trying to read your input number!");
                System.err.println("Try to input again.");
            }

            if (x == 0.0) {
                break;
            }
            y = gamma2(x);
            System.out.printf(" Gamma(%f) = (%f)! = %30f\n", x, x - 1, y);

            System.out.printf(" Gamma.logGamma(%f) = log((%f)!) = %f\n", x, x - 1, Gamma.logGamma(x));
            System.out.printf(" Math.exp(Gamm.logGamma(%f)) = (%f)! = %30f\n", x, x - 1, Math.exp(Gamma.logGamma(x)));
        }
    }
}

/*
100! =  93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Execute: java -classpath .:commons-math3-3.0.jar TestLanczos_02
Result:
Input x(0 to exit): 101
 Gamma(101.000000) = (100.000000)! = 93326215443944150000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
 Gamma.logGamma(101.000000) = log((100.000000)!) = 363.739376
 Math.exp(Gamm.logGamma(101.000000)) = (100.000000)! = 93326215443942250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
Input x(0 to exit): 0
*/