Gammai.java

/*
 * $Id: Gammai.java,v 1.24 2008/07/16 15:40:00 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

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


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

  /**
   * <code>a</code> の成分について、不完全ガンマ関数の値( <code>x</code>まで 積分した値)を返します。
   * 
   * <p>不完全ガンマ関数は、 <code>0</code>から<code>x</code>まで積分した関数 <pre><code>
   * 
   * gammai(a,x) = 1/gammac(a) integral(0,x) exp(-t)*t&circ;(a-1)
   * 
   * </code></pre> です。ただし、<code>gammac()</code>は完全ガンマ関数です。
   * 
   * @param a データ
   * @param x 積分区間
   * @return 不完全ガンマ関数の値 (incomplete gamma function)
   */
  public static DoubleMatrix gammai(final DoubleMatrix a, final DoubleNumber x) {
    if (x.isLessThan(0) || a.compareElementWise(".<=", 0).anyTrue()) { //$NON-NLS-1$
      throw new IllegalArgumentException(Messages.getString("Gammai.0")); //$NON-NLS-1$
    }

    DoubleMatrix gam = a.createZero(a.getRowSize(), a.getColumnSize());
    int m = a.getRowSize();
    int n = a.getColumnSize();

    for (int i = 1; i <= m; i++) {
      for (int j = 1; j <= n; j++) {
        if (a.getElement(i, j).add(1).isGreaterThan(x)) {
          List<Double> tmp = gammaiGser(a.getElement(i, j), x);
          double gamser = tmp.get(0).doubleValue();
          gam.setElement(i, j, gamser);
        } else {
          List<Double> tmp = gammaiGcf(a.getElement(i, j), x);
          double gammcf = tmp.get(0).doubleValue();
          gam.setElement(i, j, 1 - gammcf);
        }
      }
    }

    return gam;
  }

  /**
   * @param a データ
   * @param x 積分区間
   * @return result
   */
  private static List<Double> gammaiGser(final DoubleNumber a, final DoubleNumber x) {
    int itmax = 100;

    DoubleMatrix gg = Gammac.gammac(a.createGrid(new DoubleNumber[] {a}));
    DoubleNumber gln = gg.getElement(1, 1).log();

    if (x.isLessThanOrEquals(0)) {
      if (x.isLessThan(0)) {
        throw new IllegalArgumentException(Messages.getString("Gammai.1")); //$NON-NLS-1$
      }
      DoubleNumber gamser = a.createZero();
      return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gamser.doubleValue()), Double.valueOf(gln.doubleValue())}));
    }

    DoubleNumber ap = a;
    DoubleNumber del = a.inverse();
    DoubleNumber sum_ = a.inverse();
    for (int n = 1; n <= itmax; n++) {
      ap = ap.add(1);
      del = del.multiply(x.divide(ap));
      sum_ = sum_.add(del);
      if (del.abs().isLessThan(sum_.abs().multiply(a.getMachineEpsilon()))) {
        DoubleNumber gamser = sum_.multiply(x.unaryMinus().add(a.multiply(x.log())).subtract(gln).exp());
        return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gamser.doubleValue()), Double.valueOf(gln.doubleValue())}));
      }
    }
    throw new IllegalArgumentException(Messages.getString("Gammai.2")); //$NON-NLS-1$
  }

  /**
   * @param a データ
   * @param x 積分区間
   * @return result
   */
  private static List<Double> gammaiGcf(final DoubleNumber a, final DoubleNumber x) {
    int itmax = 100;
    DoubleNumber gold = a.createZero();
    DoubleNumber fac = a.createUnit();
    DoubleNumber b1 = a.createUnit();
    DoubleNumber b0 = a.createZero();
    DoubleNumber a0 = a.createUnit();

    DoubleMatrix gg = Gammac.gammac(a.createGrid(new DoubleNumber[] {a}));
    DoubleNumber gln = gg.getElement(1, 1).log();

    DoubleNumber a1 = x;
    for (int n = 1; n <= itmax; n++) {
      DoubleNumber ana = a.unaryMinus().add(n);
      a0 = a1.add(a0.multiply(ana)).multiply(fac);
      b0 = b1.add(b0.multiply(ana)).multiply(fac);
      DoubleNumber anf = fac.multiply(n);
      a1 = x.multiply(a0).add(anf.multiply(a1));
      b1 = x.multiply(b0).add(anf.multiply(b1));
      if (a1.isZero() == false) {
        fac = a1.inverse();
        DoubleNumber g = b1.multiply(fac);
        if (g.subtract(gold).divide(g).abs().isLessThan(a.getMachineEpsilon())) {
          DoubleNumber gammcf = x.unaryMinus().add(a.multiply(x.log())).subtract(gln).exp().multiply(g);
          return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gammcf.doubleValue()), Double.valueOf(gln.doubleValue())}));
        }
        gold = g;
      }
    }
    throw new IllegalArgumentException(Messages.getString("Gammai.3")); //$NON-NLS-1$
  }

  /**
   * <code>a</code> の成分について、不完全ガンマ関数の値( <code>x</code>まで 積分した値)を返します。
   * 
   * <p>不完全ガンマ関数は、 <code>0</code>から<code>x</code>まで積分した関数 <pre><code>
   * 
   * gammai(a,x) = 1/gammac(a) integral(0,x) exp(-t)*t&circ;(a-1)
   * 
   * </code></pre> です。ただし、<code>gammac()</code>は完全ガンマ関数です。
   * 
   * @param a データ
   * @param x 積分区間
   * @return 不完全ガンマ関数の値 (incomplete gamma function)
   */
  public static DoubleMatrix gammai(final DoubleMatrix a, final double x) {
    if (x < 0 || a.compareElementWise(".<=", 0).anyTrue()) { //$NON-NLS-1$
      throw new IllegalArgumentException(Messages.getString("Gammai.4")); //$NON-NLS-1$
    }

    DoubleMatrix gam = new DoubleMatrix(a.getRowSize(), a.getColumnSize());
    int m = a.getRowSize();
    int n = a.getColumnSize();

    for (int i = 1; i <= m; i++) {
      for (int j = 1; j <= n; j++) {
        if (x < a.getDoubleElement(i, j) + 1) {
          List<Double> tmp = gammaiGser(a.getDoubleElement(i, j), x);
          double gamser = tmp.get(0).doubleValue();
          gam.setElement(i, j, gamser);
        } else {
          List<Double> tmp = gammaiGcf(a.getDoubleElement(i, j), x);
          double gammcf = tmp.get(0).doubleValue();
          gam.setElement(i, j, 1 - gammcf);
        }
      }
    }

    return gam;
  }

  /**
   * @param a データ
   * @param x 積分区間
   * @return result
   */
  private static List<Double> gammaiGser(final double a, final double x) {
    int itmax = 100;

    DoubleMatrix gg = Gammac.gammac(new DoubleMatrix(new double[] {a}));
    double gln = Math.log(gg.getDoubleElement(1, 1));

    if (x <= 0.0) {
      if (x < 0.0) {
        throw new IllegalArgumentException(Messages.getString("Gammai.5")); //$NON-NLS-1$
      }
      double gamser = 0;
      return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gamser), Double.valueOf(gln)}));
    }

    double ap = a;
    double del = 1 / a;
    double sum_ = 1 / a;
    for (int n = 1; n <= itmax; n++) {
      ap = ap + 1;
      del = del * (x / ap);
      sum_ = sum_ + del;
      if (Math.abs(del) < Math.abs(sum_) * DoubleNumberUtil.EPS) {
        double gamser = sum_ * Math.exp(-x + a * Math.log(x) - gln);
        return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gamser), Double.valueOf(gln)}));
      }
    }
    throw new IllegalArgumentException(Messages.getString("Gammai.6")); //$NON-NLS-1$
  }

  /**
   * @param a データ
   * @param x 積分区間
   * @return result
   */
  private static List<Double> gammaiGcf(final double a, final double x) {
    int itmax = 100;
    double gold = 0;
    double fac = 1;
    double b1 = 1;
    double b0 = 0;
    double a0 = 1;

    DoubleMatrix gg = Gammac.gammac(new DoubleMatrix(new double[] {a}));
    double gln = Math.log(gg.getDoubleElement(1, 1));

    double a1 = x;
    for (int n = 1; n <= itmax; n++) {
      double ana = n - a;
      a0 = (a1 + a0 * ana) * fac;
      b0 = (b1 + b0 * ana) * fac;
      double anf = n * fac;
      a1 = x * a0 + anf * a1;
      b1 = x * b0 + anf * b1;
      if (a1 != 0) {
        fac = 1 / a1;
        double g = b1 * fac;
        if (Math.abs((g - gold) / g) < DoubleNumberUtil.EPS) {
          double gammcf = Math.exp(-x + a * Math.log(x) - gln) * g;
          return new ArrayList<>(Arrays.asList(new Double[] {Double.valueOf(gammcf), Double.valueOf(gln)}));
        }
        gold = g;
      }
    }
    throw new IllegalArgumentException(Messages.getString("Gammai.7")); //$NON-NLS-1$
  }

}