417 lignes
9 Kio
MQL5
417 lignes
9 Kio
MQL5
//+------------------------------------------------------------------+
|
|
//| 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
|