//+------------------------------------------------------------------+
//|                                                        gamma.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "config.mqh"
#include "error.mqh"
#include "const.mqh"
#include "polevl.mqh"
#include "trig.mqh"
namespace special
{
namespace cephes
{
namespace detail
{
 double gamma_P[] = {1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2,
                              4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1,
                              9.99999999999999996796E-1
                             };

 double gamma_Q[] = {-2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3,
                              1.18139785222060435552E-2,  3.58236398605498653373E-2, -2.34591795718243348568E-1,
                              7.14304917030273074085E-2,  1.00000000000000000320E0
                             };

/* Stirling's formula for the Gamma function */
 double gamma_STIR[5] =
  {
   7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3,
   3.47222221605458667310E-3, 8.33333333333482257126E-2,
  };

 double MAXSTIR = 143.01608;

/* Gamma function computed by Stirling's formula.
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
 inline double stirf(double x)
  {
   double y, w, v;

   if(x >= MAXGAM)
     {
      return (infinity());
     }
   w = 1.0 / x;
   w = 1.0 + w * special::cephes::polevl(w, gamma_STIR, 4);
   y = exp(x);
   if(x > MAXSTIR)    /* Avoid overflow in pow() */
     {
      v = pow(x, 0.5 * x - 0.25);
      y = v * (v / y);
     }
   else
     {
      y = pow(x, x - 0.5) / y;
     }
   y = SQRTPI * y * w;
   return (y);
  }
} // namespace detail

 inline double Gamma(double x)
  {
   double p, q, z;
   int i;
   int sgngam = 1;

   if(!isfinite(x))
     {
      return x;
     }
   q = abs(x);

   if(q > 33.0)
     {
      if(x < 0.0)
        {
         p = floor(q);
         if(p == q)
           {
            set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
            return (infinity());
           }
         i = (int)p;
         if((i & 1) == 0)
           {
            sgngam = -1;
           }
         z = q - p;
         if(z > 0.5)
           {
            p += 1.0;
            z = q - p;
           }
         z = q * sinpi(z);
         if(z == 0.0)
           {
            return (sgngam * infinity());
           }
         z = abs(z);
         z = M_PI / (z * detail::stirf(q));
        }
      else
        {
         z = detail::stirf(x);
        }
      return (sgngam * z);
     }

   z = 1.0;
   while(x >= 3.0)
     {
      x -= 1.0;
      z *= x;
     }

   while(x < 0.0)
     {
      if(x > -1.E-9)
        {
         if(x == 0.0)
     {
      set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
      return (infinity());
     }
   else
      return (z / ((1.0 + 0.5772156649015329 * x) * x));
        }
      z /= x;
      x += 1.0;
     }

   while(x < 2.0)
     {
      if(x < 1.e-9)
        {
         if(x == 0.0)
     {
      set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
      return (infinity());
     }
   else
      return (z / ((1.0 + 0.5772156649015329 * x) * x));
        }
      z /= x;
      x += 1.0;
     }

   if(x == 2.0)
     {
      return (z);
     }

   x -= 2.0;
   p = polevl(x, detail::gamma_P, 6);
   q = polevl(x, detail::gamma_Q, 7);
   return (z * p / q);

   if(x == 0.0)
     {
      set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
      return (infinity());
     }
   else
      return (z / ((1.0 + 0.5772156649015329 * x) * x));
  }

namespace detail
{
/* A[]: Stirling's formula expansion of log Gamma
 * B[], C[]: log Gamma function between 2 and 3
 */
 double gamma_A[] = {8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
                              -2.77777777730099687205E-3, 8.33333333333331927722E-2
                             };

 double gamma_B[] = {-1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
                              -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5
                             };

 double gamma_C[] =
  {
   /* 1.00000000000000000000E0, */
   -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
      -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6
     };

/* log( sqrt( 2*pi ) ) */
 double LS2PI = 0.91893853320467274178;

 double MAXLGM = 2.556348e305;

/* Disable optimizations for this function on 32 bit systems when compiling with GCC.
 * We've found that enabling optimizations can result in degraded precision
 * for this asymptotic approximation in that case. */
 inline double lgam_large_x(double x)
  {
   double q = (x - 0.5) * log(x) - x + LS2PI;
   if(x > 1.0e8)
     {
      return (q);
     }
   double p = 1.0 / (x * x);
   p = ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x;
   return q + p;
  }

 inline double lgam_sgn(double x, int &sign)
  {
   double p, q, u, w, z;
   int i;

   sign = 1;

   if(!isfinite(x))
     {
      return x;
     }

   if(x < -34.0)
     {
      q = -x;
      w = lgam_sgn(q, sign);
      p = floor(q);
      if(p == q)
        {
         set_error("lgam", SF_ERROR_SINGULAR, NULL);
         return (infinity());
        }
      i = (int)p;
      if((i & 1) == 0)
        {
         sign = -1;
        }
      else
        {
         sign = 1;
        }
      z = q - p;
      if(z > 0.5)
        {
         p += 1.0;
         z = p - q;
        }
      z = q * sinpi(z);
      if(z == 0.0)
        {
         set_error("lgam", SF_ERROR_SINGULAR, NULL);
         return (infinity());
        }
      /*     z = log(M_PI) - log( z ) - w; */
      z = LOGPI - log(z) - w;
      return (z);
     }

   if(x < 13.0)
     {
      z = 1.0;
      p = 0.0;
      u = x;
      while(u >= 3.0)
        {
         p -= 1.0;
         u = x + p;
         z *= u;
        }
      while(u < 2.0)
        {
         if(u == 0.0)
           {
            set_error("lgam", SF_ERROR_SINGULAR, NULL);
         return (infinity());
           }
         z /= u;
         p += 1.0;
         u = x + p;
        }
      if(z < 0.0)
        {
         sign = -1;
         z = -z;
        }
      else
        {
         sign = 1;
        }
      if(u == 2.0)
        {
         return (log(z));
        }
      p -= 2.0;
      x = x + p;
      p = x * polevl(x, gamma_B, 5) / p1evl(x, gamma_C, 6);
      return (log(z) + p);
     }

   if(x > MAXLGM)
     {
      return (sign * infinity());
     }

   if(x >= 1000.0)
     {
      return lgam_large_x(x);
     }

   q = (x - 0.5) * log(x) - x + LS2PI;
   p = 1.0 / (x * x);
   return q + polevl(p, gamma_A, 4) / x;
  }
} // namespace detail

/* Logarithm of Gamma function */
 inline double lgam(double x)
  {
   int sign;
   return detail::lgam_sgn(x, sign);
  }

/* Sign of the Gamma function */
 inline double gammasgn(double x)
  {
   double fx;

   if(isnan(x))
     {
      return x;
     }
   if(x > 0)
     {
      return 1.0;
     }
   else
     {
      fx = floor(x);
      if(x - fx == 0.0)
        {
         return 0.0;
        }
      else
         if(fmod(int(fx), 2)!=0)
           {
            return -1.0;
           }
         else
           {
            return 1.0;
           }
     }
  }

} // namespace cephes
} // namespace special