Gammac.java

/*
 * $Id: Gammac.java,v 1.13 2008/03/15 00:23:40 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.scalar.DoubleNumber;


/**
 * 完全ガンマ関数を表すクラスです。
 * 
 * <p> Complete Gamma function
 * 
 * @author koga
 * @version $Revision: 1.13 $
 * @see org.mklab.tool.matrix.Gammai
 */
public class Gammac {

  /**
   * <code>x</code>の成分について完全ガンマ関数の値を返します。
   * 
   * @param x データ
   * @return 完全ガンマ関数の値 (value of function)
   */
  public static DoubleMatrix gammac(DoubleMatrix x) {
    DoubleMatrix gam = x.createZero(x.getRowSize(), x.getColumnSize());
    int m = x.getRowSize();
    int n = x.getColumnSize();

    for (int i = 1; i <= m; i++) {
      for (int j = 1; j <= n; j++) {
        gam.setElement(i, j, gammacGammln(x.getElement(i, j)).exp());
      }
    }

    return gam;
  }

//  /**
//   * <code>x</code>の成分について完全ガンマ関数の値を返します。
//   * 
//   * @param x データ
//   * @return 完全ガンマ関数の値 (value of function)
//   */
//  public static DoubleMatrix gammac(DoubleMatrix x) {
//    DoubleMatrix gam = new DoubleMatrix(x.getRowSize(), x.getColumnSize());
//    int m = x.getRowSize();
//    int n = x.getColumnSize();
//
//    for (int i = 1; i <= m; i++) {
//      for (int j = 1; j <= n; j++) {
//        gam.setElement(i, j, Math.exp(gammacGammln(x.getDoubleElement(i, j))));
//      }
//    }
//
//    return gam;
//  }

  /**
   * @param xx データ
   * @return result
   */
  private static DoubleNumber gammacGammln(final DoubleNumber xx) {
    DoubleMatrix cof = new DoubleMatrix(new double[] {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5});

    DoubleNumber y = xx;
    DoubleNumber x = xx;
    DoubleNumber tmp = x.add(5.5);
    tmp = tmp.subtract((x.add(0.5).multiply(tmp.log())));
    DoubleNumber ser = new DoubleNumber(1.000000000190015);
    for (int i = 1; i <= 6; i++) {
      y = y.add(1);
      ser = ser.add(cof.getElement(i).divide(y));
    }

    return ser.multiply(2.5066282746310005).divide(x).log().subtract(tmp);
  }

  /**
   * @param xx データ
   * @return result
   */
  private static double gammacGammln(final double xx) {
    DoubleMatrix cof = new DoubleMatrix(new double[] {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5});

    double y = xx;
    double x = xx;
    double tmp = x + 5.5;
    tmp = tmp - ((x + 0.5) * Math.log(tmp));
    double ser = 1.000000000190015;
    for (int i = 1; i <= 6; i++) {
      y = y + 1;
      ser = ser + cof.getDoubleElement(i) / y;
    }

    return -tmp + Math.log(2.5066282746310005 * ser / x);
  }
}