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