public static double stirlingFormula(double x){
  double STIR[]={7.87311395793093628397E-4,-2.29549961613378126380E-4,-2.68132617805781232825E-3,3.47222221605458667310E-3,8.33333333333482257126E-2};
  double MAXSTIR=143.01608;
  double w=1.0 / x;
  double y=Math.exp(x);
  w=1.0 + w * polevl(w,STIR,4);
  if (x > MAXSTIR) {
    double v=Math.pow(x,0.5 * x - 0.25);
    y=v * (v / y);
  }
 else {
    y=Math.pow(x,x - 0.5) / y;
  }
  y=SQTPI * y * w;
  return y;
}

public static double stirlingFormula ( double x ) { double STIR [ ] = { 7.87311395793093628397E-4 , - 2.29549961613378126380E-4 , - 2.68132617805781232825E-3 , 3.47222221605458667310E-3 , 8.33333333333482257126E-2 } ; double MAXSTIR = 143.01608 ; double w = 1.0 / x ; double y = Math . exp ( x ) ; w = 1.0 + w * polevl ( w , STIR , 4 ) ; if ( x > MAXSTIR ) { double v = Math . pow ( x , 0.5 * x - 0.25 ) ; y = v * ( v / y ) ; } else { y = Math . pow ( x , x - 0.5 ) / y ; } y = SQTPI * y * w ; return y ; }
returns the gamma function computed by stirling ' s formula . the polynomial stir is valid for 33 < = x < = 172 .
