// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.

#ifndef __S12_GAMMA_HPP
#define __S12_GAMMA_HPP

#include <cmath>
#include <sstream>
#include <iostream>
#include <stdexcept>
#include <cfloat>

namespace SNL {
namespace S12 {

// Note that the functions Gamma and LogGamma are mutually dependent.
double LogGamma(double);
double Gamma(double);


// Note that the functions Gamma and LogGamma are mutually dependent.

inline double Gamma
(
  double x    // We require x > 0
)
{
  if (x <= 0.0)
  {
    std::stringstream os;
    os << "Invalid input argument " << x <<  ". Argument must be positive.";
    throw std::invalid_argument(os.str());
  }

  // Split the function domain into three intervals:
  // (0, 0.001), [0.001, 12), and (12, infinity)

  ///////////////////////////////////////////////////////////////////////////
  // First interval: (0, 0.001)
  //
  // For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
  // So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
  // The relative error over this interval is less than 6e-7.

  const double gamma = 0.577215664901532860606512090; // Euler's gamma constant

  if (x < 0.001)
    return 1.0 / (x * (1.0 + gamma * x));

  ///////////////////////////////////////////////////////////////////////////
  // Second interval: [0.001, 12)

  if (x < 12.0)
  {
    // The algorithm directly approximates gamma over (1,2) and uses
    // reduction identities to reduce other arguments to this interval.

    double y = x;
    int n = 0;
    bool arg_was_less_than_one = (y < 1.0);

    // Add or subtract integers as necessary to bring y into (1,2)
    // Will correct for this below
    if (arg_was_less_than_one)
    {
      y += 1.0;
    }
    else
    {
      n = static_cast<int>(floor(y)) - 1;   // will use n later
      y -= n;
    }

    // numerator coefficients for approximation over the interval (1,2)
    static const double p[] =
    {
      -1.71618513886549492533811E+0,
      2.47656508055759199108314E+1,
      -3.79804256470945635097577E+2,
      6.29331155312818442661052E+2,
      8.66966202790413211295064E+2,
      -3.14512729688483675254357E+4,
      -3.61444134186911729807069E+4,
      6.64561438202405440627855E+4
    };

    // denominator coefficients for approximation over the interval (1,2)
    static const double q[] =
    {
      -3.08402300119738975254353E+1,
      3.15350626979604161529144E+2,
      -1.01515636749021914166146E+3,
      -3.10777167157231109440444E+3,
      2.25381184209801510330112E+4,
      4.75584627752788110767815E+3,
      -1.34659959864969306392456E+5,
      -1.15132259675553483497211E+5
    };

    double num = 0.0;
    double den = 1.0;
    int i;

    double z = y - 1;

    for (i = 0; i < 8; i++)
    {
      num = (num + p[i]) * z;
      den = den * z + q[i];
    }

    double result = num / den + 1.0;

    // Apply correction if argument was not initially in (1,2)
    if (arg_was_less_than_one)
    {
      // Use identity gamma(z) = gamma(z+1)/z
      // The variable "result" now holds gamma of the original y + 1
      // Thus we use y-1 to get back the orginal y.
      result /= (y - 1.0);
    }
    else
    {
      // Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
      for (i = 0; i < n; i++)
        result *= y++;
    }

    return result;
  }

  ///////////////////////////////////////////////////////////////////////////
  // Third interval: [12, infinity)

  if (x > 171.624)
  {
    // Correct answer too large to display. Force +infinity.
    double temp = DBL_MAX;
    return temp * 2.0;
  }

  return exp(LogGamma(x));
}

inline double LogGamma
(
  double x    // x must be positive
)
{
  if (x <= 0.0)
  {
    std::stringstream os;
    os << "Invalid input argument " << x <<  ". Argument must be positive.";
    throw std::invalid_argument(os.str());
  }

  if (x < 12.0)
  {
    return log(fabs(Gamma(x)));
  }

  // Abramowitz and Stegun 6.1.41
  // Asymptotic series should be good to at least 11 or 12 figures
  // For error analysis, see Whittiker and Watson
  // A Course in Modern Analysis (1927), page 252

  static const double c[8] =
  {
    1.0 / 12.0,
    -1.0 / 360.0,
    1.0 / 1260.0,
    -1.0 / 1680.0,
    1.0 / 1188.0,
    -691.0 / 360360.0,
    1.0 / 156.0,
    -3617.0 / 122400.0
  };
  double z = 1.0 / (x * x);
  double sum = c[7];

  for (int i = 6; i >= 0; i--)
  {
    sum *= z;
    sum += c[i];
  }

  double series = sum / x;

  static const double halfLogTwoPi = 0.91893853320467274178032973640562;
  double logGamma = (x - 0.5) * log(x) - x + halfLogTwoPi + series;
  return logGamma;
}

} // namespace S12
} // namespace SNL

#endif
