dng/AlgLib_ver3.19
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
AlgLib_ver3.19/specialfunctions.mqh

6535 lines
236 KiB
MQL5
Raw Permalink Normal View History

convert
2025-05-30 14:39:48 +02:00
//+------------------------------------------------------------------+
//| specialfunctions.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
#include "ap.mqh"
#include "alglibinternal.mqh"
//+------------------------------------------------------------------+
//| Gamma function |
//+------------------------------------------------------------------+
class CGammaFunc
{
public:
static double GammaFunc(double x);
static double LnGamma(double x,double &sgngam);
private:
static double GammaStirlFunc(double x);
};
//+------------------------------------------------------------------+
//| Gamma function |
//| Input parameters: |
//| X - argument |
//| Domain: |
//| 0 < X < 171.6 |
//| -170 < X < 0, X is not an integer. |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -170,-33 20000 2.3e-15 3.3e-16 |
//| IEEE -33, 33 20000 9.4e-16 2.2e-16 |
//| IEEE 33, 171.6 20000 2.3e-15 3.2e-16 |
//+------------------------------------------------------------------+
double CGammaFunc::GammaFunc(double x)
{
//--- create variables
double p;
double pp;
double q=MathAbs(x);
double qq;
double z;
int i;
double sgngam=1;
//--- check
if(q>33.0)
{
//--- check
if(x<0.0)
{
p=(int)MathFloor(q);
i=(int)MathRound(p);
//--- check
if(i%2==0)
sgngam=-1;
z=q-p;
//--- check
if(z>0.5)
{
p++;
z=q-p;
}
//--- calculation z
z=q*MathSin(M_PI*z);
z=MathAbs(z);
z=M_PI/(z*GammaStirlFunc(q));
}
else
z=GammaStirlFunc(x);
//--- return result
return(sgngam*z);
}
//--- q<=33
z=1;
while(x>=3)
{
x--;
z*=x;
}
//--- x<0
while(x<0)
{
//--- check
if(x>-0.000000001)
return(z/((1+0.5772156649015329*x)*x));
z/=x;
x++;
}
//--- 0<x<2
while(x<2)
{
//--- check
if(x<0.000000001)
return(z/((1+0.5772156649015329*x)*x));
z=z/x;
x=x+1.0;
}
//--- check
if(x==2)
return(z);
//--- calculation
x-=2.0;
pp=1.60119522476751861407E-4;
pp=1.19135147006586384913E-3+x*pp;
pp=1.04213797561761569935E-2+x*pp;
pp=4.76367800457137231464E-2+x*pp;
pp=2.07448227648435975150E-1+x*pp;
pp=4.94214826801497100753E-1+x*pp;
pp=9.99999999999999996796E-1+x*pp;
qq=-2.31581873324120129819E-5;
qq=5.39605580493303397842E-4+x*qq;
qq=-4.45641913851797240494E-3+x*qq;
qq=1.18139785222060435552E-2+x*qq;
qq=3.58236398605498653373E-2+x*qq;
qq=-2.34591795718243348568E-1+x*qq;
qq=7.14304917030273074085E-2+x*qq;
qq=1.00000000000000000320+x*qq;
//--- return result
return(z*pp/qq);
}
//+------------------------------------------------------------------+
//| Natural logarithm of gamma function |
//| Input parameters: |
//| X - argument |
//| Result: |
//| logarithm of the absolute value of the Gamma(X). |
//| Output parameters: |
//| SgnGam - sign(Gamma(X)) |
//| Domain: |
//| 0 < X < 2.55e305 |
//| -2.55e305 < X < 0, X is not an integer. |
//| ACCURACY: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 3 28000 5.4e-16 1.1e-16 |
//| IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 |
//| The error criterion was relative when the function magnitude |
//| was greater than one but absolute when it was less than one. |
//| The following test used the relative error criterion, though |
//| at certain points the relative error could be much higher than |
//| indicated. |
//| IEEE -200, -4 10000 4.8e-16 1.3e-16 |
//+------------------------------------------------------------------+
double CGammaFunc::LnGamma(double x,double &sgngam)
{
//--- create variables
double a;
double b;
double c;
double p;
double q;
double u;
double w;
double z;
int i;
double logpi;
double ls2pi;
double tmp=0;
//--- initialization
sgngam=1;
logpi=1.14472988584940017414;
ls2pi=0.91893853320467274178;
//--- check
if(x<-34.0)
{
q=-x;
w=LnGamma(q,tmp);
p=(int)MathFloor(q);
i=(int)MathRound(p);
//--- check
if(i%2==0)
sgngam=-1;
else
sgngam=1;
z=q-p;
//--- check
if(z>0.5)
{
p++;
z=p-q;
}
z=q*MathSin(M_PI*z);
//--- return result
return(logpi-MathLog(z)-w);
}
//--- x>=-34
if(x<13)
{
z=1;
p=0;
u=x;
//--- cycle
while(u>=3)
{
p--;
u=x+p;
z*=u;
}
//--- cycle
while(u<2)
{
z/=u;
p++;
u=x+p;
}
//--- check
if(z<0)
{
sgngam=-1;
z=-z;
}
else
sgngam=1;
//--- check
if(u==2)
return(MathLog(z));
//--- calculation
p-=2;
x+=p;
b=-1378.25152569120859100;
b=-38801.6315134637840924+x*b;
b=-331612.992738871184744+x*b;
b=-1162370.97492762307383+x*b;
b=-1721737.00820839662146+x*b;
b=-853555.664245765465627+x*b;
c=1;
c=-351.815701436523470549+x*c;
c=-17064.2106651881159223+x*c;
c=-220528.590553854454839+x*c;
c=-1139334.44367982507207+x*c;
c=-2532523.07177582951285+x*c;
c=-2018891.41433532773231+x*c;
p=x*b/c;
//--- return result
return(MathLog(z)+p);
}
//--- x>=13
q=(x-0.5)*MathLog(x)-x+ls2pi;
//--- check
if(x>100000000)
return(q);
//--- change value
p=1/(x*x);
//--- check
if(x>=1000.0)
q+=((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
else
{
a=8.11614167470508450300*0.0001;
a=-(5.95061904284301438324*0.0001)+p*a;
a=7.93650340457716943945*0.0001+p*a;
a=-(2.77777777730099687205*0.001)+p*a;
a=8.33333333333331927722*0.01+p*a;
q+=a/x;
}
//--- return result
return(q);
}
//+------------------------------------------------------------------+
//| Stirling function |
//+------------------------------------------------------------------+
double CGammaFunc::GammaStirlFunc(double x)
{
//--- create variables
double y;
double w;
double v;
double stir;
//--- initialization
w=1/x;
stir=7.87311395793093628397E-4;
stir=-2.29549961613378126380E-4+w*stir;
stir=-2.68132617805781232825E-3+w*stir;
stir=3.47222221605458667310E-3+w*stir;
stir=8.33333333333482257126E-2+w*stir;
w=1+w*stir;
y=MathExp(x);
//--- check
if(x>143.01608)
{
v=MathPow(x,0.5*x-0.25);
y=v*(v/y);
}
else
y=MathPow(x,x-0.5)/y;
//--- return result
return(2.50662827463100050242*y*w);
}
//+------------------------------------------------------------------+
//| Normal distribution |
//+------------------------------------------------------------------+
class CNormalDistr
{
public:
//--- methods
static double ErrorFunction(double x);
static double ErrorFunctionC(double x);
static double NormalDistribution(const double x);
static double NormalPDF(const double x);
static double NormalCDF(const double x);
static double InvErF(const double e);
static double InvNormalDistribution(double y0);
static double InvNormalCDF(const double y0);
static double BivariateNormalPDF(const double x,const double y,const double rho);
static double BivariateNormalCDF(double x,double y,const double rho);
private:
static double BVNIntegrate3(double rangea,double rangeb,double x,double y,double gw,double gx);
static double BVNIntegrate6(double rangea,double rangeb,double x,double y,double s,double gw,double gx);
};
//+------------------------------------------------------------------+
//| Error function |
//| The integral is |
//| x |
//| - |
//| 2 | | 2 |
//| erf(x) = -------- | exp( - t ) dt. |
//| sqrt(pi) | | |
//| - |
//| 0 |
//| For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise |
//| erf(x) = 1 - erfc(x). |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,1 30000 3.7e-16 1.0e-16 |
//+------------------------------------------------------------------+
double CNormalDistr::ErrorFunction(double x)
{
//--- create variables
double xsq=0;
double s=0;
double p=0;
double q=0;
//--- initialization
s=MathSign(x);
//--- change value
x=MathAbs(x);
//--- check
if(x<0.5)
{
//--- calculation
xsq=x*x;
p=0.007547728033418631287834;
p=-0.288805137207594084924010+xsq*p;
p=14.3383842191748205576712+xsq*p;
p=38.0140318123903008244444+xsq*p;
p=3017.82788536507577809226+xsq*p;
p=7404.07142710151470082064+xsq*p;
p=80437.3630960840172832162+xsq*p;
//---
q=38.0190713951939403753468+xsq;
q=658.070155459240506326937+xsq*q;
q=6379.60017324428279487120+xsq*q;
q=34216.5257924628539769006+xsq*q;
q=80437.3630960840172826266+xsq*q;
//--- return result
return(s*1.1283791670955125738961589031*x*p/q);
}
//--- check
if(x>=10)
return(s);
//--- return result
return(s*(1-ErrorFunctionC(x)));
}
//+------------------------------------------------------------------+
//| Complementary error function |
//| 1 - erf(x) = |
//| inf. |
//| - |
//| 2 | | 2 |
//| erfc(x) = -------- | exp( - t ) dt |
//| sqrt(pi) | | |
//| - |
//| x |
//| For small x, erfc(x) = 1 - erf(x); otherwise rational |
//| approximations are computed. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,26.6417 30000 5.7e-14 1.5e-14 |
//+------------------------------------------------------------------+
double CNormalDistr::ErrorFunctionC(double x)
{
//--- create variables
double p=0;
double q=0;
//--- check
if(x<0.0)
return(2-ErrorFunctionC(-x));
//--- check
if(x<0.5)
return(1.0-ErrorFunction(x));
//--- check
if(x>=10)
return(0);
//--- calculation
p=0.5641877825507397413087057563;
p=9.675807882987265400604202961+x*p;
p=77.08161730368428609781633646+x*p;
p=368.5196154710010637133875746+x*p;
p=1143.262070703886173606073338+x*p;
p=2320.439590251635247384768711+x*p;
p=2898.0293292167655611275846+x*p;
p=1826.3348842295112592168999+x*p;
q=17.14980943627607849376131193+x;
q=137.1255960500622202878443578+x*q;
q=661.7361207107653469211984771+x*q;
q=2094.384367789539593790281779+x*q;
q=4429.612803883682726711528526+x*q;
q=6089.5424232724435504633068+x*q;
q=4958.82756472114071495438422+x*q;
q=1826.3348842295112595576438+x*q;
//--- return result
return(MathExp(-CMath::Sqr(x))*p/q);
}
//+------------------------------------------------------------------+
//| Normal distribution function |
//| Returns the area under the Gaussian probability density |
//| function, integrated from minus infinity to x: |
//| x |
//| - |
//| 1 | | 2 |
//| ndtr(x) = --------- | exp( - t /2 ) dt |
//| sqrt(2pi) | | |
//| - |
//| -inf. |
//| = ( 1 + erf(z) ) / 2 |
//| = erfc(z) / 2 |
//| where z = x/sqrt(2). Computation is via the functions |
//| erf and erfc. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -13,0 30000 3.4e-14 6.7e-15 |
//+------------------------------------------------------------------+
double CNormalDistr::NormalDistribution(const double x)
{
return(0.5*(ErrorFunction(x/1.41421356237309504880)+1));
}
//+------------------------------------------------------------------+
//| Normal distribution PDF |
//| Returns Gaussian probability density function: |
//| 1 |
//| f(x) = --------- * exp(-x^2/2) |
//| sqrt(2pi) |
//| Cephes Math Library Release 2.8: June, 2000 |
//| Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier |
//+------------------------------------------------------------------+
double CNormalDistr::NormalPDF(const double x)
{
//--- check
if(!CAp::Assert(MathIsValidNumber(x),"NormalPDF: X is infinite"))
return(0);
double result=MathExp(-(x*x/2))/MathSqrt(2*M_PI);
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Normal distribution CDF |
//| Returns the area under the Gaussian probability density |
//| function, integrated from minus infinity to x: |
//| x |
//| - |
//| 1 | | 2 |
//| ndtr(x) = --------- | exp( - t /2 ) dt |
//| sqrt(2pi) | | |
//| - |
//| -inf. |
//| |
//| = ( 1 + erf(z) ) / 2 |
//| = erfc(z) / 2 |
//| where z = x/sqrt(2). Computation is via the functions erf and |
//| erfc. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -13,0 30000 3.4e-14 6.7e-15 |
//| Cephes Math Library Release 2.8: June, 2000 |
//| Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier |
//+------------------------------------------------------------------+
double CNormalDistr::NormalCDF(const double x)
{
double result=0.5*(ErrorFunction(x/1.41421356237309504880)+1);
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Inverse of the error function |
//+------------------------------------------------------------------+
double CNormalDistr::InvErF(const double e)
{
return(InvNormalDistribution(0.5*(e+1))/MathSqrt(2));
}
//+------------------------------------------------------------------+
//| Inverse of Normal distribution function |
//| Returns the argument, x, for which the area under the |
//| Gaussian probability density function (integrated from |
//| minus infinity to x) is equal to y. |
//| For small arguments 0 < y < exp(-2), the program computes |
//| z = sqrt( -2.0 * log(y) ); then the approximation is |
//| x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). |
//| There are two rational functions P/Q, one for 0 < y < exp(-32) |
//| and the other for y up to exp(-2). For larger arguments, |
//| w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0.125, 1 20000 7.2e-16 1.3e-16 |
//| IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 |
//+------------------------------------------------------------------+
double CNormalDistr::InvNormalDistribution(double y0)
{
//--- check
if(y0<=0)
return(-CMath::m_maxrealnumber);
//--- check
if(y0>=1)
return(CMath::m_maxrealnumber);
//--- create variables
double expm2=0.13533528323661269189;
double s2pi=2.50662827463100050242;
double x;
double y=y0;
double z;
double y2;
double x0;
double x1;
int code=1;
double p0;
double q0;
double p1;
double q1;
double p2;
double q2;
//--- check
if(y>1.0-expm2)
{
y=1.0-y;
code=0;
}
//--- check
if(y>expm2)
{
y-=0.5;
y2=y*y;
p0=-59.9633501014107895267;
p0=98.0010754185999661536+y2*p0;
p0=-56.6762857469070293439+y2*p0;
p0=13.9312609387279679503+y2*p0;
p0=-1.23916583867381258016+y2*p0;
q0=1;
q0=1.95448858338141759834+y2*q0;
q0=4.67627912898881538453+y2*q0;
q0=86.3602421390890590575+y2*q0;
q0=-225.462687854119370527+y2*q0;
q0=200.260212380060660359+y2*q0;
q0=-82.0372256168333339912+y2*q0;
q0=15.9056225126211695515+y2*q0;
q0=-1.18331621121330003142+y2*q0;
x=y+y*y2*p0/q0;
x*=s2pi;
//--- return result
return(x);
}
//--- calculation
x=MathSqrt(-(2.0*MathLog(y)));
x0=x-MathLog(x)/x;
z=1.0/x;
//--- check
if(x<8.0)
{
p1=4.05544892305962419923;
p1=31.5251094599893866154+z*p1;
p1=57.1628192246421288162+z*p1;
p1=44.0805073893200834700+z*p1;
p1=14.6849561928858024014+z*p1;
p1=2.18663306850790267539+z*p1;
p1=-(1.40256079171354495875*0.1)+z*p1;
p1=-(3.50424626827848203418*0.01)+z*p1;
p1=-(8.57456785154685413611*0.0001)+z*p1;
q1=1;
q1=15.7799883256466749731+z*q1;
q1=45.3907635128879210584+z*q1;
q1=41.3172038254672030440+z*q1;
q1=15.0425385692907503408+z*q1;
q1=2.50464946208309415979+z*q1;
q1=-(1.42182922854787788574*0.1)+z*q1;
q1=-(3.80806407691578277194*0.01)+z*q1;
q1=-(9.33259480895457427372*0.0001)+z*q1;
x1=z*p1/q1;
}
else
{
p2=3.23774891776946035970;
p2=6.91522889068984211695+z*p2;
p2=3.93881025292474443415+z*p2;
p2=1.33303460815807542389+z*p2;
p2=2.01485389549179081538*0.1+z*p2;
p2=1.23716634817820021358*0.01+z*p2;
p2=3.01581553508235416007*0.0001+z*p2;
p2=2.65806974686737550832*0.000001+z*p2;
p2=6.23974539184983293730*0.000000001+z*p2;
q2=1;
q2=6.02427039364742014255+z*q2;
q2=3.67983563856160859403+z*q2;
q2=1.37702099489081330271+z*q2;
q2=2.16236993594496635890*0.1+z*q2;
q2=1.34204006088543189037*0.01+z*q2;
q2=3.28014464682127739104*0.0001+z*q2;
q2=2.89247864745380683936*0.000001+z*q2;
q2=6.79019408009981274425*0.000000001+z*q2;
x1=z*p2/q2;
}
x=x0-x1;
//--- check
if(code!=0)
x=-x;
//--- return result
return(x);
}
//+------------------------------------------------------------------+
//| Inverse of Normal CDF |
//| Returns the argument, x, for which the area under the Gaussian |
//| probability density function (integrated from minus infinity to |
//| x) is equal to y. |
//| For small arguments 0 < y < exp(-2), the program computes |
//| z = sqrt( -2.0 * log(y) ); then the approximation is |
//| x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). |
//| There are two rational functions P/Q, one for 0 < y < exp(-32) |
//| and the other for y up to exp(-2). For larger arguments, |
//| w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0.125, 1 20000 7.2e-16 1.3e-16 |
//| IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 |
//| Cephes Math Library Release 2.8: June, 2000 |
//| Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier |
//+------------------------------------------------------------------+
double CNormalDistr::InvNormalCDF(const double y0)
{
//--- create variables
double result=0;
double expm2=0;
double s2pi=0;
double x=0;
double y=0;
double z=0;
double y2=0;
double x0=0;
double x1=0;
int code=0;
double p0=0;
double q0=0;
double p1=0;
double q1=0;
double p2=0;
double q2=0;
expm2=0.13533528323661269189;
s2pi=2.50662827463100050242;
if(y0<=0.0)
{
result=-CMath::m_maxrealnumber;
return(result);
}
if(y0>=1.0)
{
result=CMath::m_maxrealnumber;
return(result);
}
code=1;
y=y0;
if(y>(1.0-expm2))
{
y=1.0-y;
code=0;
}
if(y>expm2)
{
y=y-0.5;
y2=y*y;
p0=-59.9633501014107895267;
p0=98.0010754185999661536+y2*p0;
p0=-56.6762857469070293439+y2*p0;
p0=13.9312609387279679503+y2*p0;
p0=-1.23916583867381258016+y2*p0;
q0=1;
q0=1.95448858338141759834+y2*q0;
q0=4.67627912898881538453+y2*q0;
q0=86.3602421390890590575+y2*q0;
q0=-225.462687854119370527+y2*q0;
q0=200.260212380060660359+y2*q0;
q0=-82.0372256168333339912+y2*q0;
q0=15.9056225126211695515+y2*q0;
q0=-1.18331621121330003142+y2*q0;
x=y+y*y2*p0/q0;
x=x*s2pi;
result=x;
return(result);
}
x=MathSqrt(-(2.0*MathLog(y)));
x0=x-MathLog(x)/x;
z=1.0/x;
if(x<8.0)
{
p1=4.05544892305962419923;
p1=31.5251094599893866154+z*p1;
p1=57.1628192246421288162+z*p1;
p1=44.0805073893200834700+z*p1;
p1=14.6849561928858024014+z*p1;
p1=2.18663306850790267539+z*p1;
p1=-(1.40256079171354495875*0.1)+z*p1;
p1=-(3.50424626827848203418*0.01)+z*p1;
p1=-(8.57456785154685413611*0.0001)+z*p1;
q1=1;
q1=15.7799883256466749731+z*q1;
q1=45.3907635128879210584+z*q1;
q1=41.3172038254672030440+z*q1;
q1=15.0425385692907503408+z*q1;
q1=2.50464946208309415979+z*q1;
q1=-(1.42182922854787788574*0.1)+z*q1;
q1=-(3.80806407691578277194*0.01)+z*q1;
q1=-(9.33259480895457427372*0.0001)+z*q1;
x1=z*p1/q1;
}
else
{
p2=3.23774891776946035970;
p2=6.91522889068984211695+z*p2;
p2=3.93881025292474443415+z*p2;
p2=1.33303460815807542389+z*p2;
p2=2.01485389549179081538*0.1+z*p2;
p2=1.23716634817820021358*0.01+z*p2;
p2=3.01581553508235416007*0.0001+z*p2;
p2=2.65806974686737550832*0.000001+z*p2;
p2=6.23974539184983293730*0.000000001+z*p2;
q2=1;
q2=6.02427039364742014255+z*q2;
q2=3.67983563856160859403+z*q2;
q2=1.37702099489081330271+z*q2;
q2=2.16236993594496635890*0.1+z*q2;
q2=1.34204006088543189037*0.01+z*q2;
q2=3.28014464682127739104*0.0001+z*q2;
q2=2.89247864745380683936*0.000001+z*q2;
q2=6.79019408009981274425*0.000000001+z*q2;
x1=z*p2/q2;
}
x=x0-x1;
if(code!=0)
x=-x;
//--- return result
result=x;
return(result);
}
//+------------------------------------------------------------------+
//| Bivariate normal PDF |
//| Returns probability density function of the bivariate Gaussian |
//| with correlation parameter equal to Rho: |
//| 1 ( x^2 - 2*rho*x*y + y^2 )|
//| f(x,y,rho) = ----------------- * exp( - ----------------------- )|
//| 2pi*sqrt(1-rho^2) ( 2*(1-rho^2) )|
//| with -1<rho<+1 and arbitrary x, y. |
//| This function won't fail as long as Rho is in (-1,+1) range. |
//+------------------------------------------------------------------+
double CNormalDistr::BivariateNormalPDF(const double x,const double y,
const double rho)
{
//--- create variables
double result=0;
double onerho2=0;
//--- check
if(!CAp::Assert(MathIsValidNumber(x),"BivariateNormalCDF: X is infinite"))
return(result);
if(!CAp::Assert(MathIsValidNumber(y),"BivariateNormalCDF: Y is infinite"))
return(result);
if(!CAp::Assert(MathIsValidNumber(rho),"BivariateNormalCDF: Rho is infinite"))
return(result);
if(!CAp::Assert(-1<rho && rho<1.0,"BivariateNormalCDF: Rho is not in (-1,+1) range"))
return(result);
onerho2=(1-rho)*(1+rho);
result=MathExp(-((x*x+y*y-2*rho*x*y)/(2*onerho2)))/(2*M_PI*MathSqrt(onerho2));
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Bivariate normal CDF |
//| Returns the area under the bivariate Gaussian PDF with |
//| correlation parameter equal to Rho, integrated from minus |
//| infinity to (x,y): |
//| x y |
//| - - |
//| 1 | | | | |
//| bvn(x,y,rho) = ------------------- | | f(u,v,rho)*du*dv |
//| 2pi*sqrt(1-rho^2) | | | | |
//| - - |
//| -INF -INF |
//| where |
//| ( u^2 - 2*rho*u*v + v^2 ) |
//| f(u,v,rho) = exp( - ----------------------- ) |
//| ( 2*(1-rho^2) ) |
//| with -1<rho<+1 and arbitrary x, y. |
//| This subroutine uses high-precision approximation scheme proposed|
//| by Alan Genz in "Numerical Computation of Rectangular Bivariate |
//| and Trivariate Normal and t probabilities", which computes CDF |
//| with absolute error roughly equal to 1e-14. |
//| This function won't fail as long as Rho is in (-1,+1) range. |
//+------------------------------------------------------------------+
double CNormalDistr::BivariateNormalCDF(double x,double y,
const double rho)
{
//--- create variables
double result=0;
double rangea=0;
double rangeb=0;
double s=0;
double v=0;
double v0=0;
double v1=0;
double fxys=0;
double ta=0;
double tb=0;
double tc=0;
//--- check
if(!CAp::Assert(MathIsValidNumber(x),"BivariateNormalCDF: X is infinite"))
return(result);
if(!CAp::Assert(MathIsValidNumber(y),"BivariateNormalCDF: Y is infinite"))
return(result);
if(!CAp::Assert(MathIsValidNumber(rho),"BivariateNormalCDF: Rho is infinite"))
return(result);
if(!CAp::Assert(-1<rho && rho<1.0,"BivariateNormalCDF: Rho is not in (-1,+1) range"))
return(result);
if(rho==0.0)
{
result=NormalCDF(x)*NormalCDF(y);
return(result);
}
if(MathAbs(rho)<=0.8)
{
// Rho is small, compute integral using using formula (3) by Alan Genz, integrated
// by means of 10-point Gauss-Legendre quadrature
rangea=0;
rangeb=MathArcsin(rho);
v=BVNIntegrate3(rangea,rangeb,x,y,0.2491470458134028,-0.1252334085114689);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.2491470458134028,0.1252334085114689);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.2334925365383548,-0.3678314989981802);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.2334925365383548,0.3678314989981802);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.2031674267230659,-0.5873179542866175);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.2031674267230659,0.5873179542866175);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.1600783285433462,-0.7699026741943047);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.1600783285433462,0.7699026741943047);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.1069393259953184,-0.9041172563704749);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.1069393259953184,0.9041172563704749);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.0471753363865118,-0.9815606342467192);
v+=BVNIntegrate3(rangea,rangeb,x,y,0.0471753363865118,0.9815606342467192);
v=v*0.5*(rangeb-rangea)/(2*M_PI);
result=NormalCDF(x)*NormalCDF(y)+v;
}
else
{
// Rho is large, compute integral using using formula (6) by Alan Genz, integrated
// by means of 20-point Gauss-Legendre quadrature.
x=-x;
y=-y;
s=MathSign(rho);
if(s>0.0)
fxys=NormalCDF(-MathMax(x,y));
else
fxys=MathMax(0.0,NormalCDF(-x)-NormalCDF(y));
rangea=0;
rangeb=MathSqrt(1-rho*rho);
// Compute first term (analytic integral) from formula (6)
ta=rangeb;
tb=MathAbs(x-s*y);
tc=(4-s*x*y)/8;
v0=ta*(1-tc*(tb*tb-ta*ta)/3)*MathExp(-(tb*tb/(2*ta*ta)))-tb*(1-tc*tb*tb/3)*MathSqrt(2*M_PI)*NormalCDF(-(tb/ta));
v0*=MathExp(-(s*x*y/2))/(2*M_PI);
// Compute second term (numerical integral, 20-point Gauss-Legendre rule) from formula (6)
v1=BVNIntegrate6(rangea,rangeb,x,y,s,0.1527533871307258,-0.0765265211334973);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1527533871307258,0.0765265211334973);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1491729864726037,-0.2277858511416451);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1491729864726037,0.2277858511416451);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1420961093183820,-0.3737060887154195);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1420961093183820,0.3737060887154195);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1316886384491766,-0.5108670019508271);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1316886384491766,0.5108670019508271);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1181945319615184,-0.6360536807265150);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1181945319615184,0.6360536807265150);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1019301198172404,-0.7463319064601508);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.1019301198172404,0.7463319064601508);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0832767415767048,-0.8391169718222188);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0832767415767048,0.8391169718222188);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0626720483341091,-0.9122344282513259);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0626720483341091,0.9122344282513259);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0406014298003869,-0.9639719272779138);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0406014298003869,0.9639719272779138);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0176140071391521,-0.9931285991850949);
v1+=BVNIntegrate6(rangea,rangeb,x,y,s,0.0176140071391521,0.9931285991850949);
v1=v1*0.5*(rangeb-rangea)/(2*M_PI);
result=fxys-s*(v0+v1);
}
result=MathMax(result,0);
result=MathMin(result,1);
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Internal function which computes integrand of formula (3) by |
//| Alan Genz times Gaussian weights (passed by user). |
//+------------------------------------------------------------------+
double CNormalDistr::BVNIntegrate3(double rangea,double rangeb,
double x,double y,double gw,
double gx)
{
//--- create variables
double result=0;
double r=0;
double t2=0;
double dd=0;
double sinr=0;
double cosr=0;
r=(rangeb-rangea)*0.5*gx+(rangeb+rangea)*0.5;
t2=MathTan(0.5*r);
dd=1/(1+t2*t2);
sinr=2*t2*dd;
cosr=(1-t2*t2)*dd;
result=gw*MathExp(-((x*x+y*y-2*x*y*sinr)/(2*cosr*cosr)));
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Internal function which computes integrand of formula (6) by |
//| Alan Genz times Gaussian weights (passed by user). |
//+------------------------------------------------------------------+
double CNormalDistr::BVNIntegrate6(double rangea,double rangeb,
double x,double y,double s,
double gw,double gx)
{
//--- create variables
double result=0;
double r=0;
double exphsk22x2=0;
double exphsk2=0;
double sqrt1x2=0;
double exphsk1sqrt1x2=0;
r=(rangeb-rangea)*0.5*gx+(rangeb+rangea)*0.5;
exphsk22x2=MathExp(-((x-s*y)*(x-s*y)/(2*r*r)));
exphsk2=MathExp(-(x*s*y/2));
sqrt1x2=MathSqrt((1-r)*(1+r));
exphsk1sqrt1x2=MathExp(-(x*s*y/(1+sqrt1x2)));
result=gw*exphsk22x2*(exphsk1sqrt1x2/sqrt1x2-exphsk2*(1+(4-x*y*s)*r*r/8));
//--- create variables
return(result);
}
//+------------------------------------------------------------------+
//| Incomplete gamma function |
//+------------------------------------------------------------------+
class CIncGammaF
{
public:
static double IncompleteGamma(const double a,const double x);
static double IncompleteGammaC(const double a,const double x);
static double InvIncompleteGammaC(const double a,const double y0);
};
//+------------------------------------------------------------------+
//| Incomplete gamma integral |
//| The function is defined by |
//| x |
//| - |
//| 1 | | -t a-1 |
//| igam(a,x) = ----- | e t dt. |
//| - | | |
//| | (a) - |
//| 0 |
//| In this implementation both arguments must be positive. |
//| The integral is evaluated by either a power series or |
//| continued fraction expansion, depending on the relative |
//| values of a and x. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,30 200000 3.6e-14 2.9e-15 |
//| IEEE 0,100 300000 9.9e-14 1.5e-14 |
//+------------------------------------------------------------------+
double CIncGammaF::IncompleteGamma(const double a,const double x)
{
//--- check
if(x<=0 || a<=0)
return(0);
//--- check
if(x>1 && x>a)
return(1-IncompleteGammaC(a,x));
//--- create variables
double igammaepsilon=0.000000000000001;
double ans;
double ax;
double c;
double r;
double tmp=0;
//--- change value
ax=a*MathLog(x)-x-CGammaFunc::LnGamma(a,tmp);
//--- check
if(ax<-709.78271289338399)
return(0);
//--- change values
ax=MathExp(ax);
r=a;
c=1;
ans=1;
//--- cycle
do
{
r++;
c*=x/r;
ans+=c;
}
while(c/ans>igammaepsilon);
//--- return result
return(ans*ax/a);
}
//+------------------------------------------------------------------+
//| Complemented incomplete gamma integral |
//| The function is defined by |
//| igamc(a,x) = 1 - igam(a,x) |
//| inf. |
//| - |
//| 1 | | -t a-1 |
//| = ----- | e t dt. |
//| - | | |
//| | (a) - |
//| x |
//| In this implementation both arguments must be positive. |
//| The integral is evaluated by either a power series or |
//| continued fraction expansion, depending on the relative |
//| values of a and x. |
//| ACCURACY: |
//| Tested at random a, x. |
//| a x Relative error: |
//| arithmetic domain domain # trials peak rms |
//| IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15|
//| IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15|
//+------------------------------------------------------------------+
double CIncGammaF::IncompleteGammaC(const double a,const double x)
{
//--- check
if(x<=0 || a<=0)
return(1);
//--- check
if(x<1 || x<a)
return(1-IncompleteGamma(a,x));
//--- create variables
double igammaepsilon=0.000000000000001;
double igammabignumber=4503599627370496.0;
double igammabignumberinv=2.22044604925031308085*0.0000000000000001;
double ans=0;
double ax=0;
double c=0;
double yc=0;
double r=0;
double t=0;
double y=0;
double z=0;
double pk=0;
double pkm1=0;
double pkm2=0;
double qk=0;
double qkm1=0;
double qkm2=0;
double tmp=0;
//--- change value
ax=a*MathLog(x)-x-CGammaFunc::LnGamma(a,tmp);
//--- check
if(ax<-709.78271289338399)
return(0);
//--- change values
ax=MathExp(ax);
y=1-a;
z=x+y+1;
c=0;
pkm2=1;
qkm2=x;
pkm1=x+1;
qkm1=z*x;
ans=pkm1/qkm1;
//--- cycle
do
{
c++;
y++;
z+=2;
yc=y*c;
pk=pkm1*z-pkm2*yc;
qk=qkm1*z-qkm2*yc;
//--- check
if(qk!=0)
{
r=pk/qk;
t=MathAbs((ans-r)/r);
ans=r;
}
else
t=1;
//--- change values
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
//--- check
if(MathAbs(pk)>igammabignumber)
{
pkm2=pkm2*igammabignumberinv;
pkm1=pkm1*igammabignumberinv;
qkm2=qkm2*igammabignumberinv;
qkm1=qkm1*igammabignumberinv;
}
}
while(t>igammaepsilon);
//--- return result
return(ans*ax);
}
//+------------------------------------------------------------------+
//| Inverse of complemented imcomplete gamma integral |
//| Given p, the function finds x such that |
//| igamc( a, x ) = p. |
//| Starting with the approximate value |
//| 3 |
//| x = a t |
//| where |
//| t = 1 - d - ndtri(p) sqrt(d) |
//| and |
//| d = 1/9a, |
//| the routine performs up to 10 Newton iterations to find the |
//| root of igamc(a,x) - p = 0. |
//| ACCURACY: |
//| Tested at random a, p in the intervals indicated. |
//| a p Relative error: |
//| arithmetic domain domain # trials peak rms |
//| IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15|
//| IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15|
//| IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14|
//+------------------------------------------------------------------+
double CIncGammaF::InvIncompleteGammaC(const double a,const double y0)
{
//--- create variables
double igammaepsilon=0;
double iinvgammabignumber=0;
double x0=0;
double x1=0;
double x=0;
double yl=0;
double yh=0;
double y=0;
double d=0;
double lgm=0;
double dithresh=0;
int i=0;
int dir=0;
double tmp=0;
//--- initialization
igammaepsilon=0.000000000000001;
iinvgammabignumber=4503599627370496.0;
x0=iinvgammabignumber;
yl=0;
x1=0;
yh=1;
dithresh=5*igammaepsilon;
d=1/(9*a);
y=1-d-CNormalDistr::InvNormalDistribution(y0)*MathSqrt(d);
x=a*y*y*y;
lgm=CGammaFunc::LnGamma(a,tmp);
i=0;
//--- cycle
while(i<10)
{
//--- check
if(x>x0 || x<x1)
{
d=0.0625;
break;
}
//--- function call
y=IncompleteGammaC(a,x);
//--- check
if(y<yl || y>yh)
{
d=0.0625;
break;
}
//--- check
if(y<y0)
{
x0=x;
yl=y;
}
else
{
x1=x;
yh=y;
}
//--- change value
d=(a-1)*MathLog(x)-x-lgm;
//--- check
if(d<(double)(-709.78271289338399))
{
d=0.0625;
break;
}
//--- change values
d=-MathExp(d);
d=(y-y0)/d;
//--- check
if(MathAbs(d/x)<igammaepsilon)
return(x);
//--- change values
x=x-d;
i=i+1;
}
//--- check
if(x0==iinvgammabignumber)
{
//--- check
if(x<=0.0)
x=1;
//--- cycle
while(x0==iinvgammabignumber)
{
x=(1+d)*x;
y=IncompleteGammaC(a,x);
//--- check
if(y<y0)
{
x0=x;
yl=y;
break;
}
d=d+d;
}
}
//--- change values
d=0.5;
dir=0;
i=0;
//--- cycle
while(i<400)
{
x=x1+d*(x0-x1);
y=IncompleteGammaC(a,x);
lgm=(x0-x1)/(x1+x0);
//--- check
if(MathAbs(lgm)<dithresh)
break;
//--- change value
lgm=(y-y0)/y0;
//--- check
if(MathAbs(lgm)<dithresh)
break;
//--- check
if(x<=0.0)
break;
//--- check
if(y>=y0)
{
x1=x;
yh=y;
//--- check
if(dir<0)
{
dir=0;
d=0.5;
}
else
{
//--- check
if(dir>1)
d=0.5*d+0.5;
else
d=(y0-yl)/(yh-yl);
}
dir=dir+1;
}
else
{
x0=x;
yl=y;
//--- check
if(dir>0)
{
dir=0;
d=0.5;
}
else
{
//--- check
if(dir<-1)
d=0.5*d;
else
d=(y0-yl)/(yh-yl);
}
dir=dir-1;
}
i=i+1;
}
//--- return result
return(x);
}
//+------------------------------------------------------------------+
//| Airy function |
//+------------------------------------------------------------------+
class CAiryF
{
public:
static void Airy(const double x,double &ai,double &aip,double &bi,double &bip);
};
//+------------------------------------------------------------------+
//| Airy function |
//| Solution of the differential equation |
//| y"(x) = xy. |
//| The function returns the two independent solutions Ai, Bi |
//| and their first derivatives Ai'(x), Bi'(x). |
//| Evaluation is by power series summation for small x, |
//| by rational minimax approximations for large x. |
//| ACCURACY: |
//| Error criterion is absolute when function <= 1, relative |
//| when function > 1, except * denotes relative error criterion. |
//| For large negative x, the absolute error increases as x^1.5. |
//| For large positive x, the relative error increases as x^1.5. |
//| Arithmetic domain function # trials peak rms |
//| IEEE -10, 0 Ai 10000 1.6e-15 2.7e-16 |
//| IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15*|
//| IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16 |
//| IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15*|
//| IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16 |
//| IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16 |
//+------------------------------------------------------------------+
void CAiryF::Airy(const double x,double &ai,double &aip,double &bi,
double &bip)
{
//--- create variables
double z=0;
double zz=0;
double t=0;
double f=0;
double g=0;
double uf=0;
double ug=0;
double k=0;
double zeta=0;
double theta=0;
int domflg=0;
double c1=0;
double c2=0;
double sqrt3=0;
double sqpii=0;
double afn=0;
double afd=0;
double agn=0;
double agd=0;
double apfn=0;
double apfd=0;
double apgn=0;
double apgd=0;
double an=0;
double ad=0;
double apn=0;
double apd=0;
double bn16=0;
double bd16=0;
double bppn=0;
double bppd=0;
//--- initialization
ai=0;
aip=0;
bi=0;
bip=0;
sqpii=5.64189583547756286948E-1;
c1=0.35502805388781723926;
c2=0.258819403792806798405;
sqrt3=1.732050807568877293527;
domflg=0;
//--- check
if(x>25.77)
{
//--- initialization
ai=0;
aip=0;
bi=CMath::m_maxrealnumber;
bip=CMath::m_maxrealnumber;
//--- exit the function
return;
}
//--- check
if(x<-2.09)
{
//--- calculation
domflg=15;
t=MathSqrt(-x);
zeta=-(2.0*x*t/3.0);
t=MathSqrt(t);
k=sqpii/t;
z=1.0/zeta;
zz=z*z;
//--- calculation
afn=-1.31696323418331795333E-1;
afn=afn*zz-6.26456544431912369773E-1;
afn=afn*zz-6.93158036036933542233E-1;
afn=afn*zz-2.79779981545119124951E-1;
afn=afn*zz-4.91900132609500318020E-2;
afn=afn*zz-4.06265923594885404393E-3;
afn=afn*zz-1.59276496239262096340E-4;
afn=afn*zz-2.77649108155232920844E-6;
afn=afn*zz-1.67787698489114633780E-8;
//--- calculation
afd=1.00000000000000000000E0;
afd=afd*zz+1.33560420706553243746E1;
afd=afd*zz+3.26825032795224613948E1;
afd=afd*zz+2.67367040941499554804E1;
afd=afd*zz+9.18707402907259625840E0;
afd=afd*zz+1.47529146771666414581E0;
afd=afd*zz+1.15687173795188044134E-1;
afd=afd*zz+4.40291641615211203805E-3;
afd=afd*zz+7.54720348287414296618E-5;
afd=afd*zz+4.51850092970580378464E-7;
uf=1.0+zz*afn/afd;
//--- calculation
agn=1.97339932091685679179E-2;
agn=agn*zz+3.91103029615688277255E-1;
agn=agn*zz+1.06579897599595591108E0;
agn=agn*zz+9.39169229816650230044E-1;
agn=agn*zz+3.51465656105547619242E-1;
agn=agn*zz+6.33888919628925490927E-2;
agn=agn*zz+5.85804113048388458567E-3;
agn=agn*zz+2.82851600836737019778E-4;
agn=agn*zz+6.98793669997260967291E-6;
agn=agn*zz+8.11789239554389293311E-8;
agn=agn*zz+3.41551784765923618484E-10;
//--- calculation
agd=1.00000000000000000000E0;
agd=agd*zz+9.30892908077441974853E0;
agd=agd*zz+1.98352928718312140417E1;
agd=agd*zz+1.55646628932864612953E1;
agd=agd*zz+5.47686069422975497931E0;
agd=agd*zz+9.54293611618961883998E-1;
agd=agd*zz+8.64580826352392193095E-2;
agd=agd*zz+4.12656523824222607191E-3;
agd=agd*zz+1.01259085116509135510E-4;
agd=agd*zz+1.17166733214413521882E-6;
agd=agd*zz+4.91834570062930015649E-9;
//--- calculation
ug=z*agn/agd;
theta=zeta+0.25*M_PI;
f=MathSin(theta);
g=MathCos(theta);
ai=k*(f*uf-g*ug);
bi=k*(g*uf+f*ug);
apfn=1.85365624022535566142E-1;
apfn=apfn*zz+8.86712188052584095637E-1;
apfn=apfn*zz+9.87391981747398547272E-1;
apfn=apfn*zz+4.01241082318003734092E-1;
apfn=apfn*zz+7.10304926289631174579E-2;
apfn=apfn*zz+5.90618657995661810071E-3;
apfn=apfn*zz+2.33051409401776799569E-4;
apfn=apfn*zz+4.08718778289035454598E-6;
apfn=apfn*zz+2.48379932900442457853E-8;
//--- calculation
apfd=1.00000000000000000000E0;
apfd=apfd*zz+1.47345854687502542552E1;
apfd=apfd*zz+3.75423933435489594466E1;
apfd=apfd*zz+3.14657751203046424330E1;
apfd=apfd*zz+1.09969125207298778536E1;
apfd=apfd*zz+1.78885054766999417817E0;
apfd=apfd*zz+1.41733275753662636873E-1;
apfd=apfd*zz+5.44066067017226003627E-3;
apfd=apfd*zz+9.39421290654511171663E-5;
apfd=apfd*zz+5.65978713036027009243E-7;
uf=1.0+zz*apfn/apfd;
//--- calculation
apgn=-3.55615429033082288335E-2;
apgn=apgn*zz-6.37311518129435504426E-1;
apgn=apgn*zz-1.70856738884312371053E0;
apgn=apgn*zz-1.50221872117316635393E0;
apgn=apgn*zz-5.63606665822102676611E-1;
apgn=apgn*zz-1.02101031120216891789E-1;
apgn=apgn*zz-9.48396695961445269093E-3;
apgn=apgn*zz-4.60325307486780994357E-4;
apgn=apgn*zz-1.14300836484517375919E-5;
apgn=apgn*zz-1.33415518685547420648E-7;
apgn=apgn*zz-5.63803833958893494476E-10;
//--- calculation
apgd=1.00000000000000000000E0;
apgd=apgd*zz+9.85865801696130355144E0;
apgd=apgd*zz+2.16401867356585941885E1;
apgd=apgd*zz+1.73130776389749389525E1;
apgd=apgd*zz+6.17872175280828766327E0;
apgd=apgd*zz+1.08848694396321495475E0;
apgd=apgd*zz+9.95005543440888479402E-2;
apgd=apgd*zz+4.78468199683886610842E-3;
apgd=apgd*zz+1.18159633322838625562E-4;
apgd=apgd*zz+1.37480673554219441465E-6;
apgd=apgd*zz+5.79912514929147598821E-9;
ug=z*apgn/apgd;
k=sqpii*t;
aip=-(k*(g*uf+f*ug));
bip=k*(f*uf-g*ug);
//--- exit the function
return;
}
//--- check
if(x>=(double)(2.09))
{
domflg=5;
t=MathSqrt(x);
zeta=2.0*x*t/3.0;
g=MathExp(zeta);
t=MathSqrt(t);
k=2.0*t*g;
z=1.0/zeta;
//--- calculation
an=3.46538101525629032477E-1;
an=an*z+1.20075952739645805542E1;
an=an*z+7.62796053615234516538E1;
an=an*z+1.68089224934630576269E2;
an=an*z+1.59756391350164413639E2;
an=an*z+7.05360906840444183113E1;
an=an*z+1.40264691163389668864E1;
an=an*z+9.99999999999999995305E-1;
ad=5.67594532638770212846E-1;
ad=ad*z+1.47562562584847203173E1;
ad=ad*z+8.45138970141474626562E1;
ad=ad*z+1.77318088145400459522E2;
ad=ad*z+1.64234692871529701831E2;
ad=ad*z+7.14778400825575695274E1;
ad=ad*z+1.40959135607834029598E1;
ad=ad*z+1.00000000000000000470E0;
//--- calculation
f=an/ad;
ai=sqpii*f/k;
k=-(0.5*sqpii*t/g);
apn=6.13759184814035759225E-1;
apn=apn*z+1.47454670787755323881E1;
apn=apn*z+8.20584123476060982430E1;
apn=apn*z+1.71184781360976385540E2;
apn=apn*z+1.59317847137141783523E2;
apn=apn*z+6.99778599330103016170E1;
apn=apn*z+1.39470856980481566958E1;
apn=apn*z+1.00000000000000000550E0;
apd=3.34203677749736953049E-1;
apd=apd*z+1.11810297306158156705E1;
apd=apd*z+7.11727352147859965283E1;
apd=apd*z+1.58778084372838313640E2;
apd=apd*z+1.53206427475809220834E2;
apd=apd*z+6.86752304592780337944E1;
apd=apd*z+1.38498634758259442477E1;
apd=apd*z+9.99999999999999994502E-1;
f=apn/apd;
aip=f*k;
//--- check
if(x>(double)(8.3203353))
{
//--- calculation
bn16=-2.53240795869364152689E-1;
bn16=bn16*z+5.75285167332467384228E-1;
bn16=bn16*z-3.29907036873225371650E-1;
bn16=bn16*z+6.44404068948199951727E-2;
bn16=bn16*z-3.82519546641336734394E-3;
bd16=1.00000000000000000000E0;
bd16=bd16*z-7.15685095054035237902E0;
bd16=bd16*z+1.06039580715664694291E1;
bd16=bd16*z-5.23246636471251500874E0;
bd16=bd16*z+9.57395864378383833152E-1;
bd16=bd16*z-5.50828147163549611107E-2;
//--- calculation
f=z*bn16/bd16;
k=sqpii*g;
bi=k*(1.0+f)/t;
bppn=4.65461162774651610328E-1;
bppn=bppn*z-1.08992173800493920734E0;
bppn=bppn*z+6.38800117371827987759E-1;
bppn=bppn*z-1.26844349553102907034E-1;
bppn=bppn*z+7.62487844342109852105E-3;
bppd=1.00000000000000000000E0;
bppd=bppd*z-8.70622787633159124240E0;
bppd=bppd*z+1.38993162704553213172E1;
bppd=bppd*z-7.14116144616431159572E0;
bppd=bppd*z+1.34008595960680518666E0;
bppd=bppd*z-7.84273211323341930448E-2;
f=z*bppn/bppd;
bip=k*t*(1.0+f);
//--- exit the function
return;
}
}
//--- change values
f=1.0;
g=x;
t=1.0;
uf=1.0;
ug=x;
k=1.0;
z=x*x*x;
//--- cycle
while(t>CMath::m_machineepsilon)
{
//--- calculation
uf=uf*z;
k=k+1.0;
uf=uf/k;
ug=ug*z;
k=k+1.0;
ug=ug/k;
uf=uf/k;
f=f+uf;
k=k+1.0;
ug=ug/k;
g=g+ug;
t=MathAbs(uf/f);
}
//--- change values
uf=c1*f;
ug=c2*g;
//--- check
if(domflg%2==0)
ai=uf-ug;
//--- check
if(domflg/2%2==0)
bi=sqrt3*(uf+ug);
//--- change values
k=4.0;
uf=x*x/2.0;
ug=z/3.0;
f=uf;
g=1.0+ug;
uf=uf/3.0;
t=1.0;
//--- cycle
while(t>CMath::m_machineepsilon)
{
//--- calculation
uf=uf*z;
ug=ug/k;
k=k+1.0;
ug=ug*z;
uf=uf/k;
f=f+uf;
k=k+1.0;
ug=ug/k;
uf=uf/k;
g=g+ug;
k=k+1.0;
t=MathAbs(ug/g);
}
//--- change values
uf=c1*f;
ug=c2*g;
//--- check
if(domflg/4%2==0)
aip=uf-ug;
//--- check
if(domflg/8%2==0)
bip=sqrt3*(uf+ug);
}
//+------------------------------------------------------------------+
//| Bessel function |
//+------------------------------------------------------------------+
class CBessel
{
public:
static double BesselJ0(double x);
static double BesselJ1(double x);
static double BesselJN(int n,double x);
static double BesselY0(double x);
static double BesselY1(double x);
static double BesselYN(int n,double x);
static double BesselI0(double x);
static double BesselI1(double x);
static double BesselK0(double x);
static double BesselK1(double x);
static double BesselKN(int nn,double x);
private:
static void BesselMFirstCheb(const double c,double &b0,double &b1,double &b2);
static void BesselMNextCheb(const double x,const double c,double &b0,double &b1,double &b2);
static void BesselM1FirstCheb(const double c,double &b0,double &b1,double &b2);
static void BesselM1NextCheb(const double x,const double c,double &b0,double &b1,double &b2);
static void BesselAsympt0(const double x,double &pzero,double &qzero);
static void BesselAsympt1(const double x,double &pzero,double &qzero);
};
//+------------------------------------------------------------------+
//| Bessel function of order zero |
//| Returns Bessel function of order zero of the argument. |
//| The domain is divided into the intervals [0, 5] and |
//| (5, infinity). In the first interval the following rational |
//| approximation is used: |
//| 2 2 |
//| (w - r ) (w - r ) P (w) / Q (w) |
//| 1 2 3 8 |
//| 2 |
//| where w = x and the two r's are zeros of the function. |
//| In the second interval, the Hankel asymptotic expansion |
//| is employed with two rational functions of degree 6/6 |
//| and 7/7. |
//| ACCURACY: |
//| Absolute error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 60000 4.2e-16 1.1e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselJ0(double x)
{
//--- create variables
double xsq=0;
double nn=0;
double pzero=0;
double qzero=0;
double p1=0;
double q1=0;
//--- check
if(x<0.0)
x=-x;
//--- check
if(x>8.0)
{
//--- function call
BesselAsympt0(x,pzero,qzero);
nn=x-M_PI/4;
//--- return result
return(MathSqrt(2/M_PI/x)*(pzero*MathCos(nn)-qzero*MathSin(nn)));
}
//--- calculation
xsq=CMath::Sqr(x);
p1=26857.86856980014981415848441;
p1=-40504123.71833132706360663322+xsq*p1;
p1=25071582855.36881945555156435+xsq*p1;
p1=-8085222034853.793871199468171+xsq*p1;
p1=1434354939140344.111664316553+xsq*p1;
p1=-136762035308817138.6865416609+xsq*p1;
p1=6382059341072356562.289432465+xsq*p1;
p1=-117915762910761053603.8440800+xsq*p1;
p1=493378725179413356181.6813446+xsq*p1;
q1=1.0;
q1=1363.063652328970604442810507+xsq*q1;
q1=1114636.098462985378182402543+xsq*q1;
q1=669998767.2982239671814028660+xsq*q1;
q1=312304311494.1213172572469442+xsq*q1;
q1=112775673967979.8507056031594+xsq*q1;
q1=30246356167094626.98627330784+xsq*q1;
q1=5428918384092285160.200195092+xsq*q1;
q1=493378725179413356211.3278438+xsq*q1;
//--- return result
return(p1/q1);
}
//+------------------------------------------------------------------+
//| Bessel function of order one |
//| Returns Bessel function of order one of the argument. |
//| The domain is divided into the intervals [0, 8] and |
//| (8, infinity). In the first interval a 24 term Chebyshev |
//| expansion is used. In the second, the asymptotic |
//| trigonometric representation is employed using two |
//| rational functions of degree 5/5. |
//| ACCURACY: |
//| Absolute error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 2.6e-16 1.1e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselJ1(double x)
{
//--- create variables
double result=0;
double s=0;
double xsq=0;
double nn=0;
double pzero=0;
double qzero=0;
double p1=0;
double q1=0;
//--- calculation
s=MathSign(x);
//--- check
if(x<0.0)
x=-x;
//--- check
if(x>8.0)
{
//--- function call
BesselAsympt1(x,pzero,qzero);
nn=x-3*M_PI/4;
result=MathSqrt(2/M_PI/x)*(pzero*MathCos(nn)-qzero*MathSin(nn));
//--- check
if(s<0.0)
result=-result;
//--- return result
return(result);
}
//--- calculation
xsq=CMath::Sqr(x);
p1=2701.122710892323414856790990;
p1=-4695753.530642995859767162166+xsq*p1;
p1=3413234182.301700539091292655+xsq*p1;
p1=-1322983480332.126453125473247+xsq*p1;
p1=290879526383477.5409737601689+xsq*p1;
p1=-35888175699101060.50743641413+xsq*p1;
p1=2316433580634002297.931815435+xsq*p1;
p1=-66721065689249162980.20941484+xsq*p1;
p1=581199354001606143928.050809+xsq*p1;
q1=1.0;
q1=1606.931573481487801970916749+xsq*q1;
q1=1501793.594998585505921097578+xsq*q1;
q1=1013863514.358673989967045588+xsq*q1;
q1=524371026216.7649715406728642+xsq*q1;
q1=208166122130760.7351240184229+xsq*q1;
q1=60920613989175217.46105196863+xsq*q1;
q1=11857707121903209998.37113348+xsq*q1;
q1=1162398708003212287858.529400+xsq*q1;
//--- return result
return(s*x*p1/q1);
}
//+------------------------------------------------------------------+
//| Bessel function of integer order |
//| Returns Bessel function of order n, where n is a |
//| (possibly negative) integer. |
//| The ratio of jn(x) to j0(x) is computed by backward |
//| recurrence. First the ratio jn/jn-1 is found by a |
//| continued fraction expansion. Then the recurrence |
//| relating successive orders is applied until j0 or j1 is |
//| reached. |
//| If n = 0 or 1 the routine for j0 or j1 is called |
//| directly. |
//| ACCURACY: |
//| Absolute error: |
//| arithmetic range # trials peak rms |
//| IEEE 0, 30 5000 4.4e-16 7.9e-17 |
//| Not suitable for large n or x. Use jv() (fractional order) |
//| instead. |
//+------------------------------------------------------------------+
double CBessel::BesselJN(int n,double x)
{
//--- create variables
double result=0;
double pkm2=0;
double pkm1=0;
double pk=0;
double xk=0;
double r=0;
double ans=0;
int k=0;
int sg=0;
//--- check
if(n<0)
{
n=-n;
//--- check
if(n%2==0)
sg=1;
else
sg=-1;
}
else
sg=1;
//--- check
if(x<0.0)
{
//--- check
if(n%2!=0)
sg=-sg;
x=-x;
}
//--- check
if(n==0)
return(sg*BesselJ0(x));
//--- check
if(n==1)
return(sg*BesselJ1(x));
//--- check
if(n==2)
{
//--- check
if(x==0.0)
return(0);
else
return(sg*(2.0*BesselJ1(x)/x-BesselJ0(x)));
}
//--- check
if(x<CMath::m_machineepsilon)
return(0);
//--- change values
k=53;
pk=2*(n+k);
ans=pk;
xk=x*x;
//--- cycle
do
{
pk=pk-2.0;
ans=pk-xk/ans;
k=k-1;
}
while(k!=0);
//--- change values
ans=x/ans;
pk=1.0;
pkm1=1.0/ans;
k=n-1;
r=2*k;
//--- cycle
do
{
pkm2=(pkm1*r-pk*x)/x;
pk=pkm1;
pkm1=pkm2;
r=r-2.0;
k=k-1;
}
while(k!=0);
//--- check
if(MathAbs(pk)>MathAbs(pkm1))
ans=BesselJ1(x)/pk;
else
ans=BesselJ0(x)/pkm1;
//--- return result
return(sg*ans);
}
//+------------------------------------------------------------------+
//| Bessel function of the second kind, order zero |
//| Returns Bessel function of the second kind, of order |
//| zero, of the argument. |
//| The domain is divided into the intervals [0, 5] and |
//| (5, infinity). In the first interval a rational approximation |
//| R(x) is employed to compute |
//| y0(x) = R(x) + 2 * log(x) * j0(x) / PI. |
//| Thus a call to j0() is required. |
//| In the second interval, the Hankel asymptotic expansion |
//| is employed with two rational functions of degree 6/6 |
//| and 7/7. |
//| ACCURACY: |
//| Absolute error, when y0(x) < 1; else relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 1.3e-15 1.6e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselY0(double x)
{
//--- create variables
double nn=0;
double xsq=0;
double pzero=0;
double qzero=0;
double p4=0;
double q4=0;
//--- check
if(x>8.0)
{
//--- function call
BesselAsympt0(x,pzero,qzero);
nn=x-M_PI/4;
//--- return result
return(MathSqrt(2/M_PI/x)*(pzero*MathSin(nn)+qzero*MathCos(nn)));
}
//--- calculation
xsq=CMath::Sqr(x);
p4=-41370.35497933148554125235152;
p4=59152134.65686889654273830069+xsq*p4;
p4=-34363712229.79040378171030138+xsq*p4;
p4=10255208596863.94284509167421+xsq*p4;
p4=-1648605817185729.473122082537+xsq*p4;
p4=137562431639934407.8571335453+xsq*p4;
p4=-5247065581112764941.297350814+xsq*p4;
p4=65874732757195549259.99402049+xsq*p4;
p4=-27502866786291095837.01933175+xsq*p4;
q4=1.0;
q4=1282.452772478993804176329391+xsq*q4;
q4=1001702.641288906265666651753+xsq*q4;
q4=579512264.0700729537480087915+xsq*q4;
q4=261306575504.1081249568482092+xsq*q4;
q4=91620380340751.85262489147968+xsq*q4;
q4=23928830434997818.57439356652+xsq*q4;
q4=4192417043410839973.904769661+xsq*q4;
q4=372645883898616588198.9980+xsq*q4;
//--- return result
return(p4/q4+2/M_PI*BesselJ0(x)*MathLog(x));
}
//+------------------------------------------------------------------+
//| Bessel function of second kind of order one |
//| Returns Bessel function of the second kind of order one |
//| of the argument. |
//| The domain is divided into the intervals [0, 8] and |
//| (8, infinity). In the first interval a 25 term Chebyshev |
//| expansion is used, and a call to j1() is required. |
//| In the second, the asymptotic trigonometric representation |
//| is employed using two rational functions of degree 5/5. |
//| ACCURACY: |
//| Absolute error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 1.0e-15 1.3e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselY1(double x)
{
//--- create variables
double nn=0;
double xsq=0;
double pzero=0;
double qzero=0;
double p4=0;
double q4=0;
//--- check
if(x>8.0)
{
//--- function call
BesselAsympt1(x,pzero,qzero);
nn=x-3*M_PI/4;
//--- return result
return(MathSqrt(2/M_PI/x)*(pzero*MathSin(nn)+qzero*MathCos(nn)));
}
//--- calculation
xsq=CMath::Sqr(x);
p4=-2108847.540133123652824139923;
p4=3639488548.124002058278999428+xsq*p4;
p4=-2580681702194.450950541426399+xsq*p4;
p4=956993023992168.3481121552788+xsq*p4;
p4=-196588746272214065.8820322248+xsq*p4;
p4=21931073399177975921.11427556+xsq*p4;
p4=-1212297555414509577913.561535+xsq*p4;
p4=26554738314348543268942.48968+xsq*p4;
p4=-99637534243069222259967.44354+xsq*p4;
q4=1.0;
q4=1612.361029677000859332072312+xsq*q4;
q4=1563282.754899580604737366452+xsq*q4;
q4=1128686837.169442121732366891+xsq*q4;
q4=646534088126.5275571961681500+xsq*q4;
q4=297663212564727.6729292742282+xsq*q4;
q4=108225825940881955.2553850180+xsq*q4;
q4=29549879358971486742.90758119+xsq*q4;
q4=5435310377188854170800.653097+xsq*q4;
q4=508206736694124324531442.4152+xsq*q4;
//--- return result
return(x*p4/q4+2/M_PI*(BesselJ1(x)*MathLog(x)-1/x));
}
//+------------------------------------------------------------------+
//| Bessel function of second kind of integer order |
//| Returns Bessel function of order n, where n is a |
//| (possibly negative) integer. |
//| The function is evaluated by forward recurrence on |
//| n, starting with values computed by the routines |
//| y0() and y1(). |
//| If n = 0 or 1 the routine for y0 or y1 is called |
//| directly. |
//| ACCURACY: |
//| Absolute error, except relative |
//| when y > 1: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 3.4e-15 4.3e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselYN(int n,double x)
{
//--- create variables
double a=0;
double b=0;
double tmp=0;
double s=0;
//--- initialization
s=1;
//--- check
if(n<0)
{
n=-n;
//--- check
if(n%2!=0)
s=-1;
}
//--- check
if(n==0)
return(BesselY0(x));
//--- check
if(n==1)
return(s*BesselY1(x));
//--- calculation
a=BesselY0(x);
b=BesselY1(x);
for(int i=1; i<=n-1; i++)
{
tmp=b;
b=2*i/x*b-a;
a=tmp;
}
//--- return result
return(s*b);
}
//+------------------------------------------------------------------+
//| Modified Bessel function of order zero |
//| Returns modified Bessel function of order zero of the |
//| argument. |
//| The function is defined as i0(x) = j0( ix ). |
//| The range is partitioned into the two intervals [0,8] and |
//| (8, infinity). Chebyshev polynomial expansions are employed |
//| in each interval. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,30 30000 5.8e-16 1.4e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselI0(double x)
{
//--- create variables
double y=0;
double v=0;
double z=0;
double b0=0;
double b1=0;
double b2=0;
//--- check
if(x<0.0)
x=-x;
//--- check
if(x<=8.0)
{
y=x/2.0-2.0;
//--- function calls
BesselMFirstCheb(-4.41534164647933937950E-18,b0,b1,b2);
BesselMNextCheb(y,3.33079451882223809783E-17,b0,b1,b2);
BesselMNextCheb(y,-2.43127984654795469359E-16,b0,b1,b2);
BesselMNextCheb(y,1.71539128555513303061E-15,b0,b1,b2);
BesselMNextCheb(y,-1.16853328779934516808E-14,b0,b1,b2);
BesselMNextCheb(y,7.67618549860493561688E-14,b0,b1,b2);
BesselMNextCheb(y,-4.85644678311192946090E-13,b0,b1,b2);
BesselMNextCheb(y,2.95505266312963983461E-12,b0,b1,b2);
BesselMNextCheb(y,-1.72682629144155570723E-11,b0,b1,b2);
BesselMNextCheb(y,9.67580903537323691224E-11,b0,b1,b2);
BesselMNextCheb(y,-5.18979560163526290666E-10,b0,b1,b2);
BesselMNextCheb(y,2.65982372468238665035E-9,b0,b1,b2);
BesselMNextCheb(y,-1.30002500998624804212E-8,b0,b1,b2);
BesselMNextCheb(y,6.04699502254191894932E-8,b0,b1,b2);
BesselMNextCheb(y,-2.67079385394061173391E-7,b0,b1,b2);
BesselMNextCheb(y,1.11738753912010371815E-6,b0,b1,b2);
BesselMNextCheb(y,-4.41673835845875056359E-6,b0,b1,b2);
BesselMNextCheb(y,1.64484480707288970893E-5,b0,b1,b2);
BesselMNextCheb(y,-5.75419501008210370398E-5,b0,b1,b2);
BesselMNextCheb(y,1.88502885095841655729E-4,b0,b1,b2);
BesselMNextCheb(y,-5.76375574538582365885E-4,b0,b1,b2);
BesselMNextCheb(y,1.63947561694133579842E-3,b0,b1,b2);
BesselMNextCheb(y,-4.32430999505057594430E-3,b0,b1,b2);
BesselMNextCheb(y,1.05464603945949983183E-2,b0,b1,b2);
BesselMNextCheb(y,-2.37374148058994688156E-2,b0,b1,b2);
BesselMNextCheb(y,4.93052842396707084878E-2,b0,b1,b2);
BesselMNextCheb(y,-9.49010970480476444210E-2,b0,b1,b2);
BesselMNextCheb(y,1.71620901522208775349E-1,b0,b1,b2);
BesselMNextCheb(y,-3.04682672343198398683E-1,b0,b1,b2);
BesselMNextCheb(y,6.76795274409476084995E-1,b0,b1,b2);
//--- calculation
v=0.5*(b0-b2);
//--- return result
return(MathExp(x)*v);
}
//--- change value
z=32.0/x-2.0;
//--- function calls
BesselMFirstCheb(-7.23318048787475395456E-18,b0,b1,b2);
BesselMNextCheb(z,-4.83050448594418207126E-18,b0,b1,b2);
BesselMNextCheb(z,4.46562142029675999901E-17,b0,b1,b2);
BesselMNextCheb(z,3.46122286769746109310E-17,b0,b1,b2);
BesselMNextCheb(z,-2.82762398051658348494E-16,b0,b1,b2);
BesselMNextCheb(z,-3.42548561967721913462E-16,b0,b1,b2);
BesselMNextCheb(z,1.77256013305652638360E-15,b0,b1,b2);
BesselMNextCheb(z,3.81168066935262242075E-15,b0,b1,b2);
BesselMNextCheb(z,-9.55484669882830764870E-15,b0,b1,b2);
BesselMNextCheb(z,-4.15056934728722208663E-14,b0,b1,b2);
BesselMNextCheb(z,1.54008621752140982691E-14,b0,b1,b2);
BesselMNextCheb(z,3.85277838274214270114E-13,b0,b1,b2);
BesselMNextCheb(z,7.18012445138366623367E-13,b0,b1,b2);
BesselMNextCheb(z,-1.79417853150680611778E-12,b0,b1,b2);
BesselMNextCheb(z,-1.32158118404477131188E-11,b0,b1,b2);
BesselMNextCheb(z,-3.14991652796324136454E-11,b0,b1,b2);
BesselMNextCheb(z,1.18891471078464383424E-11,b0,b1,b2);
BesselMNextCheb(z,4.94060238822496958910E-10,b0,b1,b2);
BesselMNextCheb(z,3.39623202570838634515E-9,b0,b1,b2);
BesselMNextCheb(z,2.26666899049817806459E-8,b0,b1,b2);
BesselMNextCheb(z,2.04891858946906374183E-7,b0,b1,b2);
BesselMNextCheb(z,2.89137052083475648297E-6,b0,b1,b2);
BesselMNextCheb(z,6.88975834691682398426E-5,b0,b1,b2);
BesselMNextCheb(z,3.36911647825569408990E-3,b0,b1,b2);
BesselMNextCheb(z,8.04490411014108831608E-1,b0,b1,b2);
//--- calculation
v=0.5*(b0-b2);
//--- return result
return(MathExp(x)*v/MathSqrt(x));
}
//+------------------------------------------------------------------+
//| Modified Bessel function of order one |
//| Returns modified Bessel function of order one of the |
//| argument. |
//| The function is defined as i1(x) = -i j1( ix ). |
//| The range is partitioned into the two intervals [0,8] and |
//| (8, infinity). Chebyshev polynomial expansions are employed |
//| in each interval. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 1.9e-15 2.1e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselI1(double x)
{
//--- create variables
double y=0;
double z=0;
double v=0;
double b0=0;
double b1=0;
double b2=0;
//--- initialization
z=MathAbs(x);
//--- check
if(z<=8.0)
{
//--- initialization
y=z/2.0-2.0;
//--- function call
BesselM1FirstCheb(2.77791411276104639959E-18,b0,b1,b2);
BesselM1NextCheb(y,-2.11142121435816608115E-17,b0,b1,b2);
BesselM1NextCheb(y,1.55363195773620046921E-16,b0,b1,b2);
BesselM1NextCheb(y,-1.10559694773538630805E-15,b0,b1,b2);
BesselM1NextCheb(y,7.60068429473540693410E-15,b0,b1,b2);
BesselM1NextCheb(y,-5.04218550472791168711E-14,b0,b1,b2);
BesselM1NextCheb(y,3.22379336594557470981E-13,b0,b1,b2);
BesselM1NextCheb(y,-1.98397439776494371520E-12,b0,b1,b2);
BesselM1NextCheb(y,1.17361862988909016308E-11,b0,b1,b2);
BesselM1NextCheb(y,-6.66348972350202774223E-11,b0,b1,b2);
BesselM1NextCheb(y,3.62559028155211703701E-10,b0,b1,b2);
BesselM1NextCheb(y,-1.88724975172282928790E-9,b0,b1,b2);
BesselM1NextCheb(y,9.38153738649577178388E-9,b0,b1,b2);
BesselM1NextCheb(y,-4.44505912879632808065E-8,b0,b1,b2);
BesselM1NextCheb(y,2.00329475355213526229E-7,b0,b1,b2);
BesselM1NextCheb(y,-8.56872026469545474066E-7,b0,b1,b2);
BesselM1NextCheb(y,3.47025130813767847674E-6,b0,b1,b2);
BesselM1NextCheb(y,-1.32731636560394358279E-5,b0,b1,b2);
BesselM1NextCheb(y,4.78156510755005422638E-5,b0,b1,b2);
BesselM1NextCheb(y,-1.61760815825896745588E-4,b0,b1,b2);
BesselM1NextCheb(y,5.12285956168575772895E-4,b0,b1,b2);
BesselM1NextCheb(y,-1.51357245063125314899E-3,b0,b1,b2);
BesselM1NextCheb(y,4.15642294431288815669E-3,b0,b1,b2);
BesselM1NextCheb(y,-1.05640848946261981558E-2,b0,b1,b2);
BesselM1NextCheb(y,2.47264490306265168283E-2,b0,b1,b2);
BesselM1NextCheb(y,-5.29459812080949914269E-2,b0,b1,b2);
BesselM1NextCheb(y,1.02643658689847095384E-1,b0,b1,b2);
BesselM1NextCheb(y,-1.76416518357834055153E-1,b0,b1,b2);
BesselM1NextCheb(y,2.52587186443633654823E-1,b0,b1,b2);
//--- calculation
v=0.5*(b0-b2);
z=v*z*MathExp(z);
}
else
{
//--- initialization
y=32.0/z-2.0;
//--- function calls
BesselM1FirstCheb(7.51729631084210481353E-18,b0,b1,b2);
BesselM1NextCheb(y,4.41434832307170791151E-18,b0,b1,b2);
BesselM1NextCheb(y,-4.65030536848935832153E-17,b0,b1,b2);
BesselM1NextCheb(y,-3.20952592199342395980E-17,b0,b1,b2);
BesselM1NextCheb(y,2.96262899764595013876E-16,b0,b1,b2);
BesselM1NextCheb(y,3.30820231092092828324E-16,b0,b1,b2);
BesselM1NextCheb(y,-1.88035477551078244854E-15,b0,b1,b2);
BesselM1NextCheb(y,-3.81440307243700780478E-15,b0,b1,b2);
BesselM1NextCheb(y,1.04202769841288027642E-14,b0,b1,b2);
BesselM1NextCheb(y,4.27244001671195135429E-14,b0,b1,b2);
BesselM1NextCheb(y,-2.10154184277266431302E-14,b0,b1,b2);
BesselM1NextCheb(y,-4.08355111109219731823E-13,b0,b1,b2);
BesselM1NextCheb(y,-7.19855177624590851209E-13,b0,b1,b2);
BesselM1NextCheb(y,2.03562854414708950722E-12,b0,b1,b2);
BesselM1NextCheb(y,1.41258074366137813316E-11,b0,b1,b2);
BesselM1NextCheb(y,3.25260358301548823856E-11,b0,b1,b2);
BesselM1NextCheb(y,-1.89749581235054123450E-11,b0,b1,b2);
BesselM1NextCheb(y,-5.58974346219658380687E-10,b0,b1,b2);
BesselM1NextCheb(y,-3.83538038596423702205E-9,b0,b1,b2);
BesselM1NextCheb(y,-2.63146884688951950684E-8,b0,b1,b2);
BesselM1NextCheb(y,-2.51223623787020892529E-7,b0,b1,b2);
BesselM1NextCheb(y,-3.88256480887769039346E-6,b0,b1,b2);
BesselM1NextCheb(y,-1.10588938762623716291E-4,b0,b1,b2);
BesselM1NextCheb(y,-9.76109749136146840777E-3,b0,b1,b2);
BesselM1NextCheb(y,7.78576235018280120474E-1,b0,b1,b2);
//--- calculation
v=0.5*(b0-b2);
z=v*MathExp(z)/MathSqrt(z);
}
//--- check
if(x<0.0)
z=-z;
//--- return result
return(z);
}
//+------------------------------------------------------------------+
//| Modified Bessel function, second kind, order zero |
//| Returns modified Bessel function of the second kind |
//| of order zero of the argument. |
//| The range is partitioned into the two intervals [0,8] and |
//| (8, infinity). Chebyshev polynomial expansions are employed |
//| in each interval. |
//| ACCURACY: |
//| Tested at 2000 random points between 0 and 8. Peak absolute |
//| error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 1.2e-15 1.6e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselK0(double x)
{
//--- create variables
double result=0;
double y=0;
double z=0;
double v=0;
double b0=0;
double b1=0;
double b2=0;
//--- check
if(!CAp::Assert(x>0.0,__FUNCTION__+": x<=0"))
return(EMPTY_VALUE);
//--- check
if(x<=2.0)
{
//--- initialization
y=x*x-2.0;
//--- function calls
BesselMFirstCheb(1.37446543561352307156E-16,b0,b1,b2);
BesselMNextCheb(y,4.25981614279661018399E-14,b0,b1,b2);
BesselMNextCheb(y,1.03496952576338420167E-11,b0,b1,b2);
BesselMNextCheb(y,1.90451637722020886025E-9,b0,b1,b2);
BesselMNextCheb(y,2.53479107902614945675E-7,b0,b1,b2);
BesselMNextCheb(y,2.28621210311945178607E-5,b0,b1,b2);
BesselMNextCheb(y,1.26461541144692592338E-3,b0,b1,b2);
BesselMNextCheb(y,3.59799365153615016266E-2,b0,b1,b2);
BesselMNextCheb(y,3.44289899924628486886E-1,b0,b1,b2);
BesselMNextCheb(y,-5.35327393233902768720E-1,b0,b1,b2);
v=0.5*(b0-b2);
v=v-MathLog(0.5*x)*BesselI0(x);
}
else
{
//--- initialization
z=8.0/x-2.0;
//--- function call
BesselMFirstCheb(5.30043377268626276149E-18,b0,b1,b2);
BesselMNextCheb(z,-1.64758043015242134646E-17,b0,b1,b2);
BesselMNextCheb(z,5.21039150503902756861E-17,b0,b1,b2);
BesselMNextCheb(z,-1.67823109680541210385E-16,b0,b1,b2);
BesselMNextCheb(z,5.51205597852431940784E-16,b0,b1,b2);
BesselMNextCheb(z,-1.84859337734377901440E-15,b0,b1,b2);
BesselMNextCheb(z,6.34007647740507060557E-15,b0,b1,b2);
BesselMNextCheb(z,-2.22751332699166985548E-14,b0,b1,b2);
BesselMNextCheb(z,8.03289077536357521100E-14,b0,b1,b2);
BesselMNextCheb(z,-2.98009692317273043925E-13,b0,b1,b2);
BesselMNextCheb(z,1.14034058820847496303E-12,b0,b1,b2);
BesselMNextCheb(z,-4.51459788337394416547E-12,b0,b1,b2);
BesselMNextCheb(z,1.85594911495471785253E-11,b0,b1,b2);
BesselMNextCheb(z,-7.95748924447710747776E-11,b0,b1,b2);
BesselMNextCheb(z,3.57739728140030116597E-10,b0,b1,b2);
BesselMNextCheb(z,-1.69753450938905987466E-9,b0,b1,b2);
BesselMNextCheb(z,8.57403401741422608519E-9,b0,b1,b2);
BesselMNextCheb(z,-4.66048989768794782956E-8,b0,b1,b2);
BesselMNextCheb(z,2.76681363944501510342E-7,b0,b1,b2);
BesselMNextCheb(z,-1.83175552271911948767E-6,b0,b1,b2);
BesselMNextCheb(z,1.39498137188764993662E-5,b0,b1,b2);
BesselMNextCheb(z,-1.28495495816278026384E-4,b0,b1,b2);
BesselMNextCheb(z,1.56988388573005337491E-3,b0,b1,b2);
BesselMNextCheb(z,-3.14481013119645005427E-2,b0,b1,b2);
BesselMNextCheb(z,2.44030308206595545468E0,b0,b1,b2);
//--- calculation
v=0.5*(b0-b2);
v=v*MathExp(-x)/MathSqrt(x);
}
//--- return result
return(v);
}
//+------------------------------------------------------------------+
//| Modified Bessel function, second kind, order one |
//| Computes the modified Bessel function of the second kind |
//| of order one of the argument. |
//| The range is partitioned into the two intervals [0,2] and |
//| (2, infinity). Chebyshev polynomial expansions are employed |
//| in each interval. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 30000 1.2e-15 1.6e-16 |
//+------------------------------------------------------------------+
double CBessel::BesselK1(double x)
{
//--- create variables
double result=0;
double y=0;
double z=0;
double v=0;
double b0=0;
double b1=0;
double b2=0;
//--- initialization
z=0.5*x;
//--- check
if(!CAp::Assert(z>0.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- check
if(x<=2.0)
{
//--- initialization
y=x*x-2.0;
//--- function calls
BesselM1FirstCheb(-7.02386347938628759343E-18,b0,b1,b2);
BesselM1NextCheb(y,-2.42744985051936593393E-15,b0,b1,b2);
BesselM1NextCheb(y,-6.66690169419932900609E-13,b0,b1,b2);
BesselM1NextCheb(y,-1.41148839263352776110E-10,b0,b1,b2);
BesselM1NextCheb(y,-2.21338763073472585583E-8,b0,b1,b2);
BesselM1NextCheb(y,-2.43340614156596823496E-6,b0,b1,b2);
BesselM1NextCheb(y,-1.73028895751305206302E-4,b0,b1,b2);
BesselM1NextCheb(y,-6.97572385963986435018E-3,b0,b1,b2);
BesselM1NextCheb(y,-1.22611180822657148235E-1,b0,b1,b2);
BesselM1NextCheb(y,-3.53155960776544875667E-1,b0,b1,b2);
BesselM1NextCheb(y,1.52530022733894777053E0,b0,b1,b2);
v=0.5*(b0-b2);
//--- get result
result=MathLog(z)*BesselI1(x)+v/x;
}
else
{
//--- initialization
y=8.0/x-2.0;
//--- function calls
BesselM1FirstCheb(-5.75674448366501715755E-18,b0,b1,b2);
BesselM1NextCheb(y,1.79405087314755922667E-17,b0,b1,b2);
BesselM1NextCheb(y,-5.68946255844285935196E-17,b0,b1,b2);
BesselM1NextCheb(y,1.83809354436663880070E-16,b0,b1,b2);
BesselM1NextCheb(y,-6.05704724837331885336E-16,b0,b1,b2);
BesselM1NextCheb(y,2.03870316562433424052E-15,b0,b1,b2);
BesselM1NextCheb(y,-7.01983709041831346144E-15,b0,b1,b2);
BesselM1NextCheb(y,2.47715442448130437068E-14,b0,b1,b2);
BesselM1NextCheb(y,-8.97670518232499435011E-14,b0,b1,b2);
BesselM1NextCheb(y,3.34841966607842919884E-13,b0,b1,b2);
BesselM1NextCheb(y,-1.28917396095102890680E-12,b0,b1,b2);
BesselM1NextCheb(y,5.13963967348173025100E-12,b0,b1,b2);
BesselM1NextCheb(y,-2.12996783842756842877E-11,b0,b1,b2);
BesselM1NextCheb(y,9.21831518760500529508E-11,b0,b1,b2);
BesselM1NextCheb(y,-4.19035475934189648750E-10,b0,b1,b2);
BesselM1NextCheb(y,2.01504975519703286596E-9,b0,b1,b2);
BesselM1NextCheb(y,-1.03457624656780970260E-8,b0,b1,b2);
BesselM1NextCheb(y,5.74108412545004946722E-8,b0,b1,b2);
BesselM1NextCheb(y,-3.50196060308781257119E-7,b0,b1,b2);
BesselM1NextCheb(y,2.40648494783721712015E-6,b0,b1,b2);
BesselM1NextCheb(y,-1.93619797416608296024E-5,b0,b1,b2);
BesselM1NextCheb(y,1.95215518471351631108E-4,b0,b1,b2);
BesselM1NextCheb(y,-2.85781685962277938680E-3,b0,b1,b2);
BesselM1NextCheb(y,1.03923736576817238437E-1,b0,b1,b2);
BesselM1NextCheb(y,2.72062619048444266945E0,b0,b1,b2);
v=0.5*(b0-b2);
//--- get result
result=MathExp(-x)*v/MathSqrt(x);
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Modified Bessel function, second kind, integer order |
//| Returns modified Bessel function of the second kind |
//| of order n of the argument. |
//| The range is partitioned into the two intervals [0,9.55] and |
//| (9.55, infinity). An ascending power series is used in the |
//| low range, and an asymptotic expansion in the high range. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,30 90000 1.8e-8 3.0e-10 |
//| Error is high only near the crossover point x = 9.55 |
//| between the two expansions used. |
//+------------------------------------------------------------------+
double CBessel::BesselKN(int nn,double x)
{
//--- create variables
double k=0;
double kf=0;
double nk1f=0;
double nkf=0;
double zn=0;
double t=0;
double s=0;
double z0=0;
double z=0;
double ans=0;
double fn=0;
double pn=0;
double pk=0;
double zmn=0;
double tlg=0;
double tox=0;
int i=0;
int n=0;
double eul=0;
//--- initialization
eul=5.772156649015328606065e-1;
//--- check
if(nn<0)
n=-nn;
else
n=nn;
//--- check
if(!CAp::Assert(n<=31,__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(x>0.0,__FUNCTION__+":domain error"))
return(EMPTY_VALUE);
//--- check
if(x<=9.55)
{
//--- initialization
ans=0.0;
z0=0.25*x*x;
fn=1.0;
pn=0.0;
zmn=1.0;
tox=2.0/x;
//--- check
if(n>0)
{
//--- calculation
pn=-eul;
k=1.0;
for(i=1; i<=n-1; i++)
{
pn=pn+1.0/k;
k=k+1.0;
fn=fn*k;
}
zmn=tox;
//--- check
if(n==1)
ans=1.0/x;
else
{
//--- change values
nk1f=fn/n;
kf=1.0;
s=nk1f;
z=-z0;
zn=1.0;
//--- calculation
for(i=1; i<=n-1; i++)
{
nk1f=nk1f/(n-i);
kf=kf*i;
zn=zn*z;
t=nk1f*zn/kf;
s=s+t;
//--- check
if(!CAp::Assert(CMath::m_maxrealnumber-MathAbs(t)>MathAbs(s),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(!(tox>1.0 && CMath::m_maxrealnumber/tox<zmn),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
zmn=zmn*tox;
}
//--- change values
s=s*0.5;
t=MathAbs(s);
//--- check
if(!CAp::Assert(!(zmn>1.0 && CMath::m_maxrealnumber/zmn<t),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(!(t>1.0 && CMath::m_maxrealnumber/t<zmn),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
ans=s*zmn;
}
}
//--- change values
tlg=2.0*MathLog(0.5*x);
pk=-eul;
//--- check
if(n==0)
{
pn=pk;
t=1.0;
}
else
{
pn=pn+1.0/n;
t=1.0/fn;
}
s=(pk+pn-tlg)*t;
k=1.0;
//--- cycle
do
{
t=t*(z0/(k*(k+n)));
pk=pk+1.0/k;
pn=pn+1.0/(k+n);
s=s+(pk+pn-tlg)*t;
k=k+1.0;
}
while(MathAbs(t/s)>CMath::m_machineepsilon);
s=0.5*s/zmn;
//--- check
if(n%2!=0)
{
s=-s;
}
ans=ans+s;
//--- return result
return(ans);
}
//--- check
if(x>MathLog(CMath::m_maxrealnumber))
return(0);
//--- change values
k=n;
pn=4.0*k*k;
pk=1.0;
z0=8.0*x;
fn=1.0;
t=1.0;
s=t;
nkf=CMath::m_maxrealnumber;
i=0;
//--- cycle
do
{
z=pn-pk*pk;
t=t*z/(fn*z0);
nk1f=MathAbs(t);
//--- check
if(i>=n && nk1f>(double)(nkf))
break;
//--- change values
nkf=nk1f;
s=s+t;
fn=fn+1.0;
pk=pk+2.0;
i=i+1;
}
while(MathAbs(t/s)>CMath::m_machineepsilon);
//--- return result
return(MathExp(-x)*MathSqrt(M_PI/(2.0*x))*s);
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselMFirstCheb(const double c,double &b0,double &b1,
double &b2)
{
//--- change values
b0=c;
b1=0.0;
b2=0.0;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselMNextCheb(const double x,const double c,double &b0,
double &b1,double &b2)
{
//--- change values
b2=b1;
b1=b0;
b0=x*b1-b2+c;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselM1FirstCheb(const double c,double &b0,double &b1,
double &b2)
{
//--- change values
b0=c;
b1=0.0;
b2=0.0;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselM1NextCheb(const double x,const double c,double &b0,
double &b1,double &b2)
{
//--- change values
b2=b1;
b1=b0;
b0=x*b1-b2+c;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselAsympt0(const double x,double &pzero,double &qzero)
{
//--- create variables
double xsq=0;
double p2=0;
double q2=0;
double p3=0;
double q3=0;
//--- initialization
pzero=0;
qzero=0;
//--- calculation
xsq=64.0/(x*x);
p2=0.0;
p2=2485.271928957404011288128951+xsq*p2;
p2=153982.6532623911470917825993+xsq*p2;
p2=2016135.283049983642487182349+xsq*p2;
p2=8413041.456550439208464315611+xsq*p2;
p2=12332384.76817638145232406055+xsq*p2;
p2=5393485.083869438325262122897+xsq*p2;
q2=1.0;
q2=2615.700736920839685159081813+xsq*q2;
q2=156001.7276940030940592769933+xsq*q2;
q2=2025066.801570134013891035236+xsq*q2;
q2=8426449.050629797331554404810+xsq*q2;
q2=12338310.22786324960844856182+xsq*q2;
q2=5393485.083869438325560444960+xsq*q2;
p3=-0.0;
p3=-4.887199395841261531199129300+xsq*p3;
p3=-226.2630641933704113967255053+xsq*p3;
p3=-2365.956170779108192723612816+xsq*p3;
p3=-8239.066313485606568803548860+xsq*p3;
p3=-10381.41698748464093880530341+xsq*p3;
p3=-3984.617357595222463506790588+xsq*p3;
q3=1.0;
q3=408.7714673983499223402830260+xsq*q3;
q3=15704.89191515395519392882766+xsq*q3;
q3=156021.3206679291652539287109+xsq*q3;
q3=533291.3634216897168722255057+xsq*q3;
q3=666745.4239319826986004038103+xsq*q3;
q3=255015.5108860942382983170882+xsq*q3;
//--- get result
pzero=p2/q2;
qzero=8*p3/q3/x;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CBessel::BesselAsympt1(const double x,double &pzero,double &qzero)
{
//--- create variables
double xsq=0;
double p2=0;
double q2=0;
double p3=0;
double q3=0;
//--- initialization
pzero=0;
qzero=0;
//--- calculation
xsq=64.0/(x*x);
p2=-1611.616644324610116477412898;
p2=-109824.0554345934672737413139+xsq*p2;
p2=-1523529.351181137383255105722+xsq*p2;
p2=-6603373.248364939109255245434+xsq*p2;
p2=-9942246.505077641195658377899+xsq*p2;
p2=-4435757.816794127857114720794+xsq*p2;
q2=1.0;
q2=-1455.009440190496182453565068+xsq*q2;
q2=-107263.8599110382011903063867+xsq*q2;
q2=-1511809.506634160881644546358+xsq*q2;
q2=-6585339.479723087072826915069+xsq*q2;
q2=-9934124.389934585658967556309+xsq*q2;
q2=-4435757.816794127856828016962+xsq*q2;
p3=35.26513384663603218592175580;
p3=1706.375429020768002061283546+xsq*p3;
p3=18494.26287322386679652009819+xsq*p3;
p3=66178.83658127083517939992166+xsq*p3;
p3=85145.16067533570196555001171+xsq*p3;
p3=33220.91340985722351859704442+xsq*p3;
q3=1.0;
q3=863.8367769604990967475517183+xsq*q3;
q3=37890.22974577220264142952256+xsq*q3;
q3=400294.4358226697511708610813+xsq*q3;
q3=1419460.669603720892855755253+xsq*q3;
q3=1819458.042243997298924553839+xsq*q3;
q3=708712.8194102874357377502472+xsq*q3;
//--- get result
pzero=p2/q2;
qzero=8*p3/q3/x;
}
//+------------------------------------------------------------------+
//| Beta function |
//+------------------------------------------------------------------+
class CBetaF
{
public:
static double Beta(const double a,const double b);
};
//+------------------------------------------------------------------+
//| Beta function |
//| - - |
//| | (a) | (b) |
//| beta( a, b ) = -----------. |
//| - |
//| | (a+b) |
//| For large arguments the logarithm of the function is |
//| evaluated using lgam(), then exponentiated. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,30 30000 8.1e-14 1.1e-14 |
//+------------------------------------------------------------------+
double CBetaF::Beta(const double a,const double b)
{
//--- create variables
double y=0;
double sg=0;
double s=0;
//--- initialization
sg=1;
//--- check
if(!CAp::Assert(a>0.0 || a!=(double)((int)MathFloor(a)),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(b>0.0 || b!=(double)((int)MathFloor(b)),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
y=a+b;
//--- check
if(MathAbs(y)>171.624376956302725)
{
//--- calculation
y=CGammaFunc::LnGamma(y,s);
sg=sg*s;
y=CGammaFunc::LnGamma(b,s)-y;
sg=sg*s;
y=CGammaFunc::LnGamma(a,s)+y;
sg=sg*s;
//--- check
if(!CAp::Assert(y<=(double)(MathLog(CMath::m_maxrealnumber)),__FUNCTION__+": overflow"))
return(EMPTY_VALUE);
//--- return result
return(sg*MathExp(y));
}
//--- function call
y=CGammaFunc::GammaFunc(y);
//--- check
if(!CAp::Assert(y!=0.0,"Overflow in Beta"))
return(EMPTY_VALUE);
//--- check
if(a>b)
{
//--- calculation
y=CGammaFunc::GammaFunc(a)/y;
y=y*CGammaFunc::GammaFunc(b);
}
else
{
//--- calculation
y=CGammaFunc::GammaFunc(b)/y;
y=y*CGammaFunc::GammaFunc(a);
}
//--- return result
return(y);
}
//+------------------------------------------------------------------+
//| Incomplete beta function |
//+------------------------------------------------------------------+
class CIncBetaF
{
public:
static double IncompleteBeta(double a,double b,double x);
static double InvIncompleteBeta(const double a,double b,double y);
private:
static double IncompleteBetaFracExpans(const double a,const double b,const double x,const double big,const double biginv);
static double IncompleteBetaFracExpans2(const double a,const double b,const double x,const double big,const double biginv);
static double IncompleteBetaPowSeries(const double a,const double b,const double x,const double maxgam);
};
//+------------------------------------------------------------------+
//| Incomplete beta integral |
//| Returns incomplete beta integral of the arguments, evaluated |
//| from zero to x. The function is defined as |
//| x |
//| - - |
//| | (a+b) | | a-1 b-1 |
//| ----------- | t (1-t) dt. |
//| - - | | |
//| | (a) | (b) - |
//| 0 |
//| The domain of definition is 0 <= x <= 1. In this |
//| implementation a and b are restricted to positive values. |
//| The integral from x to 1 may be obtained by the symmetry |
//| relation |
//| 1 - incbet( a, b, x ) = incbet( b, a, 1-x ). |
//| The integral is evaluated by a continued fraction expansion |
//| or, when b*x is small, by a power series. |
//| ACCURACY: |
//| Tested at uniformly distributed random points (a,b,x) with a and |
//| b in "domain" and x between 0 and 1. |
//| Relative error |
//| arithmetic domain # trials peak rms |
//| IEEE 0,5 10000 6.9e-15 4.5e-16 |
//| IEEE 0,85 250000 2.2e-13 1.7e-14 |
//| IEEE 0,1000 30000 5.3e-12 6.3e-13 |
//| IEEE 0,10000 250000 9.3e-11 7.1e-12 |
//| IEEE 0,100000 10000 8.7e-10 4.8e-11 |
//| Outputs smaller than the IEEE gradual underflow threshold |
//| were excluded from these statistics. |
//+------------------------------------------------------------------+
double CIncBetaF::IncompleteBeta(double a,double b,double x)
{
//--- check
if(!CAp::Assert(a>0 && b>0,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(x>=0 && x<=1,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- check
if(x==0)
return(0);
//--- check
if(x==1)
return(1);
//--- create variables
double t;
double xc;
double w;
double y;
int flag=0;
double sg=0;
double big=4.503599627370496e15;
double biginv=2.22044604925031308085e-16;
double maxgam=171.624376956302725;
double minlog=MathLog(CMath::m_minrealnumber);
double maxlog=MathLog(CMath::m_maxrealnumber);
//--- check
if(b*x<=1.0 && x<=0.95)
return(IncompleteBetaPowSeries(a,b,x,maxgam));
//--- change value
w=1.0-x;
//--- check
if(x>a/(a+b))
{
flag=1;
t=a;
a=b;
b=t;
xc=x;
x=w;
}
else
xc=w;
//--- check
if(flag==1 && b*x<=1.0 && x<=0.95)
{
t=IncompleteBetaPowSeries(a,b,x,maxgam);
//--- check
if(t<=CMath::m_machineepsilon)
return(1.0-CMath::m_machineepsilon);
else
return(1.0-t);
}
//--- change value
y=x*(a+b-2.0)-(a-1.0);
//--- check
if(y<0.0)
w=IncompleteBetaFracExpans(a,b,x,big,biginv);
else
w=IncompleteBetaFracExpans2(a,b,x,big,biginv)/xc;
//--- change values
y=a*MathLog(x);
t=b*MathLog(xc);
//--- check
if(a+b<maxgam && MathAbs(y)<maxlog && MathAbs(t)<maxlog)
{
t=MathPow(xc,b);
t*=MathPow(x,a);
t/=a;
t*=w;
t*=CGammaFunc::GammaFunc(a+b)/(CGammaFunc::GammaFunc(a)*CGammaFunc::GammaFunc(b));
//--- check
if(flag==1)
{
//--- check
if(t<=CMath::m_machineepsilon)
return(1.0-CMath::m_machineepsilon);
else
return(1.0-t);
}
else
return(t);
}
//--- change values
y+=t+CGammaFunc::LnGamma(a+b,sg)-CGammaFunc::LnGamma(a,sg)-CGammaFunc::LnGamma(b,sg);
y+=MathLog(w/a);
//--- check
if(y<minlog)
t=0.0;
else
t=MathExp(y);
//--- check
if(flag==1)
{
if(t<=CMath::m_machineepsilon)
t=1.0-CMath::m_machineepsilon;
else
t=1.0-t;
}
//--- return result
return(t);
}
//+------------------------------------------------------------------+
//| Inverse of imcomplete beta integral |
//| Given y, the function finds x such that |
//| incbet( a, b, x ) = y . |
//| The routine performs interval halving or Newton iterations to |
//| find the root of incbet(a,b,x) - y = 0. |
//| ACCURACY: |
//| Relative error: |
//| x a,b |
//| arithmetic domain domain # trials peak rms |
//| IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13 |
//| IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15 |
//| IEEE 0,1 0,5 50000 1.1e-12 5.5e-15 |
//| With a and b constrained to half-integer or integer values: |
//| IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13 |
//| IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16 |
//| With a = .5, b constrained to half-integer or integer values: |
//| IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11 |
//+------------------------------------------------------------------+
double CIncBetaF::InvIncompleteBeta(const double a,double b,double y)
{
//--- check
if(!CAp::Assert(y>=0 && y<=1,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- check
if(y==0)
return(0);
//--- check
if(y==1.0)
return(1);
//--- create variables
double aaa=0;
double bbb=0;
double y0=0;
double d=0;
double yyy=0;
double x=0;
double x0=0.0;
double x1=1.0;
double lgm=0;
double yp=0;
double di=0;
double dithresh=0;
double yl=0.0;
double yh=1.0;
double xt=0;
int i=0;
int rflg=0;
int dir=0;
int nflg=0;
double s=0;
int mainlooppos=0;
int ihalve=1;
int ihalvecycle=2;
int newt=3;
int newtcycle=4;
int breaknewtcycle=5;
int breakihalvecycle=6;
//--- main
while(true)
{
//--- start
if(mainlooppos==0)
{
if(a<=1.0 || b<=1.0)
{
dithresh=1.0e-6;
rflg=0;
aaa=a;
bbb=b;
y0=y;
x=aaa/(aaa+bbb);
yyy=IncompleteBeta(aaa,bbb,x);
mainlooppos=ihalve;
continue;
}
else
dithresh=1.0e-4;
//--- change value
yp=-CNormalDistr::InvNormalDistribution(y);
//--- check
if(y>0.5)
{
rflg=1;
aaa=b;
bbb=a;
y0=1.0-y;
yp=-yp;
}
else
{
rflg=0;
aaa=a;
bbb=b;
y0=y;
}
//--- change values
lgm=(yp*yp-3.0)/6.0;
x=2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
d=yp*MathSqrt(x+lgm)/x-(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))*(lgm+5.0/6.0-2.0/(3.0*x));
d=2.0*d;
//--- check
if(d<MathLog(CMath::m_minrealnumber))
{
x=0;
break;
}
//--- change values
x=aaa/(aaa+bbb*MathExp(d));
yyy=IncompleteBeta(aaa,bbb,x);
yp=(yyy-y0)/y0;
//--- check
if(MathAbs(yp)<0.2)
{
mainlooppos=newt;
continue;
}
mainlooppos=ihalve;
continue;
}
//--- check
if(mainlooppos==ihalve)
{
dir=0;
di=0.5;
i=0;
mainlooppos=ihalvecycle;
continue;
}
//--- check
if(mainlooppos==ihalvecycle)
{
//--- check
if(i<=99)
{
//--- check
if(i!=0)
{
x=x0+di*(x1-x0);
//--- check
if(x==1.0)
x=1.0-CMath::m_machineepsilon;
//--- check
if(x==0.0)
{
di=0.5;
x=x0+di*(x1-x0);
//--- check
if(x==0.0)
break;
}
//--- change values
yyy=IncompleteBeta(aaa,bbb,x);
yp=(x1-x0)/(x1+x0);
//--- check
if(MathAbs(yp)<dithresh)
{
mainlooppos=newt;
continue;
}
yp=(yyy-y0)/y0;
//--- check
if(MathAbs(yp)<dithresh)
{
mainlooppos=newt;
continue;
}
}
//--- check
if(yyy<y0)
{
x0=x;
yl=yyy;
//--- check
if(dir<0)
{
dir=0;
di=0.5;
}
else
{
//--- check
if(dir>3)
di=1.0-(1.0-di)*(1.0-di);
else
{
//--- check
if(dir>1)
di=0.5*di+0.5;
else
di=(y0-yyy)/(yh-yl);
}
}
dir++;
//--- check
if(x0>0.75)
{
//--- check
if(rflg==1)
{
rflg=0;
aaa=a;
bbb=b;
y0=y;
}
else
{
rflg=1;
aaa=b;
bbb=a;
y0=1.0-y;
}
//--- change values
x=1.0-x;
yyy=IncompleteBeta(aaa,bbb,x);
x0=0.0;
yl=0.0;
x1=1.0;
yh=1.0;
mainlooppos=ihalve;
continue;
}
}
else
{
x1=x;
//--- check
if(rflg==1 && x1<CMath::m_machineepsilon)
{
x=0.0;
break;
}
yh=yyy;
//--- check
if(dir>0)
{
dir=0;
di=0.5;
}
else
{
//--- check
if(dir<-3)
di*=di;
else
{
//--- check
if(dir<-1)
di*=0.5;
else
di=(yyy-y0)/(yh-yl);
}
}
dir--;
}
i++;
mainlooppos=ihalvecycle;
continue;
}
else
{
mainlooppos=breakihalvecycle;
continue;
}
}
//--- check
if(mainlooppos==breakihalvecycle)
{
//--- check
if(x0>=1.0)
{
x=1.0-CMath::m_machineepsilon;
break;
}
//--- check
if(x<=0.0)
{
x=0.0;
break;
}
mainlooppos=newt;
continue;
}
//--- check
if(mainlooppos==newt)
{
if(nflg!=0)
break;
//--- change values
nflg=1;
lgm=CGammaFunc::LnGamma(aaa+bbb,s)-CGammaFunc::LnGamma(aaa,s)-CGammaFunc::LnGamma(bbb,s);
i=0;
mainlooppos=newtcycle;
continue;
}
//--- check
if(mainlooppos==newtcycle)
{
//--- check
if(i<=7)
{
//--- check
if(i!=0)
yyy=IncompleteBeta(aaa,bbb,x);
//--- check
if(yyy<yl)
{
x=x0;
yyy=yl;
}
else
{
//--- check
if(yyy>yh)
{
x=x1;
yyy=yh;
}
else
{
//--- check
if(yyy<y0)
{
x0=x;
yl=yyy;
}
else
{
x1=x;
yh=yyy;
}
}
}
//--- check
if(x==1.0 || x==0.0)
{
mainlooppos=breaknewtcycle;
continue;
}
//--- change value
d=(aaa-1.0)*MathLog(x)+(bbb-1.0)*MathLog(1.0-x)+lgm;
//--- check
if(d<MathLog(CMath::m_minrealnumber))
break;
//--- check
if(d>MathLog(CMath::m_maxrealnumber))
{
mainlooppos=breaknewtcycle;
continue;
}
//--- change values
d=MathExp(d);
d=(yyy-y0)/d;
xt=x-d;
//--- check
if(xt<=x0)
{
yyy=(x-x0)/(x1-x0);
xt=x0+0.5*yyy*(x-x0);
//--- check
if(xt<=0.0)
{
mainlooppos=breaknewtcycle;
continue;
}
}
//--- check
if(xt>=x1)
{
yyy=(x1-x)/(x1-x0);
xt=x1-0.5*yyy*(x1-x);
//--- check
if(xt>=1.0)
{
mainlooppos=breaknewtcycle;
continue;
}
}
x=xt;
//--- check
if(MathAbs(d/x)<128.0*CMath::m_machineepsilon)
break;
//--- change values
i=i+1;
mainlooppos=newtcycle;
continue;
}
else
{
mainlooppos=breaknewtcycle;
continue;
}
}
//--- check
if(mainlooppos==breaknewtcycle)
{
dithresh=256.0*CMath::m_machineepsilon;
mainlooppos=ihalve;
continue;
}
}
//--- get result
if(rflg!=0)
{
if(x<=CMath::m_machineepsilon)
x=1.0-CMath::m_machineepsilon;
else
x=1.0-x;
}
//--- return result
return(x);
}
//+------------------------------------------------------------------+
//| Continued fraction expansion #1 for incomplete beta integral |
//+------------------------------------------------------------------+
double CIncBetaF::IncompleteBetaFracExpans(const double a,const double b,
const double x,const double big,
const double biginv)
{
//--- create variables
double xk=0;
double pk=0;
double pkm1=1.0;
double pkm2=0.0;
double qk=0;
double qkm1=1.0;
double qkm2=1.0;
double k1=a;
double k2=a+b;
double k3=a;
double k4=a+1.0;
double k5=1.0;
double k6=b-1.0;
double k7=k4;
double k8=a+2.0;
double r=1.0;
double t=0;
double ans=1.0;
double thresh=3.0*CMath::m_machineepsilon;
int n=0;
//--- cycle
do
{
xk=-(x*k1*k2/(k3*k4));
pk=pkm1+pkm2*xk;
qk=qkm1+qkm2*xk;
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
xk=x*k5*k6/(k7*k8);
pk=pkm1+pkm2*xk;
qk=qkm1+qkm2*xk;
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
//--- check
if(qk!=0)
r=pk/qk;
//--- check
if(r!=0)
{
t=MathAbs((ans-r)/r);
ans=r;
}
else
t=1.0;
//--- check
if(t<thresh)
break;
//--- change values
k1+=1.0;
k2+=1.0;
k3+=2.0;
k4+=2.0;
k5+=1.0;
k6-=1.0;
k7+=2.0;
k8+=2.0;
//--- check
if(MathAbs(qk)+MathAbs(pk)>big)
{
pkm2*=biginv;
pkm1*=biginv;
qkm2*=biginv;
qkm1*=biginv;
}
//--- check
if(MathAbs(qk)<biginv || MathAbs(pk)<biginv)
{
pkm2*=big;
pkm1*=big;
qkm2*=big;
qkm1*=big;
}
n++;
}
while(n!=300);
//--- return result
return(ans);
}
//+------------------------------------------------------------------+
//| Continued fraction expansion #2 |
//| for incomplete beta integral |
//+------------------------------------------------------------------+
double CIncBetaF::IncompleteBetaFracExpans2(const double a,const double b,
const double x,const double big,
const double biginv)
{
//--- create variables
double xk=0;
double pk=0;
double pkm1=1.0;
double pkm2=0.0;
double qk=0;
double qkm1=1.0;
double qkm2=1.0;
double k1=a;
double k2=b-1.0;
double k3=a;
double k4=a+1.0;
double k5=1.0;
double k6=a+b;
double k7=a+1.0;
double k8=a+2.0;
double r=1.0;
double t=0;
double ans=1.0;
double z=x/(1.0-x);
double thresh=3.0*CMath::m_machineepsilon;
int n=0;
//--- cycle
do
{
xk=-(z*k1*k2/(k3*k4));
pk=pkm1+pkm2*xk;
qk=qkm1+qkm2*xk;
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
xk=z*k5*k6/(k7*k8);
pk=pkm1+pkm2*xk;
qk=qkm1+qkm2*xk;
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
//--- check
if(qk!=0)
r=pk/qk;
//--- check
if(r!=0)
{
t=MathAbs((ans-r)/r);
ans=r;
}
else
t=1.0;
//--- check
if(t<thresh)
break;
//--- change values
k1+=1.0;
k2-=1.0;
k3+=2.0;
k4+=2.0;
k5+=1.0;
k6+=1.0;
k7+=2.0;
k8+=2.0;
//--- check
if(MathAbs(qk)+MathAbs(pk)>big)
{
pkm2*=biginv;
pkm1*=biginv;
qkm2*=biginv;
qkm1*=biginv;
}
//--- check
if(MathAbs(qk)<biginv || MathAbs(pk)<biginv)
{
pkm2*=big;
pkm1*=big;
qkm2*=big;
qkm1*=big;
}
n++;
}
while(n!=300);
//--- return result
return(ans);
}
//+------------------------------------------------------------------+
//| Power series for incomplete beta integral. |
//| Use when b*x is small and x not too close to 1. |
//+------------------------------------------------------------------+
double CIncBetaF::IncompleteBetaPowSeries(const double a,const double b,
const double x,const double maxgam)
{
//--- create variables
double s=0.0;
double u=(1.0-b)*x;
double t=u;
double v=u/(a+1.0);
double n=2.0;
double t1=v;
double ai=1.0/a;
double z=CMath::m_machineepsilon*ai;
double sg=0;
//--- cycle
while(MathAbs(v)>z)
{
u=(n-b)*x/n;
t=t*u;
v=t/(a+n);
s=s+v;
n=n+1.0;
}
//--- change values
s+=t1;
s+=ai;
u=a*MathLog(x);
//--- check
if(a+b<maxgam && MathAbs(u)<MathLog(CMath::m_maxrealnumber))
{
t=CGammaFunc::GammaFunc(a+b)/(CGammaFunc::GammaFunc(a)*CGammaFunc::GammaFunc(b));
s*=t*MathPow(x,a);
}
else
{
t=CGammaFunc::LnGamma(a+b,sg)-CGammaFunc::LnGamma(a,sg)-CGammaFunc::LnGamma(b,sg)+u+MathLog(s);
//--- check
if(t<MathLog(CMath::m_minrealnumber))
s=0.0;
else
s=MathExp(t);
}
//--- return result
return(s);
}
//+------------------------------------------------------------------+
//| Binomial distribution |
//+------------------------------------------------------------------+
class CBinomialDistr
{
public:
static double BinomialDistribution(const int k,const int n,const double p);
static double BinomialComplDistribution(const int k,const int n,const double p);
static double InvBinomialDistribution(const int k,const int n,const double y);
};
//+------------------------------------------------------------------+
//| Binomial distribution |
//| Returns the sum of the terms 0 through k of the Binomial |
//| probability density: |
//| k |
//| -- ( n ) j n-j |
//| > ( ) p (1-p) |
//| -- ( j ) |
//| j=0 |
//| The terms are not summed directly; instead the incomplete |
//| beta integral is employed, according to the formula |
//| y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ). |
//| The arguments must be positive, with p ranging from 0 to 1. |
//| ACCURACY: |
//| Tested at random points (a,b,p), with p between 0 and 1. |
//| a,b Relative error: |
//| arithmetic domain # trials peak rms |
//| For p between 0.001 and 1: |
//| IEEE 0,100 100000 4.3e-15 2.6e-16 |
//+------------------------------------------------------------------+
double CBinomialDistr::BinomialDistribution(const int k,const int n,
const double p)
{
//--- check
if(!CAp::Assert(p>=0 && p<=1,__FUNCTION__+": the eroor variable"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(k>=-1 && k<=n,__FUNCTION__+": the eroor variable"))
return(EMPTY_VALUE);
//--- check
if(k==-1)
return(0);
//--- check
if(k==n)
return(1);
//--- create variables
double dk;
double dn=n-k;
//--- check
if(k==0)
dk=MathPow(1.0-p,dn);
else
{
dk=k+1;
dk=CIncBetaF::IncompleteBeta(dn,dk,1.0-p);
}
//--- return result
return(dk);
}
//+------------------------------------------------------------------+
//| Complemented binomial distribution |
//| Returns the sum of the terms k+1 through n of the Binomial |
//| probability density: |
//| n |
//| -- ( n ) j n-j |
//| > ( ) p (1-p) |
//| -- ( j ) |
//| j=k+1 |
//| The terms are not summed directly; instead the incomplete |
//| beta integral is employed, according to the formula |
//| y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ). |
//| The arguments must be positive, with p ranging from 0 to 1. |
//| ACCURACY: |
//| Tested at random points (a,b,p). |
//| a,b Relative error: |
//| arithmetic domain # trials peak rms |
//| For p between 0.001 and 1: |
//| IEEE 0,100 100000 6.7e-15 8.2e-16 |
//| For p between 0 and .001: |
//| IEEE 0,100 100000 1.5e-13 2.7e-15 |
//+------------------------------------------------------------------+
double CBinomialDistr::BinomialComplDistribution(const int k,
const int n,
const double p)
{
//--- check
if(!CAp::Assert(p>=0 && p<=1,__FUNCTION__+": the eroor variable"))
return(EMPTY_VALUE);
//--- check
if(!CAp::Assert(k>=-1 && k<=n,__FUNCTION__+": the eroor variable"))
return(EMPTY_VALUE);
//--- check
if(k==-1)
return(1.0);
//--- check
if(k==n)
return(0);
//--- create variables
double dk;
double dn=n-k;
//--- check
if(k==0)
{
if(p<0.01)
dk=-CNearUnitYUnit::NUExp1m(dn*CNearUnitYUnit::NULog1p(-p));
else
dk=1.0-MathPow(1.0-p,dn);
}
else
{
dk=k+1;
dk=CIncBetaF::IncompleteBeta(dk,dn,p);
}
//--- return result
return(dk);
}
//+------------------------------------------------------------------+
//| Inverse binomial distribution |
//| Finds the event probability p such that the sum of the |
//| terms 0 through k of the Binomial probability density |
//| is equal to the given cumulative probability y. |
//| This is accomplished using the inverse beta integral |
//| function and the relation |
//| 1 - p = incbi( n-k, k+1, y ). |
//| ACCURACY: |
//| Tested at random points (a,b,p). |
//| a,b Relative error: |
//| arithmetic domain # trials peak rms |
//| For p between 0.001 and 1: |
//| IEEE 0,100 100000 2.3e-14 6.4e-16 |
//| IEEE 0,10000 100000 6.6e-12 1.2e-13 |
//| For p between 10^-6 and 0.001: |
//| IEEE 0,100 100000 2.0e-12 1.3e-14 |
//| IEEE 0,10000 100000 1.5e-12 3.2e-14 |
//+------------------------------------------------------------------+
double CBinomialDistr::InvBinomialDistribution(const int k,const int n,
const double y)
{
//--- create variables
double dk=0;
double dn=0;
double p=0;
//--- check
if(!CAp::Assert(k>=0 && k<n,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- initialization
dn=n-k;
//--- check
if(k==0)
{
//--- check
if(y>0.8)
p=-CNearUnitYUnit::NUExp1m(CNearUnitYUnit::NULog1p(y-1.0)/dn);
else
p=1.0-MathPow(y,1.0/dn);
}
else
{
dk=k+1;
//--- function call
p=CIncBetaF::IncompleteBeta(dn,dk,0.5);
//--- check
if(p>0.5)
p=CIncBetaF::InvIncompleteBeta(dk,dn,1.0-y);
else
p=1.0-CIncBetaF::InvIncompleteBeta(dn,dk,y);
}
//--- return result
return(p);
}
//+------------------------------------------------------------------+
//| Chebyshev polynomials |
//+------------------------------------------------------------------+
class CChebyshev
{
public:
static double ChebyshevCalculate(const int r,const int n,const double x);
static double ChebyshevSum(double &c[],const int r,const int n,const double x);
static void ChebyshevCoefficients(const int n,double &c[]);
static void FromChebyshev(double &a[],const int n,double &b[]);
};
//+------------------------------------------------------------------+
//| Calculation of the value of the Chebyshev polynomials of the |
//| first and second kinds. |
//| Parameters: |
//| r - polynomial kind, either 1 or 2. |
//| n - degree, n>=0 |
//| x - argument, -1 <= x <= 1 |
//| Result: |
//| the value of the Chebyshev polynomial at x |
//+------------------------------------------------------------------+
double CChebyshev::ChebyshevCalculate(const int r,const int n,
const double x)
{
//--- create variables
double result=0;
int i=0;
double a=0;
double b=0;
//--- Prepare A and B
if(r==1)
{
a=1;
b=x;
}
else
{
a=1;
b=2*x;
}
//--- Special cases: N=0 or N=1
if(n==0)
return(a);
//--- check
if(n==1)
return(b);
//--- General case: N>=2
for(i=2; i<=n; i++)
{
result=2*x*b-a;
a=b;
b=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Summation of Chebyshev polynomials using Clenshaw?s recurrence |
//| formula. |
//| This routine calculates |
//| c[0]*T0(x) + c[1]*T1(x) + ... + c[N]*TN(x) |
//| or |
//| c[0]*U0(x) + c[1]*U1(x) + ... + c[N]*UN(x) |
//| depending on the R. |
//| Parameters: |
//| r - polynomial kind, either 1 or 2. |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Chebyshev polynomial at x |
//+------------------------------------------------------------------+
double CChebyshev::ChebyshevSum(double &c[],const int r,const int n,
const double x)
{
//--- create variables
double result=0;
double b1=0;
double b2=0;
int i=0;
//--- initialization
b1=0;
b2=0;
//--- calculation
for(i=n; i>=1; i--)
{
result=2*x*b1-b2+c[i];
b2=b1;
b1=result;
}
//--- check
if(r==1)
result=-b2+x*b1+c[0];
else
result=-b2+2*x*b1+c[0];
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Representation of Tn as C[0] + C[1]*X + ... + C[N]*X^N |
//| Input parameters: |
//| N - polynomial degree, n>=0 |
//| Output parameters: |
//| C - coefficients |
//+------------------------------------------------------------------+
void CChebyshev::ChebyshevCoefficients(const int n,double &c[])
{
//--- allocation
ArrayResize(c,n+1);
//--- initialization
ArrayInitialize(c,0.0);
//--- check
if(n==0 || n==1)
c[n]=1;
else
{
c[n]=MathExp((n-1)*MathLog(2));
//--- calculation
for(int i=0; i<=n/2-1; i++)
c[n-2*(i+1)]=-(c[n-2*i]*(n-2*i)*(n-2*i-1)/4/(i+1)/(n-i-1));
}
}
//+------------------------------------------------------------------+
//| Conversion of a series of Chebyshev polynomials to a power |
//| series. |
//| Represents A[0]*T0(x) + A[1]*T1(x) + ... + A[N]*Tn(x) as |
//| B[0] + B[1]*X + ... + B[N]*X^N. |
//| Input parameters: |
//| A - Chebyshev series coefficients |
//| N - degree, N>=0 |
//| Output parameters |
//| B - power series coefficients |
//+------------------------------------------------------------------+
void CChebyshev::FromChebyshev(double &a[],const int n,double &b[])
{
//--- create variables
int i=0;
int k=0;
double e=0;
double d=0;
//--- allocation
ArrayResize(b,n+1);
//--- initialization
ArrayInitialize(b,0.0);
//--- calculation
d=0;
i=0;
do
{
k=i;
//--- cycle
do
{
e=b[k];
b[k]=0;
//--- check
if(i<=1 && k==i)
b[k]=1;
else
{
//--- check
if(i!=0)
b[k]=2*d;
//--- check
if(k>i+1)
b[k]=b[k]-b[k-2];
}
//--- change values
d=e;
k=k+1;
}
while(k<=n);
//--- change values
d=b[i];
e=0;
k=i;
//--- calculation
while(k<=n)
{
e=e+b[k]*a[k];
k=k+2;
}
//--- change values
b[i]=e;
i=i+1;
}
while(i<=n);
}
//+------------------------------------------------------------------+
//| Chi-square distribution |
//+------------------------------------------------------------------+
class CChiSquareDistr
{
public:
static double ChiSquareDistribution(const double v,const double x);
static double ChiSquareComplDistribution(const double v,const double x);
static double InvChiSquareDistribution(const double v,const double y);
};
//+------------------------------------------------------------------+
//| Chi-square distribution |
//| Returns the area under the left hand tail (from 0 to x) |
//| of the Chi square probability density function with |
//| v degrees of freedom. |
//| x |
//| - |
//| 1 | | v/2-1 -t/2 |
//| P( x | v ) = ----------- | t e dt |
//| v/2 - | | |
//| 2 | (v/2) - |
//| 0 |
//| where x is the Chi-square variable. |
//| The incomplete gamma integral is used, according to the |
//| formula |
//| y = chdtr( v, x ) = igam( v/2.0, x/2.0 ). |
//| The arguments must both be positive. |
//| ACCURACY: |
//| See incomplete gamma function |
//+------------------------------------------------------------------+
double CChiSquareDistr::ChiSquareDistribution(const double v,
const double x)
{
//--- check
if(!CAp::Assert(x>=0 && v>=1,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- return result
return(CIncGammaF::IncompleteGamma(v/2.0,x/2.0));
}
//+------------------------------------------------------------------+
//| Complemented Chi-square distribution |
//| Returns the area under the right hand tail (from x to |
//| infinity) of the Chi square probability density function |
//| with v degrees of freedom: |
//| inf. |
//| - |
//| 1 | | v/2-1 -t/2 |
//| P( x | v ) = ----------- | t e dt |
//| v/2 - | | |
//| 2 | (v/2) - |
//| x |
//| where x is the Chi-square variable. |
//| The incomplete gamma integral is used, according to the |
//| formula |
//| y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ). |
//| The arguments must both be positive. |
//| ACCURACY: |
//| See incomplete gamma function |
//+------------------------------------------------------------------+
double CChiSquareDistr::ChiSquareComplDistribution(const double v,
const double x)
{
//--- check
if(!CAp::Assert(x>=0.0 && v>=1.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- return result
return(CIncGammaF::IncompleteGammaC(v/2.0,x/2.0));
}
//+------------------------------------------------------------------+
//| Inverse of complemented Chi-square distribution |
//| Finds the Chi-square argument x such that the integral |
//| from x to infinity of the Chi-square density is equal |
//| to the given cumulative probability y. |
//| This is accomplished using the inverse gamma integral |
//| function and the relation |
//| x/2 = igami( df/2, y ); |
//| ACCURACY: |
//| See inverse incomplete gamma function |
//+------------------------------------------------------------------+
double CChiSquareDistr::InvChiSquareDistribution(const double v,
const double y)
{
//--- check
if(!CAp::Assert((y>=0.0 && y<=1.0) && v>=1.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- return result
return(2*CIncGammaF::InvIncompleteGammaC(0.5*v,y));
}
//+------------------------------------------------------------------+
//| Dawson integral |
//+------------------------------------------------------------------+
class CDawson
{
public:
static double DawsonIntegral(double x);
};
//+------------------------------------------------------------------+
//| Dawson's Integral |
//| Approximates the integral |
//| x |
//| - |
//| 2 | | 2 |
//| dawsn(x) = exp( -x ) | exp( t ) dt |
//| | | |
//| - |
//| 0 |
//| Three different rational approximations are employed, for |
//| the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,10 10000 6.9e-16 1.0e-16 |
//+------------------------------------------------------------------+
double CDawson::DawsonIntegral(double x)
{
//--- create variables
int sg=1;
double x2=0;
double y=0;
double an=0;
double ad=0;
double bn=0;
double bd=0;
double cn=0;
double cd=0;
//--- check
if(x<0.0)
{
sg=-1;
x=-x;
}
//--- check
if(x<3.25)
{
//--- calculation
x2=x*x;
an=1.13681498971755972054E-11;
an=an*x2+8.49262267667473811108E-10;
an=an*x2+1.94434204175553054283E-8;
an=an*x2+9.53151741254484363489E-7;
an=an*x2+3.07828309874913200438E-6;
an=an*x2+3.52513368520288738649E-4;
an=an*x2+-8.50149846724410912031E-4;
an=an*x2+4.22618223005546594270E-2;
an=an*x2+-9.17480371773452345351E-2;
an=an*x2+9.99999999999999994612E-1;
ad=2.40372073066762605484E-11;
ad=ad*x2+1.48864681368493396752E-9;
ad=ad*x2+5.21265281010541664570E-8;
ad=ad*x2+1.27258478273186970203E-6;
ad=ad*x2+2.32490249820789513991E-5;
ad=ad*x2+3.25524741826057911661E-4;
ad=ad*x2+3.48805814657162590916E-3;
ad=ad*x2+2.79448531198828973716E-2;
ad=ad*x2+1.58874241960120565368E-1;
ad=ad*x2+5.74918629489320327824E-1;
ad=ad*x2+1.00000000000000000539E0;
y=x*an/ad;
//--- return result
return(sg*y);
}
x2=1.0/(x*x);
//--- check
if(x<6.25)
{
//--- calculation
bn=5.08955156417900903354E-1;
bn=bn*x2-2.44754418142697847934E-1;
bn=bn*x2+9.41512335303534411857E-2;
bn=bn*x2-2.18711255142039025206E-2;
bn=bn*x2+3.66207612329569181322E-3;
bn=bn*x2-4.23209114460388756528E-4;
bn=bn*x2+3.59641304793896631888E-5;
bn=bn*x2-2.14640351719968974225E-6;
bn=bn*x2+9.10010780076391431042E-8;
bn=bn*x2-2.40274520828250956942E-9;
bn=bn*x2+3.59233385440928410398E-11;
bd=1.00000000000000000000E0;
bd=bd*x2-6.31839869873368190192E-1;
bd=bd*x2+2.36706788228248691528E-1;
bd=bd*x2-5.31806367003223277662E-2;
bd=bd*x2+8.48041718586295374409E-3;
bd=bd*x2-9.47996768486665330168E-4;
bd=bd*x2+7.81025592944552338085E-5;
bd=bd*x2-4.55875153252442634831E-6;
bd=bd*x2+1.89100358111421846170E-7;
bd=bd*x2-4.91324691331920606875E-9;
bd=bd*x2+7.18466403235734541950E-11;
y=1.0/x+x2*bn/(bd*x);
//--- return result
return(sg*0.5*y);
}
//--- check
if(x>1.0E9)
return(sg*0.5/x);
//--- calculation
cn=-5.90592860534773254987E-1;
cn=cn*x2+6.29235242724368800674E-1;
cn=cn*x2-1.72858975380388136411E-1;
cn=cn*x2+1.64837047825189632310E-2;
cn=cn*x2-4.86827613020462700845E-4;
cd=1.00000000000000000000E0;
cd=cd*x2-2.69820057197544900361E0;
cd=cd*x2+1.73270799045947845857E0;
cd=cd*x2-3.93708582281939493482E-1;
cd=cd*x2+3.44278924041233391079E-2;
cd=cd*x2-9.73655226040941223894E-4;
y=1.0/x+x2*cn/(cd*x);
//--- return result
return(sg*0.5*y);
}
//+------------------------------------------------------------------+
//| Elliptic integral |
//+------------------------------------------------------------------+
class CElliptic
{
public:
static double EllipticIntegralK(const double m);
static double EllipticIntegralKhighPrecision(const double m1);
static double IncompleteEllipticIntegralK(double phi,const double m);
static double EllipticIntegralE(double m);
static double IncompleteEllipticIntegralE(const double phi,const double m);
};
//+------------------------------------------------------------------+
//| Complete elliptic integral of the first kind |
//| Approximates the integral |
//| pi/2 |
//| - |
//| | | |
//| | dt |
//| K(m) = | ------------------ |
//| | 2 |
//| | | sqrt( 1 - m sin t ) |
//| - |
//| 0 |
//| using the approximation |
//| P(x) - log x Q(x). |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,1 30000 2.5e-16 6.8e-17 |
//+------------------------------------------------------------------+
double CElliptic::EllipticIntegralK(const double m)
{
//--- return result
return(EllipticIntegralKhighPrecision(1.0-m));
}
//+------------------------------------------------------------------+
//| Complete elliptic integral of the first kind |
//| Approximates the integral |
//| pi/2 |
//| - |
//| | | |
//| | dt |
//| K(m) = | ------------------ |
//| | 2 |
//| | | sqrt( 1 - m sin t ) |
//| - |
//| 0 |
//| where m = 1 - m1, using the approximation |
//| P(x) - log x Q(x). |
//| The argument m1 is used rather than m so that the logarithmic |
//| singularity at m = 1 will be shifted to the origin; this |
//| preserves maximum accuracy. |
//| K(0) = pi/2. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,1 30000 2.5e-16 6.8e-17 |
//+------------------------------------------------------------------+
double CElliptic::EllipticIntegralKhighPrecision(const double m1)
{
//--- create variables
double result=0;
double p=0;
double q=0;
//--- check
if(m1<=CMath::m_machineepsilon)
result=1.3862943611198906188E0-0.5*MathLog(m1);
else
{
//--- calculation
p=1.37982864606273237150E-4;
p=p*m1+2.28025724005875567385E-3;
p=p*m1+7.97404013220415179367E-3;
p=p*m1+9.85821379021226008714E-3;
p=p*m1+6.87489687449949877925E-3;
p=p*m1+6.18901033637687613229E-3;
p=p*m1+8.79078273952743772254E-3;
p=p*m1+1.49380448916805252718E-2;
p=p*m1+3.08851465246711995998E-2;
p=p*m1+9.65735902811690126535E-2;
p=p*m1+1.38629436111989062502E0;
q=2.94078955048598507511E-5;
q=q*m1+9.14184723865917226571E-4;
q=q*m1+5.94058303753167793257E-3;
q=q*m1+1.54850516649762399335E-2;
q=q*m1+2.39089602715924892727E-2;
q=q*m1+3.01204715227604046988E-2;
q=q*m1+3.73774314173823228969E-2;
q=q*m1+4.88280347570998239232E-2;
q=q*m1+7.03124996963957469739E-2;
q=q*m1+1.24999999999870820058E-1;
q=q*m1+4.99999999999999999821E-1;
//--- get result
result=p-q*MathLog(m1);
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Incomplete elliptic integral of the first kind F(phi|m) |
//| Approximates the integral |
//| phi |
//| - |
//| | | |
//| | dt |
//| F(phi_\m) = | ------------------ |
//| | 2 |
//| | | sqrt( 1 - m sin t ) |
//| - |
//| 0 |
//| of amplitude phi and modulus m, using the arithmetic - |
//| geometric mean algorithm. |
//| ACCURACY: |
//| Tested at random points with m in [0, 1] and phi as indicated. |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -10,10 200000 7.4e-16 1.0e-16 |
//+------------------------------------------------------------------+
double CElliptic::IncompleteEllipticIntegralK(double phi,const double m)
{
//--- create variables
double pio2=1.57079632679489661923;
double a=0;
double b=0;
double c=0;
double e=0;
double temp=0;
double t=0;
double k=0;
int d=0;
int md=0;
int s=0;
int npio2=0;
//--- check
if(m==0.0)
return(phi);
//--- initialization
a=1.0-m;
//--- check
if(a==0.0)
return(MathLog(MathTan(0.5*(pio2+phi))));
//--- initialization
npio2=(int)MathFloor(phi/pio2);
//--- check
if(npio2%2!=0)
npio2=npio2+1;
//--- check
if(npio2!=0)
{
k=EllipticIntegralK(1-a);
phi=phi-npio2*pio2;
}
else
k=0;
//--- check
if(phi<0.0)
{
phi=-phi;
s=-1;
}
else
s=0;
//--- calculation
b=MathSqrt(a);
t=MathTan(phi);
//--- check
if(MathAbs(t)>10)
{
e=1.0/(b*t);
//--- check
if(MathAbs(e)<10)
{
e=MathArctan(e);
//--- check
if(npio2==0)
k=EllipticIntegralK(1-a);
//--- function call
temp=k-IncompleteEllipticIntegralK(e,m);
//--- check
if(s<0)
temp=-temp;
//--- return result
return(temp+npio2*k);
}
}
//--- change values
a=1.0;
c=MathSqrt(m);
d=1;
md=0;
//--- cycle
while(MathAbs(c/a)>CMath::m_machineepsilon)
{
//--- calculation
temp=b/a;
phi=phi+MathArctan(t*temp)+md*M_PI;
md=(int)((phi+pio2)/M_PI);
t=t*(1.0+temp)/(1.0-temp*t*t);
c=0.5*(a-b);
temp=MathSqrt(a*b);
a=0.5*(a+b);
b=temp;
d=d+d;
}
//--- change value
temp=(MathArctan(t)+md*M_PI)/(d*a);
//--- check
if(s<0)
temp=-temp;
//--- return result
return(temp+npio2*k);
}
//+------------------------------------------------------------------+
//| Complete elliptic integral of the second kind |
//| Approximates the integral |
//| pi/2 |
//| - |
//| | | 2 |
//| E(m) = | sqrt( 1 - m sin t ) dt |
//| | | |
//| - |
//| 0 |
//| using the approximation |
//| P(x) - x log x Q(x). |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 1 10000 2.1e-16 7.3e-17 |
//+------------------------------------------------------------------+
double CElliptic::EllipticIntegralE(double m)
{
//--- create variables
double p=0;
double q=0;
//--- check
if(!CAp::Assert(m>=0.0 && m<=1.0,__FUNCTION__+": m<0 or m>1"))
return(EMPTY_VALUE);
//--- change value
m=1.0-m;
//--- check
if(m==0.0)
return(1.0);
//--- calculation
p=1.53552577301013293365E-4;
p=p*m+2.50888492163602060990E-3;
p=p*m+8.68786816565889628429E-3;
p=p*m+1.07350949056076193403E-2;
p=p*m+7.77395492516787092951E-3;
p=p*m+7.58395289413514708519E-3;
p=p*m+1.15688436810574127319E-2;
p=p*m+2.18317996015557253103E-2;
p=p*m+5.68051945617860553470E-2;
p=p*m+4.43147180560990850618E-1;
p=p*m+1.00000000000000000299E0;
q=3.27954898576485872656E-5;
q=q*m+1.00962792679356715133E-3;
q=q*m+6.50609489976927491433E-3;
q=q*m+1.68862163993311317300E-2;
q=q*m+2.61769742454493659583E-2;
q=q*m+3.34833904888224918614E-2;
q=q*m+4.27180926518931511717E-2;
q=q*m+5.85936634471101055642E-2;
q=q*m+9.37499997197644278445E-2;
q=q*m+2.49999999999888314361E-1;
//--- return result
return(p-q*m*MathLog(m));
}
//+------------------------------------------------------------------+
//| Incomplete elliptic integral of the second kind |
//| Approximates the integral |
//| phi |
//| - |
//| | | |
//| | 2 |
//| E(phi_\m) = | sqrt( 1 - m sin t ) dt |
//| | |
//| | | |
//| - |
//| 0 |
//| of amplitude phi and modulus m, using the arithmetic - |
//| geometric mean algorithm. |
//| ACCURACY: |
//| Tested at random arguments with phi in [-10, 10] and m in |
//| [0, 1]. |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -10,10 150000 3.3e-15 1.4e-16 |
//+------------------------------------------------------------------+
double CElliptic::IncompleteEllipticIntegralE(const double phi,
const double m)
{
//--- create variables
double pio2=1.57079632679489661923;
double a=0;
double b=0;
double c=0;
double e=0;
double temp=0;
double lphi=0;
double t=0;
double ebig=0;
int d=0;
int md=0;
int npio2=0;
int s=0;
//--- check
if(m==0.0)
return(phi);
//--- change values
lphi=phi;
npio2=(int)MathFloor(lphi/pio2);
//--- check
if(npio2%2!=0)
npio2=npio2+1;
lphi=lphi-npio2*pio2;
//--- check
if(lphi<0.0)
{
lphi=-lphi;
s=-1;
}
else
s=1;
//--- change values
a=1.0-m;
ebig=EllipticIntegralE(m);
//--- check
if(a==0.0)
{
temp=MathSin(lphi);
//--- check
if(s<0)
temp=-temp;
//--- return result
return(temp+npio2*ebig);
}
//--- calculation
t=MathTan(lphi);
b=MathSqrt(a);
//--- Thanks to Brian Fitzgerald <fitzgb@mml0.meche.rpi.edu>
//--- for pointing out an instability near odd multiples of pi/2
if(MathAbs(t)>10)
{
//--- Transform the amplitude
e=1.0/(b*t);
//--- ... but avoid multiple recursions.
if(MathAbs(e)<10)
{
e=MathArctan(e);
temp=ebig+m*MathSin(lphi)*MathSin(e)-IncompleteEllipticIntegralE(e,m);
//--- check
if(s<0)
temp=-temp;
//--- return result
return(temp+npio2*ebig);
}
}
//--- change values
c=MathSqrt(m);
a=1.0;
d=1;
e=0.0;
md=0;
//--- cycle
while(MathAbs(c/a)>CMath::m_machineepsilon)
{
//--- calculation
temp=b/a;
lphi=lphi+MathArctan(t*temp)+md*M_PI;
md=(int)((lphi+pio2)/M_PI);
t=t*(1.0+temp)/(1.0-temp*t*t);
c=0.5*(a-b);
temp=MathSqrt(a*b);
a=0.5*(a+b);
b=temp;
d=d+d;
e=e+c*MathSin(lphi);
}
//--- change value
temp=ebig/EllipticIntegralK(m);
temp=temp*((MathArctan(t)+md*M_PI)/(d*a));
temp=temp+e;
//--- check
if(s<0)
temp=-temp;
//--- return result
return(temp+npio2*ebig);
}
//+------------------------------------------------------------------+
//| Exponential integral |
//+------------------------------------------------------------------+
class CExpIntegrals
{
public:
static double ExponentialIntegralEi(const double x);
static double ExponentialIntegralEn(const double x,const int n);
};
//+------------------------------------------------------------------+
//| Exponential integral Ei(x) |
//| x |
//| - t |
//| | | e |
//| Ei(x) = -|- --- dt . |
//| | | t |
//| - |
//| -inf |
//| Not defined for x <= 0. |
//| See also expn.c. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0,100 50000 8.6e-16 1.3e-16 |
//+------------------------------------------------------------------+
double CExpIntegrals::ExponentialIntegralEi(const double x)
{
//--- create variables
double eul=0.5772156649015328606065;
double f=0;
double f1=0;
double f2=0;
double w=0;
//--- check
if(x<=0.0)
return(0);
//--- check
if(x<2.0)
{
//--- calculation
f1=-5.350447357812542947283;
f1=f1*x+218.5049168816613393830;
f1=f1*x-4176.572384826693777058;
f1=f1*x+55411.76756393557601232;
f1=f1*x-331338.1331178144034309;
f1=f1*x+1592627.163384945414220;
f2=1.000000000000000000000;
f2=f2*x-52.50547959112862969197;
f2=f2*x+1259.616186786790571525;
f2=f2*x-17565.49581973534652631;
f2=f2*x+149306.2117002725991967;
f2=f2*x-729494.9239640527645655;
f2=f2*x+1592627.163384945429726;
f=f1/f2;
//--- return result
return(eul+MathLog(x)+x*f);
}
//--- check
if(x<4.0)
{
//--- calculation
w=1/x;
f1=1.981808503259689673238E-2;
f1=f1*w-1.271645625984917501326;
f1=f1*w-2.088160335681228318920;
f1=f1*w+2.755544509187936721172;
f1=f1*w-4.409507048701600257171E-1;
f1=f1*w+4.665623805935891391017E-2;
f1=f1*w-1.545042679673485262580E-3;
f1=f1*w+7.059980605299617478514E-5;
f2=1.000000000000000000000;
f2=f2*w+1.476498670914921440652;
f2=f2*w+5.629177174822436244827E-1;
f2=f2*w+1.699017897879307263248E-1;
f2=f2*w+2.291647179034212017463E-2;
f2=f2*w+4.450150439728752875043E-3;
f2=f2*w+1.727439612206521482874E-4;
f2=f2*w+3.953167195549672482304E-5;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//--- check
if(x<8.0)
{
//--- calculation
w=1/x;
f1=-1.373215375871208729803;
f1=f1*w-7.084559133740838761406E-1;
f1=f1*w+1.580806855547941010501;
f1=f1*w-2.601500427425622944234E-1;
f1=f1*w+2.994674694113713763365E-2;
f1=f1*w-1.038086040188744005513E-3;
f1=f1*w+4.371064420753005429514E-5;
f1=f1*w+2.141783679522602903795E-6;
f2=1.000000000000000000000;
f2=f2*w+8.585231423622028380768E-1;
f2=f2*w+4.483285822873995129957E-1;
f2=f2*w+7.687932158124475434091E-2;
f2=f2*w+2.449868241021887685904E-2;
f2=f2*w+8.832165941927796567926E-4;
f2=f2*w+4.590952299511353531215E-4;
f2=f2*w+-4.729848351866523044863E-6;
f2=f2*w+2.665195537390710170105E-6;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//--- check
if(x<16.0)
{
//--- calculation
w=1/x;
f1=-2.106934601691916512584;
f1=f1*w+1.732733869664688041885;
f1=f1*w-2.423619178935841904839E-1;
f1=f1*w+2.322724180937565842585E-2;
f1=f1*w+2.372880440493179832059E-4;
f1=f1*w-8.343219561192552752335E-5;
f1=f1*w+1.363408795605250394881E-5;
f1=f1*w-3.655412321999253963714E-7;
f1=f1*w+1.464941733975961318456E-8;
f1=f1*w+6.176407863710360207074E-10;
f2=1.000000000000000000000;
f2=f2*w-2.298062239901678075778E-1;
f2=f2*w+1.105077041474037862347E-1;
f2=f2*w-1.566542966630792353556E-2;
f2=f2*w+2.761106850817352773874E-3;
f2=f2*w-2.089148012284048449115E-4;
f2=f2*w+1.708528938807675304186E-5;
f2=f2*w-4.459311796356686423199E-7;
f2=f2*w+1.394634930353847498145E-8;
f2=f2*w+6.150865933977338354138E-10;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//--- check
if(x<32.0)
{
//--- calculation
w=1/x;
f1=-2.458119367674020323359E-1;
f1=f1*w-1.483382253322077687183E-1;
f1=f1*w+7.248291795735551591813E-2;
f1=f1*w-1.348315687380940523823E-2;
f1=f1*w+1.342775069788636972294E-3;
f1=f1*w-7.942465637159712264564E-5;
f1=f1*w+2.644179518984235952241E-6;
f1=f1*w-4.239473659313765177195E-8;
f2=1.000000000000000000000;
f2=f2*w-1.044225908443871106315E-1;
f2=f2*w-2.676453128101402655055E-1;
f2=f2*w+9.695000254621984627876E-2;
f2=f2*w-1.601745692712991078208E-2;
f2=f2*w+1.496414899205908021882E-3;
f2=f2*w-8.462452563778485013756E-5;
f2=f2*w+2.728938403476726394024E-6;
f2=f2*w-4.239462431819542051337E-8;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//--- check
if(x<64)
{
//--- calculation
w=1/x;
f1=1.212561118105456670844E-1;
f1=f1*w-5.823133179043894485122E-1;
f1=f1*w+2.348887314557016779211E-1;
f1=f1*w-3.040034318113248237280E-2;
f1=f1*w+1.510082146865190661777E-3;
f1=f1*w-2.523137095499571377122E-5;
f2=1.000000000000000000000;
f2=f2*w-1.002252150365854016662;
f2=f2*w+2.928709694872224144953E-1;
f2=f2*w-3.337004338674007801307E-2;
f2=f2*w+1.560544881127388842819E-3;
f2=f2*w-2.523137093603234562648E-5;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//--- calculation
w=1/x;
f1=-7.657847078286127362028E-1;
f1=f1*w+6.886192415566705051750E-1;
f1=f1*w-2.132598113545206124553E-1;
f1=f1*w+3.346107552384193813594E-2;
f1=f1*w-3.076541477344756050249E-3;
f1=f1*w+1.747119316454907477380E-4;
f1=f1*w-6.103711682274170530369E-6;
f1=f1*w+1.218032765428652199087E-7;
f1=f1*w-1.086076102793290233007E-9;
f2=1.000000000000000000000;
f2=f2*w-1.888802868662308731041;
f2=f2*w+1.066691687211408896850;
f2=f2*w-2.751915982306380647738E-1;
f2=f2*w+3.930852688233823569726E-2;
f2=f2*w-3.414684558602365085394E-3;
f2=f2*w+1.866844370703555398195E-4;
f2=f2*w-6.345146083130515357861E-6;
f2=f2*w+1.239754287483206878024E-7;
f2=f2*w-1.086076102793126632978E-9;
f=f1/f2;
//--- return result
return(MathExp(x)*w*(1+w*f));
}
//+------------------------------------------------------------------+
//| Exponential integral En(x) |
//| Evaluates the exponential integral |
//| inf. |
//| - |
//| | | -xt |
//| | e |
//| E (x) = | ---- dt. |
//| n | n |
//| | | t |
//| - |
//| 1 |
//| Both n and x must be nonnegative. |
//| The routine employs either a power series, a continued |
//| fraction, or an asymptotic formula depending on the |
//| relative values of n and x. |
//| ACCURACY: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE 0, 30 10000 1.7e-15 3.6e-16 |
//+------------------------------------------------------------------+
double CExpIntegrals::ExponentialIntegralEn(const double x,const int n)
{
//--- create variables
double eul=0.57721566490153286060;
double big=1.44115188075855872*MathPow(10,17);
double result=0;
double r=0;
double t=0;
double yk=0;
double xk=0;
double pk=0;
double pkm1=0;
double pkm2=0;
double qk=0;
double qkm1=0;
double qkm2=0;
double psi=0;
double z=0;
int i=0;
int k=0;
//--- check
if(((n<0 || x<0.0) || x>170) || (x==0.0 && n<2))
return(-1);
//--- check
if(x==0.0)
return(1.0/(double)(n-1));
//--- check
if(n==0)
return(MathExp(-x)/x);
//--- check
if(n>5000)
{
//--- calculation
xk=x+n;
yk=1/(xk*xk);
t=n;
result=yk*t*(6*x*x-8*t*x+t*t);
result=yk*(result+t*(t-2.0*x));
result=yk*(result+t);
result=(result+1)*MathExp(-x)/xk;
//--- return result
return(result);
}
//--- check
if(x<=1.0)
{
//--- calculation
psi=-eul-MathLog(x);
for(i=1; i<=n-1; i++)
psi=psi+1.0/(double)i;
z=-x;
xk=0;
yk=1;
pk=1-n;
//--- check
if(n==1)
result=0.0;
else
result=1.0/pk;
do
{
//--- change values
xk=xk+1;
yk=yk*z/xk;
pk=pk+1;
//--- check
if(pk!=0.0)
result=result+yk/pk;
//--- check
if(result!=0.0)
t=MathAbs(yk/result);
else
t=1;
}
while(t>=CMath::m_machineepsilon);
t=1;
for(i=1; i<=n-1; i++)
t=t*z/i;
//--- return result
return(psi*t-result);
}
else
{
//--- calculation
k=1;
pkm2=1;
qkm2=x;
pkm1=1.0;
qkm1=x+n;
result=pkm1/qkm1;
do
{
k=k+1;
//--- check
if(k%2==1)
{
yk=1;
xk=n+(double)(k-1)/2.0;
}
else
{
yk=x;
xk=(double)k/2.0;
}
//--- change values
pk=pkm1*yk+pkm2*xk;
qk=qkm1*yk+qkm2*xk;
//--- check
if(qk!=0.0)
{
r=pk/qk;
t=MathAbs((result-r)/r);
result=r;
}
else
t=1;
//--- change values
pkm2=pkm1;
pkm1=pk;
qkm2=qkm1;
qkm1=qk;
//--- check
if(MathAbs(pk)>big)
{
//--- change values
pkm2=pkm2/big;
pkm1=pkm1/big;
qkm2=qkm2/big;
qkm1=qkm1/big;
}
}
while(t>=CMath::m_machineepsilon);
result=result*MathExp(-x);
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| F distribution |
//+------------------------------------------------------------------+
class CFDistr
{
public:
static double FDistribution(const int a,const int b,const double x);
static double FComplDistribution(const int a,const int b,const double x);
static double InvFDistribution(const int a,const int b,const double y);
};
//+------------------------------------------------------------------+
//| F distribution |
//| Returns the area from zero to x under the F density |
//| function (also known as Snedcor's density or the |
//| variance ratio density). This is the density |
//| of x = (u1/df1)/(u2/df2), where u1 and u2 are random |
//| variables having Chi square distributions with df1 |
//| and df2 degrees of freedom, respectively. |
//| The incomplete beta integral is used, according to the |
//| formula |
//| P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ). |
//| The arguments a and b are greater than zero, and x is |
//| nonnegative. |
//| ACCURACY: |
//| Tested at random points (a,b,x). |
//| x a,b Relative error: |
//| arithmetic domain domain # trials peak rms |
//| IEEE 0,1 0,100 100000 9.8e-15 1.7e-15 |
//| IEEE 1,5 0,100 100000 6.5e-15 3.5e-16 |
//| IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12 |
//| IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13 |
//+------------------------------------------------------------------+
double CFDistr::FDistribution(const int a,const int b,const double x)
{
//--- check
if(!CAp::Assert(a>=1 && b>=1 && x>=0,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- calculation
double w=a*x;
w/=b+w;
//--- return result
return(CIncBetaF::IncompleteBeta(0.5*a,0.5*b,w));
}
//+------------------------------------------------------------------+
//| Complemented F distribution |
//| Returns the area from x to infinity under the F density |
//| function (also known as Snedcor's density or the |
//| variance ratio density). |
//| inf. |
//| - |
//| 1 | | a-1 b-1 |
//| 1-P(x) = ------ | t (1-t) dt |
//| B(a,b) | | |
//| - |
//| x |
//| The incomplete beta integral is used, according to the |
//| formula |
//| P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ). |
//| ACCURACY: |
//| Tested at random points (a,b,x) in the indicated intervals. |
//| x a,b Relative error: |
//| arithmetic domain domain # trials peak rms |
//| IEEE 0,1 1,100 100000 3.7e-14 5.9e-16 |
//| IEEE 1,5 1,100 100000 8.0e-15 1.6e-15 |
//| IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13 |
//| IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12 |
//+------------------------------------------------------------------+
double CFDistr::FComplDistribution(const int a,const int b,const double x)
{
//--- check
if(!CAp::Assert((a>=1 && b>=1) && x>=0.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- calculation
double w=b/(b+a*x);
//--- return result
return(CIncBetaF::IncompleteBeta(0.5*b,0.5*a,w));
}
//+------------------------------------------------------------------+
//| Inverse of complemented F distribution |
//| Finds the F density argument x such that the integral |
//| from x to infinity of the F density is equal to the |
//| given probability p. |
//| This is accomplished using the inverse beta integral |
//| function and the relations |
//| z = incbi( df2/2, df1/2, p ) |
//| x = df2 (1-z) / (df1 z). |
//| Note: the following relations hold for the inverse of |
//| the uncomplemented F distribution: |
//| z = incbi( df1/2, df2/2, p ) |
//| x = df2 z / (df1 (1-z)). |
//| ACCURACY: |
//| Tested at random points (a,b,p). |
//| a,b Relative error: |
//| arithmetic domain # trials peak rms |
//| For p between .001 and 1: |
//| IEEE 1,100 100000 8.3e-15 4.7e-16 |
//| IEEE 1,10000 100000 2.1e-11 1.4e-13 |
//| For p between 10^-6 and 10^-3: |
//| IEEE 1,100 50000 1.3e-12 8.4e-15 |
//| IEEE 1,10000 50000 3.0e-12 4.8e-14 |
//+------------------------------------------------------------------+
double CFDistr::InvFDistribution(const int a,const int b,const double y)
{
//--- create variables
double result=0;
double w=0;
//--- check
if(!CAp::Assert(((a>=1 && b>=1) && y>0.0) && y<=1.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- Compute probability for x=0.5
w=CIncBetaF::IncompleteBeta(0.5*b,0.5*a,0.5);
//--- If that is greater than y,then the solution w < .5
//--- Otherwise,solve at 1-y to remove cancellation in (b - b*w)
if(w>y || y<0.001)
{
//--- calculation
w=CIncBetaF::InvIncompleteBeta(0.5*b,0.5*a,y);
result=(b-b*w)/(a*w);
}
else
{
//--- calculation
w=CIncBetaF::InvIncompleteBeta(0.5*a,0.5*b,1.0-y);
result=b*w/(a*(1.0-w));
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Fresnel integral |
//+------------------------------------------------------------------+
class CFresnel
{
public:
static void FresnelIntegral(double x,double &c,double &s);
};
//+------------------------------------------------------------------+
//| Fresnel integral |
//| Evaluates the Fresnel integrals |
//| x |
//| - |
//| | | |
//| C(x) = | cos(pi/2 t**2) dt, |
//| | | |
//| - |
//| 0 |
//| x |
//| - |
//| | | |
//| S(x) = | sin(pi/2 t**2) dt. |
//| | | |
//| - |
//| 0 |
//| The integrals are evaluated by a power series for x < 1. |
//| For x >= 1 auxiliary functions f(x) and g(x) are employed |
//| such that |
//| C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 ) |
//| S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 ) |
//| ACCURACY: |
//| Relative error. |
//| Arithmetic function domain # trials peak rms |
//| IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16|
//| IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16|
//+------------------------------------------------------------------+
void CFresnel::FresnelIntegral(double x,double &c,double &s)
{
//--- create variables
double xxa=0;
double f=0;
double g=0;
double cc=0;
double ss=0;
double t=0;
double u=0;
double x2=0;
double sn=0;
double sd=0;
double cn=0;
double cd=0;
double fn=0;
double fd=0;
double gn=0;
double gd=0;
double mpi=0;
double mpio2=0;
//--- initialization
mpi=3.14159265358979323846;
mpio2=1.57079632679489661923;
xxa=x;
x=MathAbs(xxa);
x2=x*x;
//--- check
if(x2<2.5625)
{
//--- calculation
t=x2*x2;
sn=-2.99181919401019853726E3;
sn=sn*t+7.08840045257738576863E5;
sn=sn*t-6.29741486205862506537E7;
sn=sn*t+2.54890880573376359104E9;
sn=sn*t-4.42979518059697779103E10;
sn=sn*t+3.18016297876567817986E11;
sd=1.00000000000000000000E0;
sd=sd*t+2.81376268889994315696E2;
sd=sd*t+4.55847810806532581675E4;
sd=sd*t+5.17343888770096400730E6;
sd=sd*t+4.19320245898111231129E8;
sd=sd*t+2.24411795645340920940E10;
sd=sd*t+6.07366389490084639049E11;
cn=-4.98843114573573548651E-8;
cn=cn*t+9.50428062829859605134E-6;
cn=cn*t-6.45191435683965050962E-4;
cn=cn*t+1.88843319396703850064E-2;
cn=cn*t-2.05525900955013891793E-1;
cn=cn*t+9.99999999999999998822E-1;
cd=3.99982968972495980367E-12;
cd=cd*t+9.15439215774657478799E-10;
cd=cd*t+1.25001862479598821474E-7;
cd=cd*t+1.22262789024179030997E-5;
cd=cd*t+8.68029542941784300606E-4;
cd=cd*t+4.12142090722199792936E-2;
cd=cd*t+1.00000000000000000118E0;
s=MathSign(xxa)*x*x2*sn/sd;
c=MathSign(xxa)*x*cn/cd;
//--- exit the function
return;
}
//--- check
if(x>36974.0)
{
//--- calculation
c=MathSign(xxa)*0.5;
s=MathSign(xxa)*0.5;
//--- exit the function
return;
}
//--- calculation
x2=x*x;
t=mpi*x2;
u=1/(t*t);
t=1/t;
fn=4.21543555043677546506E-1;
fn=fn*u+1.43407919780758885261E-1;
fn=fn*u+1.15220955073585758835E-2;
fn=fn*u+3.45017939782574027900E-4;
fn=fn*u+4.63613749287867322088E-6;
fn=fn*u+3.05568983790257605827E-8;
fn=fn*u+1.02304514164907233465E-10;
fn=fn*u+1.72010743268161828879E-13;
fn=fn*u+1.34283276233062758925E-16;
fn=fn*u+3.76329711269987889006E-20;
fd=1.00000000000000000000E0;
fd=fd*u+7.51586398353378947175E-1;
fd=fd*u+1.16888925859191382142E-1;
fd=fd*u+6.44051526508858611005E-3;
fd=fd*u+1.55934409164153020873E-4;
fd=fd*u+1.84627567348930545870E-6;
fd=fd*u+1.12699224763999035261E-8;
fd=fd*u+3.60140029589371370404E-11;
fd=fd*u+5.88754533621578410010E-14;
fd=fd*u+4.52001434074129701496E-17;
fd=fd*u+1.25443237090011264384E-20;
gn=5.04442073643383265887E-1;
gn=gn*u+1.97102833525523411709E-1;
gn=gn*u+1.87648584092575249293E-2;
gn=gn*u+6.84079380915393090172E-4;
gn=gn*u+1.15138826111884280931E-5;
gn=gn*u+9.82852443688422223854E-8;
gn=gn*u+4.45344415861750144738E-10;
gn=gn*u+1.08268041139020870318E-12;
gn=gn*u+1.37555460633261799868E-15;
gn=gn*u+8.36354435630677421531E-19;
gn=gn*u+1.86958710162783235106E-22;
gd=1.00000000000000000000E0;
gd=gd*u+1.47495759925128324529E0;
gd=gd*u+3.37748989120019970451E-1;
gd=gd*u+2.53603741420338795122E-2;
gd=gd*u+8.14679107184306179049E-4;
gd=gd*u+1.27545075667729118702E-5;
gd=gd*u+1.04314589657571990585E-7;
gd=gd*u+4.60680728146520428211E-10;
gd=gd*u+1.10273215066240270757E-12;
gd=gd*u+1.38796531259578871258E-15;
gd=gd*u+8.39158816283118707363E-19;
gd=gd*u+1.86958710162783236342E-22;
//--- get result
f=1-u*fn/fd;
g=t*gn/gd;
t=mpio2*x2;
cc=MathCos(t);
ss=MathSin(t);
t=mpi*x;
c=0.5+(f*ss-g*cc)/t;
s=0.5-(f*cc+g*ss)/t;
c=c*MathSign(xxa);
s=s*MathSign(xxa);
}
//+------------------------------------------------------------------+
//| Hermite polynomial |
//+------------------------------------------------------------------+
class CHermite
{
public:
static double HermiteCalculate(const int n,const double x);
static double HermiteSum(double &c[],const int n,const double x);
static void HermiteCoefficients(const int n,double &c[]);
};
//+------------------------------------------------------------------+
//| Calculation of the value of the Hermite polynomial. |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Hermite polynomial Hn at x |
//+------------------------------------------------------------------+
double CHermite::HermiteCalculate(const int n,const double x)
{
//--- create variables
double result=0;
double a=1.0;
double b=2.0*x;
//--- Special cases: N=0 or N=1
if(n==0)
return(a);
//--- check
if(n==1)
return(b);
//--- General case: N>=2
for(int i=2; i<=n; i++)
{
result=2*x*b-2*(i-1)*a;
a=b;
b=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Summation of Hermite polynomials using Clenshaw?s recurrence |
//| formula. |
//| This routine calculates |
//| c[0]*H0(x) + c[1]*H1(x) + ... + c[N]*HN(x) |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Hermite polynomial at x |
//+------------------------------------------------------------------+
double CHermite::HermiteSum(double &c[],const int n,const double x)
{
//--- create variables
double result=0;
double b1=0;
double b2=0;
//--- calculation
for(int i=n; i>=0; i--)
{
result=2*(x*b1-(i+1)*b2)+c[i];
b2=b1;
b1=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Representation of Hn as C[0] + C[1]*X + ... + C[N]*X^N |
//| Input parameters: |
//| N - polynomial degree, n>=0 |
//| Output parameters: |
//| C - coefficients |
//+------------------------------------------------------------------+
void CHermite::HermiteCoefficients(const int n,double &c[])
{
//--- allocation
ArrayResize(c,n+1);
//--- initialization
ArrayInitialize(c,0.0);
//--- calculation
c[n]=MathExp(n*MathLog(2));
for(int i=0; i<=n/2-1; i++)
c[n-2*(i+1)]=-(c[n-2*i]*(n-2*i)*(n-2*i-1)/4/(i+1));
}
//+------------------------------------------------------------------+
//| Jacobian elliptic functions |
//+------------------------------------------------------------------+
class CJacobianElliptic
{
public:
static void JacobianEllipticFunctions(const double u,const double m,double &sn,double &cn,double &dn,double &ph);
};
//+------------------------------------------------------------------+
//| Jacobian Elliptic Functions |
//| Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), |
//| and dn(u|m) of parameter m between 0 and 1, and real |
//| argument u. |
//| These functions are periodic, with quarter-period on the |
//| real axis equal to the complete elliptic integral |
//| ellpk(1.0-m). |
//| Relation to incomplete elliptic integral: |
//| If u = ellik(phi,m), then sn(u|m) = sin(phi), |
//| and cn(u|m) = cos(phi). Phi is called the amplitude of u. |
//| Computation is by means of the arithmetic-geometric mean |
//| algorithm, except when m is within 1e-9 of 0 or 1. In the |
//| latter case with m close to 1, the approximation applies |
//| only for phi < pi/2. |
//| ACCURACY: |
//| Tested at random points with u between 0 and 10, m between |
//| 0 and 1. |
//| Absolute error (* = relative error): |
//| arithmetic function # trials peak rms |
//| IEEE phi 10000 9.2e-16* 1.4e-16* |
//| IEEE sn 50000 4.1e-15 4.6e-16 |
//| IEEE cn 40000 3.6e-15 4.4e-16 |
//| IEEE dn 10000 1.3e-12 1.8e-14 |
//| Peak error observed in consistency check using addition |
//| theorem for sn(u+v) was 4e-16 (absolute). Also tested by |
//| the above relation to the incomplete elliptic integral. |
//| Accuracy deteriorates when u is large. |
//+------------------------------------------------------------------+
void CJacobianElliptic::JacobianEllipticFunctions(const double u,
const double m,
double &sn,double &cn,
double &dn,double &ph)
{
//--- create variables
double ai=0;
double b=0;
double phi=0;
double t=0;
double twon=0;
int i=0;
//--- create arrays
double a[];
double c[];
//--- initialization
sn=0;
cn=0;
dn=0;
ph=0;
//--- check
if(!CAp::Assert(m>=0.0 && m<=1.0,__FUNCTION__+": m<0 or m>1"))
return;
//--- allocation
ArrayResize(a,9);
ArrayResize(c,9);
//--- check
if(m<1.0e-9)
{
//--- calculation
t=MathSin(u);
b=MathCos(u);
ai=0.25*m*(u-t*b);
sn=t-ai*b;
cn=b+ai*t;
ph=u-ai;
dn=1.0-0.5*m*t*t;
//--- exit the function
return;
}
//--- check
if(m>=0.9999999999)
{
//--- calculation
ai=0.25*(1.0-m);
b=MathCosh(u);
t=MathTanh(u);
phi=1.0/b;
twon=b*MathSinh(u);
sn=t+ai*(twon-u)/(b*b);
ph=2.0*MathArctan(MathExp(u))-1.57079632679489661923+ai*(twon-u)/b;
ai=ai*t*phi;
cn=phi-ai*(twon-u);
dn=phi+ai*(twon+u);
//--- exit the function
return;
}
//--- change values
a[0]=1.0;
b=MathSqrt(1.0-m);
c[0]=MathSqrt(m);
twon=1.0;
i=0;
//--- cycle
while(MathAbs(c[i]/a[i])>CMath::m_machineepsilon)
{
//--- check
if(i>7)
{
//--- check
if(!CAp::Assert(false,__FUNCTION__+": overflow"))
return;
//--- break the cycle
break;
}
//--- calculation
ai=a[i];
i=i+1;
c[i]=0.5*(ai-b);
t=MathSqrt(ai*b);
a[i]=0.5*(ai+b);
b=t;
twon=twon*2.0;
}
phi=twon*a[i]*u;
//--- cycle
do
{
//--- calculation
t=c[i]*MathSin(phi)/a[i];
b=phi;
phi=(MathArcsin(t)+phi)/2.0;
i=i-1;
}
while(i!=0);
//--- get result
sn=MathSin(phi);
t=MathCos(phi);
cn=t;
dn=t/MathCos(phi-b);
ph=phi;
}
//+------------------------------------------------------------------+
//| Laguerre polynomial |
//+------------------------------------------------------------------+
class CLaguerre
{
public:
static double LaguerreCalculate(const int n,const double x);
static double LaguerreSum(double &c[],const int n,const double x);
static void LaguerreCoefficients(const int n,double &c[]);
};
//+------------------------------------------------------------------+
//| Calculation of the value of the Laguerre polynomial. |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Laguerre polynomial Ln at x |
//+------------------------------------------------------------------+
double CLaguerre::LaguerreCalculate(const int n,const double x)
{
//--- create variables
double result=0;
double a=0;
double b=0;
double i=0;
//--- initialization
result=1;
a=1;
b=1-x;
//--- check
if(n==1)
result=b;
//--- initialization
i=2;
//--- cycle
while(i<=n)
{
//--- calculation
result=((2*i-1-x)*b-(i-1)*a)/i;
a=b;
b=result;
i=i+1;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Summation of Laguerre polynomials using Clenshaw?s recurrence |
//| formula. |
//| This routine calculates c[0]*L0(x) + c[1]*L1(x) + ... + |
//| + c[N]*LN(x) |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Laguerre polynomial at x |
//+------------------------------------------------------------------+
double CLaguerre::LaguerreSum(double &c[],const int n,const double x)
{
//--- create variables
double result=0;
double b1=0;
double b2=0;
//--- calculation
for(int i=n; i>=0; i--)
{
result=(2*i+1-x)*b1/(i+1)-(i+1)*b2/(i+2)+c[i];
b2=b1;
b1=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Representation of Ln as C[0] + C[1]*X + ... + C[N]*X^N |
//| Input parameters: |
//| N - polynomial degree, n>=0 |
//| Output parameters: |
//| C - coefficients |
//+------------------------------------------------------------------+
void CLaguerre::LaguerreCoefficients(const int n,double &c[])
{
//--- allocation
ArrayResize(c,n+1);
//--- initialization
c[0]=1;
//--- calculation
for(int i=0; i<n; i++)
c[i+1]=-(c[i]*(n-i)/(i+1)/(i+1));
}
//+------------------------------------------------------------------+
//| Legendre polynomial |
//+------------------------------------------------------------------+
class CLegendre
{
public:
static double LegendreCalculate(const int n,const double x);
static double LegendreSum(double &c[],const int n,const double x);
static void LegendreCoefficients(const int n,double &c[]);
};
//+------------------------------------------------------------------+
//| Calculation of the value of the Legendre polynomial Pn. |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Legendre polynomial Pn at x |
//+------------------------------------------------------------------+
double CLegendre::LegendreCalculate(const int n,const double x)
{
//--- create variables
double result=1.0;
double a=1.0;
double b=x;
//--- check
if(n==0)
return(a);
//--- check
if(n==1)
return(b);
//--- calculation
for(int i=2; i<=n; i++)
{
result=((2*i-1)*x*b-(i-1)*a)/i;
a=b;
b=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Summation of Legendre polynomials using Clenshaw?s recurrence |
//| formula. |
//| This routine calculates |
//| c[0]*P0(x) + c[1]*P1(x) + ... + c[N]*PN(x) |
//| Parameters: |
//| n - degree, n>=0 |
//| x - argument |
//| Result: |
//| the value of the Legendre polynomial at x |
//+------------------------------------------------------------------+
double CLegendre::LegendreSum(double &c[],const int n,const double x)
{
//--- create variables
double result=0;
double b1=0;
double b2=0;
//--- calculation
for(int i=n; i>=0; i--)
{
result=(2*i+1)*x*b1/(i+1)-(i+1)*b2/(i+2)+c[i];
b2=b1;
b1=result;
}
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| Representation of Pn as C[0] + C[1]*X + ... + C[N]*X^N |
//| Input parameters: |
//| N - polynomial degree, n>=0 |
//| Output parameters: |
//| C - coefficients |
//+------------------------------------------------------------------+
void CLegendre::LegendreCoefficients(const int n,double &c[])
{
//--- allocation
ArrayResize(c,n+1);
//--- initialization
ArrayInitialize(c,0.0);
c[n]=1;
//--- calculation
for(int i=1; i<=n; i++)
c[n]=c[n]*(n+i)/2/i;
for(int i=0; i<=n/2-1; i++)
c[n-2*(i+1)]=-(c[n-2*i]*(n-2*i)*(n-2*i-1)/2/(i+1)/(2*(n-i)-1));
}
//+------------------------------------------------------------------+
//| Poisson distribution |
//+------------------------------------------------------------------+
class CPoissonDistr
{
public:
static double PoissonDistribution(const int k,const double m);
static double PoissonComplDistribution(const int k,const double m);
static double InvPoissonDistribution(const int k,const double y);
};
//+------------------------------------------------------------------+
//| Poisson distribution |
//| Returns the sum of the first k+1 terms of the Poisson |
//| distribution: |
//| k j |
//| -- -m m |
//| > e -- |
//| -- j! |
//| j=0 |
//| The terms are not summed directly; instead the incomplete |
//| gamma integral is employed, according to the relation |
//| y = pdtr( k, m ) = igamc( k+1, m ). |
//| The arguments must both be positive. |
//| ACCURACY: |
//| See incomplete gamma function |
//+------------------------------------------------------------------+
double CPoissonDistr::PoissonDistribution(const int k,const double m)
{
//--- check
if(!CAp::Assert(k>=0 && m>0.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- return result
return(CIncGammaF::IncompleteGammaC(k+1,m));
}
//+------------------------------------------------------------------+
//| Complemented Poisson distribution |
//| Returns the sum of the terms k+1 to infinity of the Poisson |
//| distribution: |
//| inf. j |
//| -- -m m |
//| > e -- |
//| -- j! |
//| j=k+1 |
//| The terms are not summed directly; instead the incomplete |
//| gamma integral is employed, according to the formula |
//| y = pdtrc( k, m ) = igam( k+1, m ). |
//| The arguments must both be positive. |
//| ACCURACY: |
//| See incomplete gamma function |
//+------------------------------------------------------------------+
double CPoissonDistr::PoissonComplDistribution(const int k,const double m)
{
//--- check
if(!CAp::Assert(k>=0 && m>0.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- return result
return(CIncGammaF::IncompleteGamma(k+1,m));
}
//+------------------------------------------------------------------+
//| Inverse Poisson distribution |
//| Finds the Poisson variable x such that the integral |
//| from 0 to x of the Poisson density is equal to the |
//| given probability y. |
//| This is accomplished using the inverse gamma integral |
//| function and the relation |
//| m = igami( k+1, y ). |
//| ACCURACY: |
//| See inverse incomplete gamma function |
//+------------------------------------------------------------------+
double CPoissonDistr::InvPoissonDistribution(const int k,const double y)
{
//--- check
if(!CAp::Assert((k>=0 && y>=0.0) && y<1.0,__FUNCTION__+": domain error"))
return(EMPTY_VALUE);
//--- return result
return(CIncGammaF::InvIncompleteGammaC(k+1,y));
}
//+------------------------------------------------------------------+
//| Psi function |
//+------------------------------------------------------------------+
class CPsiF
{
public:
static double Psi(double x);
};
//+------------------------------------------------------------------+
//| Psi (digamma) function |
//| d - |
//| psi(x) = -- ln | (x) |
//| dx |
//| is the logarithmic derivative of the gamma function. |
//| For integer x, |
//| n-1 |
//| - |
//| psi(n) = -EUL + > 1/k. |
//| - |
//| k=1 |
//| This formula is used for 0 < n <= 10. If x is negative, it |
//| is transformed to a positive argument by the reflection |
//| formula psi(1-x) = psi(x) + pi cot(pi x). |
//| For general positive x, the argument is made greater than 10 |
//| using the recurrence psi(x+1) = psi(x) + 1/x. |
//| Then the following asymptotic expansion is applied: |
//| inf. B |
//| - 2k |
//| psi(x) = log(x) - 1/2x - > ------- |
//| - 2k |
//| k=1 2k x |
//| where the B2k are Bernoulli numbers. |
//| ACCURACY: |
//| Relative error (except absolute when |psi| < 1): |
//| arithmetic domain # trials peak rms |
//| IEEE 0,30 30000 1.3e-15 1.4e-16 |
//| IEEE -30,0 40000 1.5e-15 2.2e-16 |
//+------------------------------------------------------------------+
double CPsiF::Psi(double x)
{
//--- create variables
double p=0;
double q=0;
double nz=0;
double s=0;
double w=0;
double y=0;
double z=0;
double polv=0;
int i=0;
int n=0;
int negative=0;
//--- check
if(x<=0.0)
{
negative=1;
q=x;
p=(int)MathFloor(q);
//--- check
if(p==q)
{
//--- check
if(!CAp::Assert(false,__FUNCTION__+": singularity in Psi(x)"))
return(EMPTY_VALUE);
//--- return result
return(CMath::m_maxrealnumber);
}
nz=q-p;
//--- check
if(nz!=0.5)
{
//--- check
if(nz>0.5)
{
p=p+1.0;
nz=q-p;
}
//--- change value
nz=M_PI/MathTan(M_PI*nz);
}
else
nz=0.0;
x=1.0-x;
}
//--- check
if(x<=10.0 && x==(double)((int)MathFloor(x)))
{
y=0.0;
n=(int)MathFloor(x);
//--- calculation
for(i=1; i<=n-1; i++)
{
w=i;
y=y+1.0/w;
}
y=y-0.57721566490153286061;
}
else
{
//--- change values
s=x;
w=0.0;
//--- cycle
while(s<10.0)
{
w=w+1.0/s;
s=s+1.0;
}
//--- check
if(s<1.0E17)
{
//--- calculation
z=1.0/(s*s);
polv=8.33333333333333333333E-2;
polv=polv*z-2.10927960927960927961E-2;
polv=polv*z+7.57575757575757575758E-3;
polv=polv*z-4.16666666666666666667E-3;
polv=polv*z+3.96825396825396825397E-3;
polv=polv*z-8.33333333333333333333E-3;
polv=polv*z+8.33333333333333333333E-2;
y=z*polv;
}
else
y=0.0;
//--- change value
y=MathLog(s)-0.5/s-y-w;
}
//--- check
if(negative!=0)
y=y-nz;
//--- return result
return(y);
}
//+------------------------------------------------------------------+
//| Student's t distribution |
//+------------------------------------------------------------------+
class CStudenttDistr
{
public:
static double StudenttDistribution(const int k,const double t);
static double InvStudenttDistribution(const int k,double p);
};
//+------------------------------------------------------------------+
//| Student's t distribution |
//| Computes the integral from minus infinity to t of the Student |
//| t distribution with integer k > 0 degrees of freedom: |
//| t |
//| - |
//| | | |
//| - | 2 -(k+1)/2 |
//| | ( (k+1)/2 ) | ( x ) |
//| ---------------------- | ( 1 + --- ) dx |
//| - | ( k ) |
//| sqrt( k pi ) | ( k/2 ) | |
//| | | |
//| - |
//| -inf. |
//| Relation to incomplete beta integral: |
//| 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z ) |
//| where |
//| z = k/(k + t**2). |
//| For t < -2, this is the method of computation. For higher t, |
//| a direct method is derived from integration by parts. |
//| Since the function is symmetric about t=0, the area under the |
//| right tail of the density is found by calling the function |
//| with -t instead of t. |
//| ACCURACY: |
//| Tested at random 1<=k<=25. The "domain" refers to t. |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE -100,-2 50000 5.9e-15 1.4e-15 |
//| IEEE -2,100 500000 2.7e-15 4.9e-17 |
//+------------------------------------------------------------------+
double CStudenttDistr::StudenttDistribution(const int k,const double t)
{
//--- check
if(!CAp::Assert(k>0,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- check
if(t==0)
return(0.5);
//--- create variables
double x=0;
double rk=0;
double z=0;
double f=0;
double tz=0;
double p=0;
double xsqk=0;
int j=0;
//--- check
if(t<-2.0)
{
rk=k;
z=rk/(rk+t*t);
//--- return result
return(0.5*CIncBetaF::IncompleteBeta(0.5*rk,0.5,z));
}
//--- check
if(t<0)
x=-t;
else
x=t;
//--- change values
rk=k;
z=1.0+x*x/rk;
//--- check
if(k%2!=0)
{
xsqk=x/MathSqrt(rk);
p=MathArctan(xsqk);
//--- check
if(k>1)
{
f=1.0;
tz=1.0;
j=3;
//--- cycle
while(j<k-1 && tz/f>CMath::m_machineepsilon)
{
//--- calculation
tz*=(j-1)/(z*j);
f+=tz;
j+=2;
}
p+=f*xsqk/z;
}
p*=2.0/M_PI;
}
else
{
f=1.0;
tz=1.0;
j=2;
//--- cycle
while(j<k-1 && tz/f>CMath::m_machineepsilon)
{
//--- calculation
tz*=(j-1)/(z*j);
f+=tz;
j+=2;
}
p=f*x/MathSqrt(z*rk);
}
//--- check
if(t<0)
p=-p;
//--- return result
return(0.5+0.5*p);
}
//+------------------------------------------------------------------+
//| Functional inverse of Student's t distribution |
//| Given probability p, finds the argument t such that stdtr(k,t) |
//| is equal to p. |
//| ACCURACY: |
//| Tested at random 1<=k<=100. The "domain" refers to p: |
//| Relative error: |
//| arithmetic domain # trials peak rms |
//| IEEE .001,.999 25000 5.7e-15 8.0e-16 |
//| IEEE 10^-6,.001 25000 2.0e-12 2.9e-14 |
//+------------------------------------------------------------------+
double CStudenttDistr::InvStudenttDistribution(const int k,double p)
{
//--- create variables
double t=0;
double rk=k;
double z=0;
int rflg=0;
CIncBetaF ibetaf;
//--- check
if(!CAp::Assert(k>0 && p>0 && p<1,__FUNCTION__+": the error variable"))
return(EMPTY_VALUE);
//--- check
if(p>0.25 && p<0.75)
{
//--- check
if(p==0.5)
return(0);
//--- change values
z=1.0-2.0*p;
z=CIncBetaF::InvIncompleteBeta(0.5,0.5*rk,MathAbs(z));
t=MathSqrt(rk*z/(1.0-z));
//--- check
if(p<0.5)
t=-t;
//--- return result
return(t);
}
//--- 0<p<=0.25 or 0.75<=p<1
rflg=-1;
//--- check
if(p>=0.5)
{
p=1.0-p;
rflg=1;
}
//--- change value
z=CIncBetaF::InvIncompleteBeta(0.5*rk,0.5,2.0*p);
//--- check
if(CMath::m_maxrealnumber*z<rk)
return(rflg*CMath::m_maxrealnumber);
t=MathSqrt(rk/z-rk);
//--- return result
return(rflg*t);
}
//+------------------------------------------------------------------+
//| Trigonometric integrals |
//+------------------------------------------------------------------+
class CTrigIntegrals
{
public:
static void SineCosineIntegrals(double x,double &si,double &ci);
static void HyperbolicSineCosineIntegrals(double x,double &shi,double &chi);
private:
static void ChebIterationShiChi(const double x,const double c,double &b0,double &b1,double &b2);
};
//+------------------------------------------------------------------+
//| Sine and cosine integrals |
//| Evaluates the integrals |
//| x |
//| - |
//| | cos t - 1 |
//| Ci(x) = eul + ln x + | --------- dt, |
//| | t |
//| - |
//| 0 |
//| x |
//| - |
//| | sin t |
//| Si(x) = | ----- dt |
//| | t |
//| - |
//| 0 |
//| where eul = 0.57721566490153286061 is Euler's constant. |
//| The integrals are approximated by rational functions. |
//| For x > 8 auxiliary functions f(x) and g(x) are employed |
//| such that |
//| Ci(x) = f(x) sin(x) - g(x) cos(x) |
//| Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x) |
//| ACCURACY: |
//| Test interval = [0,50]. |
//| Absolute error, except relative when > 1: |
//| arithmetic function # trials peak rms |
//| IEEE Si 30000 4.4e-16 7.3e-17 |
//| IEEE Ci 30000 6.9e-16 5.1e-17 |
//+------------------------------------------------------------------+
void CTrigIntegrals::SineCosineIntegrals(double x,double &si,double &ci)
{
//--- create variables
double z=0;
double c=0;
double s=0;
double f=0;
double g=0;
int sg=0;
double sn=0;
double sd=0;
double cn=0;
double cd=0;
double fn=0;
double fd=0;
double gn=0;
double gd=0;
//--- initialization
si=0;
ci=0;
//--- check
if(x<0.0)
{
sg=-1;
x=-x;
}
else
sg=0;
//--- check
if(x==0.0)
{
si=0;
ci=-CMath::m_maxrealnumber;
//--- exit the function
return;
}
//--- check
if(x>1.0E9)
{
si=1.570796326794896619-MathCos(x)/x;
ci=MathSin(x)/x;
//--- exit the function
return;
}
//--- check
if(x<=4.0)
{
//--- calculation
z=x*x;
sn=-8.39167827910303881427E-11;
sn=sn*z+4.62591714427012837309E-8;
sn=sn*z-9.75759303843632795789E-6;
sn=sn*z+9.76945438170435310816E-4;
sn=sn*z-4.13470316229406538752E-2;
sn=sn*z+1.00000000000000000302E0;
sd=2.03269266195951942049E-12;
sd=sd*z+1.27997891179943299903E-9;
sd=sd*z+4.41827842801218905784E-7;
sd=sd*z+9.96412122043875552487E-5;
sd=sd*z+1.42085239326149893930E-2;
sd=sd*z+9.99999999999999996984E-1;
s=x*sn/sd;
cn=2.02524002389102268789E-11;
cn=cn*z-1.35249504915790756375E-8;
cn=cn*z+3.59325051419993077021E-6;
cn=cn*z-4.74007206873407909465E-4;
cn=cn*z+2.89159652607555242092E-2;
cn=cn*z-1.00000000000000000080E0;
cd=4.07746040061880559506E-12;
cd=cd*z+3.06780997581887812692E-9;
cd=cd*z+1.23210355685883423679E-6;
cd=cd*z+3.17442024775032769882E-4;
cd=cd*z+5.10028056236446052392E-2;
cd=cd*z+4.00000000000000000080E0;
c=z*cn/cd;
//--- check
if(sg!=0)
s=-s;
//--- calculation
si=s;
ci=0.57721566490153286061+MathLog(x)+c;
//--- exit the function
return;
}
//--- change values
s=MathSin(x);
c=MathCos(x);
z=1.0/(x*x);
//--- check
if(x<8.0)
{
//--- calculation
fn=4.23612862892216586994E0;
fn=fn*z+5.45937717161812843388E0;
fn=fn*z+1.62083287701538329132E0;
fn=fn*z+1.67006611831323023771E-1;
fn=fn*z+6.81020132472518137426E-3;
fn=fn*z+1.08936580650328664411E-4;
fn=fn*z+5.48900223421373614008E-7;
fd=1.00000000000000000000E0;
fd=fd*z+8.16496634205391016773E0;
fd=fd*z+7.30828822505564552187E0;
fd=fd*z+1.86792257950184183883E0;
fd=fd*z+1.78792052963149907262E-1;
fd=fd*z+7.01710668322789753610E-3;
fd=fd*z+1.10034357153915731354E-4;
fd=fd*z+5.48900252756255700982E-7;
f=fn/(x*fd);
gn=8.71001698973114191777E-2;
gn=gn*z+6.11379109952219284151E-1;
gn=gn*z+3.97180296392337498885E-1;
gn=gn*z+7.48527737628469092119E-2;
gn=gn*z+5.38868681462177273157E-3;
gn=gn*z+1.61999794598934024525E-4;
gn=gn*z+1.97963874140963632189E-6;
gn=gn*z+7.82579040744090311069E-9;
gd=1.00000000000000000000E0;
gd=gd*z+1.64402202413355338886E0;
gd=gd*z+6.66296701268987968381E-1;
gd=gd*z+9.88771761277688796203E-2;
gd=gd*z+6.22396345441768420760E-3;
gd=gd*z+1.73221081474177119497E-4;
gd=gd*z+2.02659182086343991969E-6;
gd=gd*z+7.82579218933534490868E-9;
g=z*gn/gd;
}
else
{
//--- calculation
fn=4.55880873470465315206E-1;
fn=fn*z+7.13715274100146711374E-1;
fn=fn*z+1.60300158222319456320E-1;
fn=fn*z+1.16064229408124407915E-2;
fn=fn*z+3.49556442447859055605E-4;
fn=fn*z+4.86215430826454749482E-6;
fn=fn*z+3.20092790091004902806E-8;
fn=fn*z+9.41779576128512936592E-11;
fn=fn*z+9.70507110881952024631E-14;
fd=1.00000000000000000000E0;
fd=fd*z+9.17463611873684053703E-1;
fd=fd*z+1.78685545332074536321E-1;
fd=fd*z+1.22253594771971293032E-2;
fd=fd*z+3.58696481881851580297E-4;
fd=fd*z+4.92435064317881464393E-6;
fd=fd*z+3.21956939101046018377E-8;
fd=fd*z+9.43720590350276732376E-11;
fd=fd*z+9.70507110881952025725E-14;
f=fn/(x*fd);
gn=6.97359953443276214934E-1;
gn=gn*z+3.30410979305632063225E-1;
gn=gn*z+3.84878767649974295920E-2;
gn=gn*z+1.71718239052347903558E-3;
gn=gn*z+3.48941165502279436777E-5;
gn=gn*z+3.47131167084116673800E-7;
gn=gn*z+1.70404452782044526189E-9;
gn=gn*z+3.85945925430276600453E-12;
gn=gn*z+3.14040098946363334640E-15;
gd=1.00000000000000000000E0;
gd=gd*z+1.68548898811011640017E0;
gd=gd*z+4.87852258695304967486E-1;
gd=gd*z+4.67913194259625806320E-2;
gd=gd*z+1.90284426674399523638E-3;
gd=gd*z+3.68475504442561108162E-5;
gd=gd*z+3.57043223443740838771E-7;
gd=gd*z+1.72693748966316146736E-9;
gd=gd*z+3.87830166023954706752E-12;
gd=gd*z+3.14040098946363335242E-15;
g=z*gn/gd;
}
si=1.570796326794896619-f*c-g*s;
//--- check
if(sg!=0)
si=-si;
//--- get result
ci=f*s-g*c;
}
//+------------------------------------------------------------------+
//| Hyperbolic sine and cosine integrals |
//| Approximates the integrals |
//| x |
//| - |
//| | | cosh t - 1 |
//| Chi(x) = eul + ln x + | ----------- dt, |
//| | | t |
//| - |
//| 0 |
//| x |
//| - |
//| | | sinh t |
//| Shi(x) = | ------ dt |
//| | | t |
//| - |
//| 0 |
//| where eul = 0.57721566490153286061 is Euler's constant. |
//| The integrals are evaluated by power series for x < 8 |
//| and by Chebyshev expansions for x between 8 and 88. |
//| For large x, both functions approach exp(x)/2x. |
//| Arguments greater than 88 in magnitude return MAXNUM. |
//| ACCURACY: |
//| Test interval 0 to 88. |
//| Relative error: |
//| arithmetic function # trials peak rms |
//| IEEE Shi 30000 6.9e-16 1.6e-16 |
//| Absolute error, except relative when |Chi| > 1: |
//| IEEE Chi 30000 8.4e-16 1.4e-16 |
//+------------------------------------------------------------------+
void CTrigIntegrals::HyperbolicSineCosineIntegrals(double x,
double &shi,
double &chi)
{
//--- create variables
double k=0;
double z=0;
double c=0;
double s=0;
double a=0;
int sg=0;
double b0=0;
double b1=0;
double b2=0;
//--- initialization
shi=0;
chi=0;
//--- check
if(x<0.0)
{
sg=-1;
x=-x;
}
else
sg=0;
//--- check
if(x==0.0)
{
shi=0;
chi=-CMath::m_maxrealnumber;
//--- exit the function
return;
}
//--- check
if(x<8.0)
{
//--- initialization
z=x*x;
a=1.0;
s=1.0;
c=0.0;
k=2.0;
do
{
//--- calculation
a=a*z/k;
c=c+a/k;
k=k+1.0;
a=a/k;
s=s+a/k;
k=k+1.0;
}
while(MathAbs(a/s)>=CMath::m_machineepsilon);
s=s*x;
}
else
{
//--- check
if(x<18.0)
{
//--- calculation
a=(576.0/x-52.0)/10.0;
k=MathExp(x)/x;
b0=1.83889230173399459482E-17;
b1=0.0;
//--- function calls
ChebIterationShiChi(a,-9.55485532279655569575E-17,b0,b1,b2);
ChebIterationShiChi(a,2.04326105980879882648E-16,b0,b1,b2);
ChebIterationShiChi(a,1.09896949074905343022E-15,b0,b1,b2);
ChebIterationShiChi(a,-1.31313534344092599234E-14,b0,b1,b2);
ChebIterationShiChi(a,5.93976226264314278932E-14,b0,b1,b2);
ChebIterationShiChi(a,-3.47197010497749154755E-14,b0,b1,b2);
ChebIterationShiChi(a,-1.40059764613117131000E-12,b0,b1,b2);
ChebIterationShiChi(a,9.49044626224223543299E-12,b0,b1,b2);
ChebIterationShiChi(a,-1.61596181145435454033E-11,b0,b1,b2);
ChebIterationShiChi(a,-1.77899784436430310321E-10,b0,b1,b2);
ChebIterationShiChi(a,1.35455469767246947469E-9,b0,b1,b2);
ChebIterationShiChi(a,-1.03257121792819495123E-9,b0,b1,b2);
ChebIterationShiChi(a,-3.56699611114982536845E-8,b0,b1,b2);
ChebIterationShiChi(a,1.44818877384267342057E-7,b0,b1,b2);
ChebIterationShiChi(a,7.82018215184051295296E-7,b0,b1,b2);
ChebIterationShiChi(a,-5.39919118403805073710E-6,b0,b1,b2);
ChebIterationShiChi(a,-3.12458202168959833422E-5,b0,b1,b2);
ChebIterationShiChi(a,8.90136741950727517826E-5,b0,b1,b2);
ChebIterationShiChi(a,2.02558474743846862168E-3,b0,b1,b2);
ChebIterationShiChi(a,2.96064440855633256972E-2,b0,b1,b2);
ChebIterationShiChi(a,1.11847751047257036625E0,b0,b1,b2);
//--- calculation
s=k*0.5*(b0-b2);
b0=-8.12435385225864036372E-18;
b1=0.0;
//--- function calls
ChebIterationShiChi(a,2.17586413290339214377E-17,b0,b1,b2);
ChebIterationShiChi(a,5.22624394924072204667E-17,b0,b1,b2);
ChebIterationShiChi(a,-9.48812110591690559363E-16,b0,b1,b2);
ChebIterationShiChi(a,5.35546311647465209166E-15,b0,b1,b2);
ChebIterationShiChi(a,-1.21009970113732918701E-14,b0,b1,b2);
ChebIterationShiChi(a,-6.00865178553447437951E-14,b0,b1,b2);
ChebIterationShiChi(a,7.16339649156028587775E-13,b0,b1,b2);
ChebIterationShiChi(a,-2.93496072607599856104E-12,b0,b1,b2);
ChebIterationShiChi(a,-1.40359438136491256904E-12,b0,b1,b2);
ChebIterationShiChi(a,8.76302288609054966081E-11,b0,b1,b2);
ChebIterationShiChi(a,-4.40092476213282340617E-10,b0,b1,b2);
ChebIterationShiChi(a,-1.87992075640569295479E-10,b0,b1,b2);
ChebIterationShiChi(a,1.31458150989474594064E-8,b0,b1,b2);
ChebIterationShiChi(a,-4.75513930924765465590E-8,b0,b1,b2);
ChebIterationShiChi(a,-2.21775018801848880741E-7,b0,b1,b2);
ChebIterationShiChi(a,1.94635531373272490962E-6,b0,b1,b2);
ChebIterationShiChi(a,4.33505889257316408893E-6,b0,b1,b2);
ChebIterationShiChi(a,-6.13387001076494349496E-5,b0,b1,b2);
ChebIterationShiChi(a,-3.13085477492997465138E-4,b0,b1,b2);
ChebIterationShiChi(a,4.97164789823116062801E-4,b0,b1,b2);
ChebIterationShiChi(a,2.64347496031374526641E-2,b0,b1,b2);
ChebIterationShiChi(a,1.11446150876699213025E0,b0,b1,b2);
c=k*0.5*(b0-b2);
}
else
{
//--- check
if(x<=88.0)
{
//--- calculation
a=(6336.0/x-212.0)/70.0;
k=MathExp(x)/x;
b0=-1.05311574154850938805E-17;
b1=0.0;
//--- function calls
ChebIterationShiChi(a,2.62446095596355225821E-17,b0,b1,b2);
ChebIterationShiChi(a,8.82090135625368160657E-17,b0,b1,b2);
ChebIterationShiChi(a,-3.38459811878103047136E-16,b0,b1,b2);
ChebIterationShiChi(a,-8.30608026366935789136E-16,b0,b1,b2);
ChebIterationShiChi(a,3.93397875437050071776E-15,b0,b1,b2);
ChebIterationShiChi(a,1.01765565969729044505E-14,b0,b1,b2);
ChebIterationShiChi(a,-4.21128170307640802703E-14,b0,b1,b2);
ChebIterationShiChi(a,-1.60818204519802480035E-13,b0,b1,b2);
ChebIterationShiChi(a,3.34714954175994481761E-13,b0,b1,b2);
ChebIterationShiChi(a,2.72600352129153073807E-12,b0,b1,b2);
ChebIterationShiChi(a,1.66894954752839083608E-12,b0,b1,b2);
ChebIterationShiChi(a,-3.49278141024730899554E-11,b0,b1,b2);
ChebIterationShiChi(a,-1.58580661666482709598E-10,b0,b1,b2);
ChebIterationShiChi(a,-1.79289437183355633342E-10,b0,b1,b2);
ChebIterationShiChi(a,1.76281629144264523277E-9,b0,b1,b2);
ChebIterationShiChi(a,1.69050228879421288846E-8,b0,b1,b2);
ChebIterationShiChi(a,1.25391771228487041649E-7,b0,b1,b2);
ChebIterationShiChi(a,1.16229947068677338732E-6,b0,b1,b2);
ChebIterationShiChi(a,1.61038260117376323993E-5,b0,b1,b2);
ChebIterationShiChi(a,3.49810375601053973070E-4,b0,b1,b2);
ChebIterationShiChi(a,1.28478065259647610779E-2,b0,b1,b2);
ChebIterationShiChi(a,1.03665722588798326712E0,b0,b1,b2);
//--- calculation
s=k*0.5*(b0-b2);
b0=8.06913408255155572081E-18;
b1=0.0;
//--- function calls
ChebIterationShiChi(a,-2.08074168180148170312E-17,b0,b1,b2);
ChebIterationShiChi(a,-5.98111329658272336816E-17,b0,b1,b2);
ChebIterationShiChi(a,2.68533951085945765591E-16,b0,b1,b2);
ChebIterationShiChi(a,4.52313941698904694774E-16,b0,b1,b2);
ChebIterationShiChi(a,-3.10734917335299464535E-15,b0,b1,b2);
ChebIterationShiChi(a,-4.42823207332531972288E-15,b0,b1,b2);
ChebIterationShiChi(a,3.49639695410806959872E-14,b0,b1,b2);
ChebIterationShiChi(a,6.63406731718911586609E-14,b0,b1,b2);
ChebIterationShiChi(a,-3.71902448093119218395E-13,b0,b1,b2);
ChebIterationShiChi(a,-1.27135418132338309016E-12,b0,b1,b2);
ChebIterationShiChi(a,2.74851141935315395333E-12,b0,b1,b2);
ChebIterationShiChi(a,2.33781843985453438400E-11,b0,b1,b2);
ChebIterationShiChi(a,2.71436006377612442764E-11,b0,b1,b2);
ChebIterationShiChi(a,-2.56600180000355990529E-10,b0,b1,b2);
ChebIterationShiChi(a,-1.61021375163803438552E-9,b0,b1,b2);
ChebIterationShiChi(a,-4.72543064876271773512E-9,b0,b1,b2);
ChebIterationShiChi(a,-3.00095178028681682282E-9,b0,b1,b2);
ChebIterationShiChi(a,7.79387474390914922337E-8,b0,b1,b2);
ChebIterationShiChi(a,1.06942765566401507066E-6,b0,b1,b2);
ChebIterationShiChi(a,1.59503164802313196374E-5,b0,b1,b2);
ChebIterationShiChi(a,3.49592575153777996871E-4,b0,b1,b2);
ChebIterationShiChi(a,1.28475387530065247392E-2,b0,b1,b2);
ChebIterationShiChi(a,1.03665693917934275131E0,b0,b1,b2);
c=k*0.5*(b0-b2);
}
else
{
//--- check
if(sg!=0)
shi=-CMath::m_maxrealnumber;
else
shi=CMath::m_maxrealnumber;
chi=CMath::m_maxrealnumber;
//--- exit the function
return;
}
}
}
//--- check
if(sg!=0)
s=-s;
//--- get result
shi=s;
chi=0.57721566490153286061+MathLog(x)+c;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CTrigIntegrals::ChebIterationShiChi(const double x,const double c,
double &b0,double &b1,double &b2)
{
//--- change values
b2=b1;
b1=b0;
b0=x*b1-b2+c;
}
//+------------------------------------------------------------------+
Reference in a new issue Copy permalink