code: public static double Stirling(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 * Special.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 = 2.50662827463100050242E0 * y * w;
          return y;
      }
docs: Gamma function as computed by Stirling's formula.
  
  @param x Value.
  @return Result.