//+------------------------------------------------------------------+
//|                                                         igam.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "error.mqh"
#include "const.mqh"
#include "gamma.mqh"
#include "igam_asymp_coeff.mqh"
#include "lanczos.mqh"
#include "ndtr.mqh"
#include "unity.mqh"
namespace special
{
namespace cephes
{

namespace detail
{

int igam_MAXITER = 2000;
int IGAM = 1;
int IGAMC = 0;
double igam_SMALL = 20;
double igam_LARGE = 200;
double igam_SMALLRATIO = 0.3;
double igam_LARGERATIO = 4.5;

double igam_big = 4.503599627370496e15;
double igam_biginv = 2.22044604925031308085e-16;

/* Compute
 *
 * x^a * exp(-x) / gamma(a)
 *
 * corrected from (15) and (16) in [2] by replacing exp(x - a) with
 * exp(a - x).
 */
inline double igam_fac(double a, double x)
  {
   double ax, fac, res, num;

   if(abs(a - x) > 0.4 * abs(a))
     {
      ax = a * log(x) - x - special::cephes::lgam(a);
      if(ax < -MAXLOG)
        {
         set_error("igam", SF_ERROR_UNDERFLOW, NULL);
         return 0.0;
        }
      return exp(ax);
     }

   fac = a + special::cephes::lanczos_g - 0.5;
   res = sqrt(fac / exp(1)) / special::cephes::lanczos_sum_expg_scaled(a);

   if((a < 200) && (x < 200))
     {
      res *= exp(a - x) * pow(x / fac, a);
     }
   else
     {
      num = x - a - special::cephes::lanczos_g + 0.5;
      res *= exp(a * special::cephes::log1pmx(num / fac) + x * (0.5 - special::cephes::lanczos_g) / fac);
     }

   return res;
  }

/* Compute igamc using DLMF 8.9.2. */
inline double igamc_continued_fraction(double a, double x)
  {
   int i;
   double ans, ax, c, yc, r, t, y, z;
   double pk, pkm1, pkm2, qk, qkm1, qkm2;

   ax = igam_fac(a, x);
   if(ax == 0.0)
     {
      return 0.0;
     }

   /* continued fraction */
   y = 1.0 - a;
   z = x + y + 1.0;
   c = 0.0;
   pkm2 = 1.0;
   qkm2 = x;
   pkm1 = x + 1.0;
   qkm1 = z * x;
   ans = pkm1 / qkm1;

   for(i = 0; i < igam_MAXITER; i++)
     {
      c += 1.0;
      y += 1.0;
      z += 2.0;
      yc = y * c;
      pk = pkm1 * z - pkm2 * yc;
      qk = qkm1 * z - qkm2 * yc;
      if(qk != 0)
        {
         r = pk / qk;
         t = abs((ans - r) / r);
         ans = r;
        }
      else
         t = 1.0;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;
      if(abs(pk) > igam_big)
        {
         pkm2 *= igam_biginv;
         pkm1 *= igam_biginv;
         qkm2 *= igam_biginv;
         qkm1 *= igam_biginv;
        }
      if(t <= MACHEP)
        {
         break;
        }
     }

   return (ans * ax);
  }

/* Compute igam using DLMF 8.11.4. */
inline double igam_series(double a, double x)
  {
   int i;
   double ans, ax, c, r;

   ax = igam_fac(a, x);
   if(ax == 0.0)
     {
      return 0.0;
     }

   /* power series */
   r = a;
   c = 1.0;
   ans = 1.0;

   for(i = 0; i < igam_MAXITER; i++)
     {
      r += 1.0;
      c *= x / r;
      ans += c;
      if(c <= MACHEP * ans)
        {
         break;
        }
     }

   return (ans * ax / a);
  }

/* Compute igamc using DLMF 8.7.3. This is related to the series in
 * igam_series but extra care is taken to avoid cancellation.
 */
inline double igamc_series(double a, double x)
  {
   int n;
   double fac = 1;
   double sum = 0;
   double term, logx;

   for(n = 1; n < igam_MAXITER; n++)
     {
      fac *= -x / n;
      term = fac / (a + n);
      sum += term;
      if(abs(term) <= MACHEP * abs(sum))
        {
         break;
        }
     }

   logx = log(x);
   term = -special::cephes::expm1(a * logx - special::cephes::lgam1p(a));
   return term - exp(a * logx - special::cephes::lgam(a)) * sum;
  }

/* Compute igam/igamc using DLMF 8.12.3/8.12.4. */
inline double asymptotic_series(double a, double x, int func)
  {
   int k, n, sgn;
   int maxpow = 0;
   double lambda = x / a;
   double sigma = (x - a) / a;
   double eta, res, ck, ckterm, term, absterm;
   double absoldterm = infinity();
   double etapow[igam_asymp_coeff_N] = {1};
   double sum = 0;
   double afac = 1;

   if(func == detail::IGAM)
     {
      sgn = -1;
     }
   else
     {
      sgn = 1;
     }

   if(lambda > 1)
     {
      eta = sqrt(-2 * special::cephes::log1pmx(sigma));
     }
   else
      if(lambda < 1)
        {
         eta = -sqrt(-2 * special::cephes::log1pmx(sigma));
        }
      else
        {
         eta = 0;
        }
   res = 0.5 * special::cephes::erfc(sgn * eta * sqrt(a / 2));

   for(k = 0; k < igam_asymp_coeff_K; k++)
     {
      ck = igam_asymp_coeff_d[k][0];
      for(n = 1; n < igam_asymp_coeff_N; n++)
        {
         if(n > maxpow)
           {
            etapow[n] = eta * etapow[n - 1];
            maxpow += 1;
           }
         ckterm = igam_asymp_coeff_d[k][n] * etapow[n];
         ck += ckterm;
         if(abs(ckterm) < MACHEP * abs(ck))
           {
            break;
           }
        }
      term = ck * afac;
      absterm = abs(term);
      if(absterm > absoldterm)
        {
         break;
        }
      sum += term;
      if(absterm < MACHEP * abs(sum))
        {
         break;
        }
      absoldterm = absterm;
      afac /= a;
     }
   res += sgn * exp(-0.5 * a * eta * eta) * sum / sqrt(2 * M_PI * a);

   return res;
  }

} // namespace detail


inline double igam(double a, double x)
  {
   double absxma_a;

   if(x < 0 || a < 0)
     {
      set_error("gammainc", SF_ERROR_DOMAIN, NULL);
      return quiet_NaN();
     }
   else
      if(a == 0)
        {
         if(x > 0)
           {
            return 1;
           }
         else
           {
            return quiet_NaN();
           }
        }
      else
         if(x == 0)
           {
            /* Zero integration limit */
            return 0;
           }
         else
            if(isinf(a))
              {
               if(isinf(x))
                 {
                  return quiet_NaN();
                 }
               return 0;
              }
            else
               if(isinf(x))
                 {
                  return 1;
                 }

   /* Asymptotic regime where a ~ x; see [2]. */
   absxma_a = abs(x - a) / a;
   if((a > detail::igam_SMALL) && (a < detail::igam_LARGE) && (absxma_a < detail::igam_SMALLRATIO))
     {
      return detail::asymptotic_series(a, x, detail::IGAM);
     }
   else
      if((a > detail::igam_LARGE) && (absxma_a < detail::igam_LARGERATIO / sqrt(a)))
        {
         return detail::asymptotic_series(a, x, detail::IGAM);
        }

   if((x > 1.0) && (x > a))
     {
      return (1.0 - igamc(a, x));
     }

   return detail::igam_series(a, x);
  }

double igamc(double a, double x)
  {
   double absxma_a;

   if(x < 0 || a < 0)
     {
      set_error("gammaincc", SF_ERROR_DOMAIN, NULL);
      return quiet_NaN();
     }
   else
      if(a == 0)
        {
         if(x > 0)
           {
            return 0;
           }
         else
           {
            return quiet_NaN();
           }
        }
      else
         if(x == 0)
           {
            return 1;
           }
         else
            if(isinf(a))
              {
               if(isinf(x))
                 {
                  return quiet_NaN();
                 }
               return 1;
              }
            else
               if(isinf(x))
                 {
                  return 0;
                 }

   /* Asymptotic regime where a ~ x; see [2]. */
   absxma_a = abs(x - a) / a;
   if((a > detail::igam_SMALL) && (a < detail::igam_LARGE) && (absxma_a < detail::igam_SMALLRATIO))
     {
      return detail::asymptotic_series(a, x, detail::IGAMC);
     }
   else
      if((a > detail::igam_LARGE) && (absxma_a < detail::igam_LARGERATIO / sqrt(a)))
        {
         return detail::asymptotic_series(a, x, detail::IGAMC);
        }

   /* Everywhere else; see [2]. */
   if(x > 1.1)
     {
      if(x < a)
        {
         return 1.0 - detail::igam_series(a, x);
        }
      else
        {
         return detail::igamc_continued_fraction(a, x);
        }
     }
   else
      if(x <= 0.5)
        {
         if(-0.4 / log(x) < a)
           {
            return 1.0 - detail::igam_series(a, x);
           }
         else
           {
            return detail::igamc_series(a, x);
           }
        }
      else
        {
         if(x * 1.1 < a)
           {
            return 1.0 - detail::igam_series(a, x);
           }
         else
           {
            return detail::igamc_series(a, x);
           }
        }
  }

} // namespace cephes
} // namespace special