public static double lnGamma(double x){
  double p, q, w, z;
  double A[]={8.11614167470508450300E-4,-5.95061904284301438324E-4,7.93650340457716943945E-4,-2.77777777730099687205E-3,8.33333333333331927722E-2};
  double B[]={-1.37825152569120859100E3,-3.88016315134637840924E4,-3.31612992738871184744E5,-1.16237097492762307383E6,-1.72173700820839662146E6,-8.53555664245765465627E5};
  double C[]={-3.51815701436523470549E2,-1.70642106651881159223E4,-2.20528590553854454839E5,-1.13933444367982507207E6,-2.53252307177582951285E6,-2.01889141433532773231E6};
  if (x < -34.0) {
    q=-x;
    w=lnGamma(q);
    p=Math.floor(q);
    if (p == q) {
      throw new ArithmeticException("lnGamma: Overflow");
    }
    z=q - p;
    if (z > 0.5) {
      p+=1.0;
      z=p - q;
    }
    z=q * Math.sin(Math.PI * z);
    if (z == 0.0) {
      throw new ArithmeticException("lnGamma: Overflow");
    }
    z=LOGPI - Math.log(z) - w;
    return z;
  }
  if (x < 13.0) {
    z=1.0;
    while (x >= 3.0) {
      x-=1.0;
      z*=x;
    }
    while (x < 2.0) {
      if (x == 0.0) {
        throw new ArithmeticException("lnGamma: Overflow");
      }
      z/=x;
      x+=1.0;
    }
    if (z < 0.0) {
      z=-z;
    }
    if (x == 2.0) {
      return Math.log(z);
    }
    x-=2.0;
    p=x * polevl(x,B,5) / p1evl(x,C,6);
    return (Math.log(z) + p);
  }
  if (x > 2.556348e305) {
    throw new ArithmeticException("lnGamma: Overflow");
  }
  q=(x - 0.5) * Math.log(x) - x + 0.91893853320467274178;
  if (x > 1.0e8) {
    return (q);
  }
  p=1.0 / (x * x);
  if (x >= 1000.0) {
    q+=((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x;
  }
 else {
    q+=polevl(p,A,4) / x;
  }
  return q;
}

public static double lnGamma ( double x ) { double p , q , w , z ; double A [ ] = { 8.11614167470508450300E-4 , - 5.95061904284301438324E-4 , 7.93650340457716943945E-4 , - 2.77777777730099687205E-3 , 8.33333333333331927722E-2 } ; double B [ ] = { - 1.37825152569120859100E3 , - 3.88016315134637840924E4 , - 3.31612992738871184744E5 , - 1.16237097492762307383E6 , - 1.72173700820839662146E6 , - 8.53555664245765465627E5 } ; double C [ ] = { - 3.51815701436523470549E2 , - 1.70642106651881159223E4 , - 2.20528590553854454839E5 , - 1.13933444367982507207E6 , - 2.53252307177582951285E6 , - 2.01889141433532773231E6 } ; if ( x < - 34.0 ) { q = - x ; w = lnGamma ( q ) ; p = Math . floor ( q ) ; if ( p == q ) { throw new ArithmeticException ( "lnGamma: Overflow" ) ; } z = q - p ; if ( z > 0.5 ) { p += 1.0 ; z = p - q ; } z = q * Math . sin ( Math . PI * z ) ; if ( z == 0.0 ) { throw new ArithmeticException ( "lnGamma: Overflow" ) ; } z = LOGPI - Math . log ( z ) - w ; return z ; } if ( x < 13.0 ) { z = 1.0 ; while ( x >= 3.0 ) { x -= 1.0 ; z *= x ; } while ( x < 2.0 ) { if ( x == 0.0 ) { throw new ArithmeticException ( "lnGamma: Overflow" ) ; } z /= x ; x += 1.0 ; } if ( z < 0.0 ) { z = - z ; } if ( x == 2.0 ) { return Math . log ( z ) ; } x -= 2.0 ; p = x * polevl ( x , B , 5 ) / p1evl ( x , C , 6 ) ; return ( Math . log ( z ) + p ) ; } if ( x > 2.556348e305 ) { throw new ArithmeticException ( "lnGamma: Overflow" ) ; } q = ( x - 0.5 ) * Math . log ( x ) - x + 0.91893853320467274178 ; if ( x > 1.0e8 ) { return ( q ) ; } p = 1.0 / ( x * x ) ; if ( x >= 1000.0 ) { q += ( ( 7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3 ) * p + 0.0833333333333333333333 ) / x ; } else { q += polevl ( p , A , 4 ) / x ; } return q ; }
returns natural logarithm of gamma function .
