Princeec13/mql5
Watch
0
0
Fork
You've already forked mql5
 1
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions
mql5/Include/Math/Alglib/integration.mqh

2980 lines
114 KiB
MQL5
RawPermalinkNormal ViewHistory

new files added
2025-08-05 01:57:33 -04:00
//+------------------------------------------------------------------+
//| integration.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"
#include "linalg.mqh"
#include "specialfunctions.mqh"
//+------------------------------------------------------------------+
//| Gauss quadrature formula |
//+------------------------------------------------------------------+
class CGaussQ
{
public:
static void GQGenerateRec(double &alpha[],double &beta[],const double mu0,const int n,int &info,double &x[],double &w[]);
static void GQGenerateGaussLobattoRec(double &calpha[],double &cbeta[],const double mu0,const double a,const double b,int n,int &info,double &x[],double &w[]);
static void GQGenerateGaussRadauRec(double &calpha[],double &cbeta[],const double mu0,const double a,int n,int &info,double &x[],double &w[]);
static void GQGenerateGaussLegendre(const int n,int &info,double &x[],double &w[]);
static void GQGenerateGaussJacobi(const int n,const double alpha,const double beta,int &info,double &x[],double &w[]);
static void GQGenerateGaussLaguerre(const int n,const double alpha,int &info,double &x[],double &w[]);
static void GQGenerateGaussHermite(const int n,int &info,double &x[],double &w[]);
};
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss quadrature formula |
//| The algorithm generates the N-point Gauss quadrature formula |
//| with weight function given by coefficients alpha and beta of a |
//| recurrence relation which generates a system of orthogonal |
//| polynomials: |
//| P-1(x) = 0 |
//| P0(x) = 1 |
//| Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x) |
//| and zeroth moment Mu0 |
//| Mu0 = integral(W(x)dx,a,b) |
//| INPUT PARAMETERS: |
//| Alpha - array[0..N-1], alpha coefficients |
//| Beta - array[0..N-1], beta coefficients |
//| Zero-indexed element is not used and may be |
//| arbitrary. Beta[I]>0. |
//| Mu0 - zeroth moment of the weight function. |
//| N - number of nodes of the quadrature formula, N>=1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -2 Beta[i]<=0 |
//| * -1 incorrect N was passed |
//| * 1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateRec(double &alpha[],double &beta[],
const double mu0,const int n,
int &info,double &x[],double &w[])
{
//--- create a variable
int i=0;
//--- create arrays
double d[];
double e[];
//--- create matrix
CMatrixDouble z;
//--- initialization
info=0;
//--- check
if(n<1)
{
info=-1;
return;
}
info=1;
//--- Initialize
ArrayResize(d,n);
ArrayResize(e,n);
//--- calculation
for(i=1; i<=n-1; i++)
{
d[i-1]=alpha[i-1];
//--- check
if(beta[i]<=0.0)
{
info=-2;
return;
}
e[i-1]=MathSqrt(beta[i]);
}
//--- change value
d[n-1]=alpha[n-1];
//--- EVD
if(!CEigenVDetect::SMatrixTdEVD(d,e,n,3,z))
{
info=-3;
return;
}
//--- Generate
ArrayResize(x,n);
ArrayResize(w,n);
//--- calculation
for(i=1; i<=n; i++)
{
x[i-1]=d[i-1];
w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
}
}
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss-Lobatto quadrature |
//| formula |
//| The algorithm generates the N-point Gauss-Lobatto quadrature |
//| formula with weight function given by coefficients alpha and beta|
//| of a recurrence which generates a system of orthogonal |
//| polynomials. |
//| P-1(x) = 0 |
//| P0(x) = 1 |
//| Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x) |
//| and zeroth moment Mu0 |
//| Mu0 = integral(W(x)dx,a,b) |
//| INPUT PARAMETERS: |
//| Alpha - array[0..N-2], alpha coefficients |
//| Beta - array[0..N-2], beta coefficients. |
//| Zero-indexed element is not used, may be |
//| arbitrary. Beta[I]>0 |
//| Mu0 - zeroth moment of the weighting function. |
//| A - left boundary of the integration interval. |
//| B - right boundary of the integration interval. |
//| N - number of nodes of the quadrature formula, N>=3 |
//| (including the left and right boundary nodes). |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -2 Beta[i]<=0 |
//| * -1 incorrect N was passed |
//| * 1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLobattoRec(double &calpha[],double &cbeta[],
const double mu0,const double a,
const double b,int n,int &info,
double &x[],double &w[])
{
//--- create variables
int i=0;
double pim1a=0;
double pia=0;
double pim1b=0;
double pib=0;
double t=0;
double a11=0;
double a12=0;
double a21=0;
double a22=0;
double b1=0;
double b2=0;
double alph=0;
double bet=0;
//--- create arrays
double d[];
double e[];
//--- create matrix
CMatrixDouble z;
//--- creating copies of arrays
double alpha[];
double beta[];
ArrayCopy(alpha,calpha);
ArrayCopy(beta,cbeta);
//--- initialization
info=0;
//--- check
if(n<=2)
{
info=-1;
return;
}
info=1;
//--- Initialize,D[1:N+1],E[1:N]
n=n-2;
//--- allocation
ArrayResize(d,n+2);
ArrayResize(e,n+1);
//--- copy
for(i=1; i<=n+1; i++)
d[i-1]=alpha[i-1];
for(i=1; i<=n; i++)
{
//--- check
if(beta[i]<=0.0)
{
info=-2;
return;
}
e[i-1]=MathSqrt(beta[i]);
}
//--- Caclulate Pn(a),Pn+1(a),Pn(b),Pn+1(b)
beta[0]=0;
pim1a=0;
pia=1;
pim1b=0;
pib=1;
for(i=1; i<=n+1; i++)
{
//--- Pi(a)
t=(a-alpha[i-1])*pia-beta[i-1]*pim1a;
pim1a=pia;
pia=t;
//--- Pi(b)
t=(b-alpha[i-1])*pib-beta[i-1]*pim1b;
pim1b=pib;
pib=t;
}
//--- Calculate alpha'(n+1),beta'(n+1)
a11=pia;
a12=pim1a;
a21=pib;
a22=pim1b;
b1=a*pia;
b2=b*pib;
//--- check
if(MathAbs(a11)>MathAbs(a21))
{
//--- change values
a22=a22-a12*a21/a11;
b2=b2-b1*a21/a11;
bet=b2/a22;
alph=(b1-bet*a12)/a11;
}
else
{
//--- change values
a12=a12-a22*a11/a21;
b1=b1-b2*a11/a21;
bet=b1/a12;
alph=(b2-bet*a22)/a21;
}
//--- check
if(bet<0.0)
{
info=-3;
return;
}
//--- change values
d[n+1]=alph;
e[n]=MathSqrt(bet);
//--- EVD
if(!CEigenVDetect::SMatrixTdEVD(d,e,n+2,3,z))
{
info=-3;
return;
}
//--- Generate
ArrayResize(x,n+2);
ArrayResize(w,n+2);
//--- get result
for(i=1; i<=n+2; i++)
{
x[i-1]=d[i-1];
w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
}
}
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss-Radau quadrature |
//| formula |
//| The algorithm generates the N-point Gauss-Radau quadrature |
//| formula with weight function given by the coefficients alpha and |
//| beta of a recurrence which generates a system of orthogonal |
//| polynomials. |
//| P-1(x) = 0 |
//| P0(x) = 1 |
//| Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x) |
//| and zeroth moment Mu0 |
//| Mu0 = integral(W(x)dx,a,b) |
//| INPUT PARAMETERS: |
//| Alpha - array[0..N-2], alpha coefficients. |
//| Beta - array[0..N-1], beta coefficients |
//| Zero-indexed element is not used. |
//| Beta[I]>0 |
//| Mu0 - zeroth moment of the weighting function. |
//| A - left boundary of the integration interval. |
//| N - number of nodes of the quadrature formula, N>=2 |
//| (including the left boundary node). |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -2 Beta[i]<=0 |
//| * -1 incorrect N was passed |
//| * 1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussRadauRec(double &calpha[],double &cbeta[],
const double mu0,const double a,
int n,int &info,double &x[],
double &w[])
{
//--- create variables
int i=0;
double polim1=0;
double poli=0;
double t=0;
//--- create arrays
double d[];
double e[];
//--- create matrix
CMatrixDouble z;
//--- creating copies of arrays
double alpha[];
double beta[];
ArrayCopy(alpha,calpha);
ArrayCopy(beta,cbeta);
//--- initialization
info=0;
//--- check
if(n<2)
{
info=-1;
return;
}
info=1;
//--- Initialize,D[1:N],E[1:N]
n=n-1;
//--- allocation
ArrayResize(d,n+1);
ArrayResize(e,n);
for(i=1; i<=n; i++)
{
d[i-1]=alpha[i-1];
//--- check
if(beta[i]<=0.0)
{
info=-2;
return;
}
e[i-1]=MathSqrt(beta[i]);
}
//--- Caclulate Pn(a),Pn-1(a),and D[N+1]
beta[0]=0;
polim1=0;
poli=1;
//--- calculation
for(i=1; i<=n; i++)
{
t=(a-alpha[i-1])*poli-beta[i-1]*polim1;
polim1=poli;
poli=t;
}
d[n]=a-beta[n]*polim1/poli;
//--- EVD
if(!CEigenVDetect::SMatrixTdEVD(d,e,n+1,3,z))
{
info=-3;
return;
}
//--- Generate
ArrayResize(x,n+1);
ArrayResize(w,n+1);
//--- get result
for(i=1; i<=n+1; i++)
{
x[i-1]=d[i-1];
w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
}
}
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] |
//| with N nodes. |
//| INPUT PARAMETERS: |
//| N - number of nodes, >=1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. N is too |
//| large to obtain weights/nodes with |
//| high enough accuracy. Try to use |
//| multiple precision version. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLegendre(const int n,int &info,
double &x[],double &w[])
{
//--- create a variable
int i=0;
//--- create arrays
double alpha[];
double beta[];
//--- initialization
info=0;
//--- check
if(n<1)
{
info=-1;
return;
}
//--- allocation
ArrayResize(alpha,n);
ArrayResize(beta,n);
//--- initialization
for(i=0; i<=n-1; i++)
alpha[i]=0;
beta[0]=2;
for(i=1; i<=n-1; i++)
beta[i]=1/(4-1/CMath::Sqr(i));
//--- function call
GQGenerateRec(alpha,beta,beta[0],n,info,x,w);
//--- test basic properties to detect errors
if(info>0)
{
//--- check
if(x[0]<-1.0 || x[n-1]>1.0)
info=-4;
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] |
//| with weight function W(x)=Power(1-x,Alpha)*Power(1+x,Beta). |
//| INPUT PARAMETERS: |
//| N - number of nodes, >=1 |
//| Alpha - power-law coefficient, Alpha>-1 |
//| Beta - power-law coefficient, Beta>-1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. Alpha or |
//| Beta are too close to -1 to obtain |
//| weights/nodes with high enough |
//| accuracy, or, may be, N is too large.|
//| Try to use multiple precision version|
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N/Alpha/Beta was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussJacobi(const int n,const double alpha,
const double beta,int &info,
double &x[],double &w[])
{
//--- create variables
double alpha2=0;
double beta2=0;
double apb=0;
double t=0;
int i=0;
double s=0;
//--- create arrays
double a[];
double b[];
//--- initialization
info=0;
//--- check
if((n<1 || alpha<=-1.0) || beta<=-1.0)
{
info=-1;
return;
}
//--- allocation
ArrayResize(a,n);
ArrayResize(b,n);
//--- initialization
apb=alpha+beta;
a[0]=(beta-alpha)/(apb+2);
t=(apb+1)*MathLog(2)+CGammaFunc::LnGamma(alpha+1,s)+CGammaFunc::LnGamma(beta+1,s)-CGammaFunc::LnGamma(apb+2,s);
//--- check
if(t>MathLog(CMath::m_maxrealnumber))
{
info=-4;
return;
}
b[0]=MathExp(t);
//--- check
if(n>1)
{
//--- change values
alpha2=CMath::Sqr(alpha);
beta2=CMath::Sqr(beta);
a[1]=(beta2-alpha2)/((apb+2)*(apb+4));
b[1]=4*(alpha+1)*(beta+1)/((apb+3)*CMath::Sqr(apb+2));
//--- calculation
for(i=2; i<=n-1; i++)
{
a[i]=0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
b[i]=0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*CMath::Sqr(1+0.5*apb/i));
}
}
//--- function call
GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
if(info>0)
{
//--- check
if(x[0]<-1.0 || x[n-1]>1.0)
info=-4;
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) |
//| with weight function W(x)=Power(x,Alpha)*Exp(-x) |
//| INPUT PARAMETERS: |
//| N - number of nodes, >=1 |
//| Alpha - power-law coefficient, Alpha>-1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. Alpha is |
//| too close to -1 to obtain |
//| weights/nodes with high enough |
//| accuracy or, may be, N is too large.|
//| Try to use multiple precision |
//| version. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N/Alpha was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLaguerre(const int n,const double alpha,
int &info,double &x[],
double &w[])
{
//--- create arrays
double a[];
double b[];
//--- create variables
double t=0;
int i=0;
double s=0;
//--- initialization
info=0;
//--- check
if(n<1 || alpha<=-1.0)
{
info=-1;
return;
}
//--- allocation
ArrayResize(a,n);
ArrayResize(b,n);
//--- initialization
a[0]=alpha+1;
t=CGammaFunc::LnGamma(alpha+1,s);
//--- check
if(t>=MathLog(CMath::m_maxrealnumber))
{
info=-4;
return;
}
b[0]=MathExp(t);
//--- check
if(n>1)
{
//--- calculation
for(i=1; i<=n-1; i++)
{
a[i]=2*i+alpha+1;
b[i]=i*(i+alpha);
}
}
//--- function call
GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
if(info>0)
{
//--- check
if(x[0]<0.0)
info=-4;
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Hermite quadrature on |
//| (-inf,+inf) with weight function W(x)=Exp(-x*x) |
//| INPUT PARAMETERS: |
//| N - number of nodes, >=1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. May be, N |
//| is too large. Try to use multiple |
//| precision version. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N/Alpha was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| W - array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussHermite(const int n,int &info,
double &x[],double &w[])
{
//--- create a variable
int i=0;
//--- create arrays
double a[];
double b[];
//--- initialization
info=0;
//--- check
if(n<1)
{
info=-1;
return;
}
//--- allocation
ArrayResize(a,n);
ArrayResize(b,n);
//--- initialization
for(i=0; i<=n-1; i++)
a[i]=0;
b[0]=MathSqrt(4*MathArctan(1));
//--- check
if(n>1)
{
for(i=1; i<=n-1; i++)
b[i]=0.5*i;
}
//--- function call
GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
if(info>0)
{
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Gauss-Kronrod quadrature formula |
//+------------------------------------------------------------------+
class CGaussKronrodQ
{
public:
static void GKQGenerateRec(double &calpha[],double &cbeta[],const double mu0,int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
static void GKQGenerateGaussLegendre(const int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
static void GKQGenerateGaussJacobi(const int n,const double alpha,const double beta,int &info,double &x[],double &wkronrod[],double &wgauss[]);
static void GKQLegendreCalc(const int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
static void GKQLegendreTbl(const int n,double &x[],double &wkronrod[],double &wgauss[],double &eps);
};
//+------------------------------------------------------------------+
//| Computation of nodes and weights of a Gauss-Kronrod quadrature |
//| formula |
//| The algorithm generates the N-point Gauss-Kronrod quadrature |
//| formula with weight function given by coefficients alpha and beta|
//| of a recurrence relation which generates a system of orthogonal |
//| polynomials: |
//| P-1(x) = 0 |
//| P0(x) = 1 |
//| Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x) |
//| and zero moment Mu0 |
//| Mu0 = integral(W(x)dx,a,b) |
//| INPUT PARAMETERS: |
//| Alpha - alpha coefficients, array[0..floor(3*K/2)]. |
//| Beta - beta coefficients, array[0..ceil(3*K/2)]. |
//| Beta[0] is not used and may be arbitrary. |
//| Beta[I]>0. |
//| Mu0 - zeroth moment of the weight function. |
//| N - number of nodes of the Gauss-Kronrod |
//| quadrature formula, |
//| N >= 3, |
//| N = 2*K+1. |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -5 no real and positive Gauss-Kronrod |
//| formula can be created for such a |
//| weight function with a given number |
//| of nodes. |
//| * -4 N is too large, task may be ill |
//| conditioned - x[i]=x[i+1] found. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -2 Beta[i]<=0 |
//| * -1 incorrect N was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| in ascending order. |
//| WKronrod - array[0..N-1] - Kronrod weights |
//| WGauss - array[0..N-1] - Gauss weights (interleaved |
//| with zeros corresponding to extended Kronrod |
//| nodes). |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateRec(double &calpha[],double &cbeta[],
const double mu0,int n,int &info,
double &x[],double &wkronrod[],
double &wgauss[])
{
//--- create variables
int i=0;
int j=0;
int wlen=0;
int woffs=0;
double u=0;
int m=0;
int l=0;
int k=0;
int i_=0;
//--- create arrays
double ta[];
double t[];
double s[];
double xgtmp[];
double wgtmp[];
//--- creating copies of arrays
double alpha[];
double beta[];
ArrayCopy(alpha,calpha);
ArrayCopy(beta,cbeta);
//--- initialization
info=0;
//--- check
if(n%2!=1 || n<3)
{
info=-1;
return;
}
for(i=0; i<=(int)MathCeil((double)(3*(n/2))/2.0); i++)
{
//--- check
if(beta[i]<=0.0)
{
info=-2;
return;
}
}
info=1;
//--- from external conventions about N/Beta/Mu0 to internal
n=n/2;
beta[0]=mu0;
//--- Calculate Gauss nodes/weights,save them for later processing
CGaussQ::GQGenerateRec(alpha,beta,mu0,n,info,xgtmp,wgtmp);
//--- check
if(info<0)
return;
//--- Resize:
//--- * A from 0..floor(3*n/2) to 0..2*n
//--- * B from 0..ceil(3*n/2) to 0..2*n
ArrayResize(ta,(int)MathFloor((double)(3*n)/2.0)+1);
//--- copy
for(i_=0; i_<=(int)MathFloor((double)(3*n)/2.0); i_++)
ta[i_]=alpha[i_];
//--- allocation
ArrayResize(alpha,2*n+1);
//--- copy
for(i_=0; i_<=(int)MathFloor((double)(3*n)/2.0); i_++)
alpha[i_]=ta[i_];
//--- initialization
for(i=(int)MathFloor((double)(3*n)/2.0)+1; i<=2*n; i++)
alpha[i]=0;
//--- allocation
ArrayResize(ta,(int)MathCeil((double)(3*n)/2.0)+1);
//--- copy
for(i_=0; i_<=(int)MathCeil((double)(3*n)/2.0); i_++)
ta[i_]=beta[i_];
//--- allocation
ArrayResize(beta,2*n+1);
//--- copy
for(i_=0; i_<=(int)MathCeil((double)(3*n)/2.0); i_++)
beta[i_]=ta[i_];
//--- initialization
for(i=(int)MathCeil((double)(3*n)/2.0)+1; i<=2*n; i++)
beta[i]=0;
//--- Initialize T,S
wlen=2+n/2;
//--- allocation
ArrayResize(t,wlen);
ArrayResize(s,wlen);
ArrayResize(ta,wlen);
woffs=1;
for(i=0; i<=wlen-1; i++)
{
t[i]=0;
s[i]=0;
}
//--- Algorithm from Dirk P. Laurie,"Calculation of Gauss-Kronrod quadrature rules",1997.
t[woffs+0]=beta[n+1];
for(m=0; m<=n-2; m++)
{
u=0;
//--- calculation
for(k=(m+1)/2; k>=0; k--)
{
l=m-k;
u=u+(alpha[k+n+1]-alpha[l])*t[woffs+k]+beta[k+n+1]*s[woffs+k-1]-beta[l]*s[woffs+k];
s[woffs+k]=u;
}
//--- copy
for(i_=0; i_<=wlen-1; i_++)
ta[i_]=t[i_];
//--- copy
for(i_=0; i_<=wlen-1; i_++)
t[i_]=s[i_];
//--- copy
for(i_=0; i_<=wlen-1; i_++)
s[i_]=ta[i_];
}
//--- copy
for(j=n/2; j>=0; j--)
s[woffs+j]=s[woffs+j-1];
//--- cycle
for(m=n-1; m<=2*n-3; m++)
{
u=0;
//--- calculation
for(k=m+1-n; k<=(m-1)/2; k++)
{
l=m-k;
j=n-1-l;
u=u-(alpha[k+n+1]-alpha[l])*t[woffs+j]-beta[k+n+1]*s[woffs+j]+beta[l]*s[woffs+j+1];
s[woffs+j]=u;
}
//--- check
if(m%2==0)
{
k=m/2;
alpha[k+n+1]=alpha[k]+(s[woffs+j]-beta[k+n+1]*s[woffs+j+1])/t[woffs+j+1];
}
else
{
k=(m+1)/2;
beta[k+n+1]=s[woffs+j]/s[woffs+j+1];
}
//--- copy
for(i_=0; i_<=wlen-1; i_++)
ta[i_]=t[i_];
//--- copy
for(i_=0; i_<=wlen-1; i_++)
t[i_]=s[i_];
//--- copy
for(i_=0; i_<=wlen-1; i_++)
s[i_]=ta[i_];
}
//--- change value
alpha[2*n]=alpha[n-1]-beta[2*n]*s[woffs+0]/t[woffs+0];
//--- calculation of Kronrod nodes and weights,unpacking of Gauss weights
CGaussQ::GQGenerateRec(alpha,beta,mu0,2*n+1,info,x,wkronrod);
//--- check
if(info==-2)
info=-5;
//--- check
if(info<0)
return;
for(i=0; i<=2*n-1; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
//--- check
if(info<0)
return;
//--- allocation
ArrayResize(wgauss,2*n+1);
//--- get result
for(i=0; i<=2*n; i++)
wgauss[i]=0;
for(i=0; i<=n-1; i++)
wgauss[2*i+1]=wgtmp[i];
}
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre |
//| quadrature with N points. |
//| GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed |
//| table) is used depending on machine precision and number of |
//| nodes. |
//| INPUT PARAMETERS: |
//| N - number of Kronrod nodes, must be odd number, |
//| >=3. |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. N is too |
//| obtain large to weights/nodes with |
//| high enough accuracy. Try to use |
//| multiple precision version. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| ordered in ascending order. |
//| WKronrod - array[0..N-1] - Kronrod weights |
//| WGauss - array[0..N-1] - Gauss weights (interleaved |
//| with zeros corresponding to extended Kronrod |
//| nodes). |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateGaussLegendre(const int n,int &info,
double &x[],double &wkronrod[],
double &wgauss[])
{
//--- create a variable
double eps=0;
//--- initialization
info=0;
//--- check
if(CMath::m_machineepsilon>1.0E-32 && (((((n==15 || n==21) || n==31) || n==41) || n==51) || n==61))
{
info=1;
GKQLegendreTbl(n,x,wkronrod,wgauss,eps);
}
else
GKQLegendreCalc(n,info,x,wkronrod,wgauss);
}
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi |
//| quadrature on [-1,1] with weight function |
//| W(x)=Power(1-x,Alpha)*Power(1+x,Beta). |
//| INPUT PARAMETERS: |
//| N - number of Kronrod nodes, must be odd number, |
//| >=3. |
//| Alpha - power-law coefficient, Alpha>-1 |
//| Beta - power-law coefficient, Beta>-1 |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -5 no real and positive Gauss-Kronrod |
//| formula can be created for such a |
//| weight function with a given number |
//| of nodes. |
//| * -4 an error was detected when |
//| calculating weights/nodes. Alpha or |
//| Beta are too close to -1 to obtain |
//| weights/nodes with high enough |
//| accuracy, or, may be, N is too large.|
//| Try to use multiple precision version|
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N was passed |
//| * +1 OK |
//| * +2 OK, but quadrature rule have exterior|
//| nodes, x[0]<-1 or x[n-1]>+1 |
//| X - array[0..N-1] - array of quadrature nodes, |
//| ordered in ascending order. |
//| WKronrod - array[0..N-1] - Kronrod weights |
//| WGauss - array[0..N-1] - Gauss weights (interleaved |
//| with zeros corresponding to extended Kronrod |
//| nodes). |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateGaussJacobi(const int n,const double alpha,
const double beta,int &info,
double &x[],double &wkronrod[],
double &wgauss[])
{
//--- create variables
int clen=0;
double alpha2=0;
double beta2=0;
double apb=0;
double t=0;
int i=0;
double s=0;
//--- create arrays
double a[];
double b[];
//--- initialization
info=0;
//--- check
if(n%2!=1 || n<3)
{
info=-1;
return;
}
//--- check
if(alpha<=-1.0 || beta<=-1.0)
{
info=-1;
return;
}
//--- change value
clen=(int)MathCeil((double)(3*(n/2))/2.0)+1;
//--- allocation
ArrayResize(a,clen);
ArrayResize(b,clen);
//--- initialization
for(i=0; i<=clen-1; i++)
a[i]=0;
apb=alpha+beta;
a[0]=(beta-alpha)/(apb+2);
t=(apb+1)*MathLog(2)+CGammaFunc::LnGamma(alpha+1,s)+CGammaFunc::LnGamma(beta+1,s)-CGammaFunc::LnGamma(apb+2,s);
//--- check
if(t>MathLog(CMath::m_maxrealnumber))
{
info=-4;
return;
}
b[0]=MathExp(t);
//--- check
if(clen>1)
{
//--- change values
alpha2=CMath::Sqr(alpha);
beta2=CMath::Sqr(beta);
a[1]=(beta2-alpha2)/((apb+2)*(apb+4));
b[1]=4*(alpha+1)*(beta+1)/((apb+3)*CMath::Sqr(apb+2));
//--- calculation
for(i=2; i<=clen-1; i++)
{
a[i]=0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
b[i]=0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*CMath::Sqr(1+0.5*apb/i));
}
}
//--- function call
GKQGenerateRec(a,b,b[0],n,info,x,wkronrod,wgauss);
//--- test basic properties to detect errors
if(info>0)
{
//--- check
if(x[0]<-1.0 || x[n-1]>1.0)
info=2;
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes for quadrature with N |
//| points. |
//| Reduction to tridiagonal eigenproblem is used. |
//| INPUT PARAMETERS: |
//| N - number of Kronrod nodes, must be odd number, |
//| >=3. |
//| OUTPUT PARAMETERS: |
//| Info - error code: |
//| * -4 an error was detected when |
//| calculating weights/nodes. N is too |
//| large to obtain weights/nodes with |
//| high enough accuracy. |
//| Try to use multiple precision |
//| version. |
//| * -3 internal eigenproblem solver hasn't |
//| converged |
//| * -1 incorrect N was passed |
//| * +1 OK |
//| X - array[0..N-1] - array of quadrature nodes, |
//| ordered in ascending order. |
//| WKronrod - array[0..N-1] - Kronrod weights |
//| WGauss - array[0..N-1] - Gauss weights (interleaved |
//| with zeros corresponding to extended Kronrod |
//| nodes). |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQLegendreCalc(const int n,int &info,double &x[],
double &wkronrod[],double &wgauss[])
{
//--- create variables
int alen=0;
int blen=0;
double mu0=0;
int k=0;
int i=0;
//--- create arrays
double alpha[];
double beta[];
//--- initialization
info=0;
//--- check
if(n%2!=1 || n<3)
{
info=-1;
return;
}
//--- initialization
mu0=2;
alen=(int)MathFloor((double)(3*(n/2))/2.0)+1;
blen=(int)MathCeil((double)(3*(n/2))/2.0)+1;
//--- allocation
ArrayResize(alpha,alen);
ArrayResize(beta,blen);
//--- initialization
for(k=0; k<=alen-1; k++)
alpha[k]=0;
beta[0]=2;
for(k=1; k<=blen-1; k++)
beta[k]=1/(4-1/CMath::Sqr(k));
//--- function call
GKQGenerateRec(alpha,beta,mu0,n,info,x,wkronrod,wgauss);
//--- test basic properties to detect errors
if(info>0)
{
//--- check
if(x[0]<-1.0 || x[n-1]>1.0)
info=-4;
for(i=0; i<=n-2; i++)
{
//--- check
if(x[i]>=x[i+1])
info=-4;
}
}
}
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes for quadrature with N |
//| points using pre-calculated table. Nodes/weights were computed |
//| with accuracy up to 1.0E-32 (if MPFR version of ALGLIB is used). |
//| In standard double precision accuracy reduces to something about |
//| 2.0E-16 (depending on your compiler's handling of long floating |
//| point constants). |
//| INPUT PARAMETERS: |
//| N - number of Kronrod nodes. |
//| N can be 15, 21, 31, 41, 51, 61. |
//| OUTPUT PARAMETERS: |
//| X - array[0..N-1] - array of quadrature nodes, |
//| ordered in ascending order. |
//| WKronrod - array[0..N-1] - Kronrod weights |
//| WGauss - array[0..N-1] - Gauss weights (interleaved |
//| with zeros corresponding to extended Kronrod |
//| nodes). |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQLegendreTbl(const int n,double &x[],
double &wkronrod[],double &wgauss[],
double &eps)
{
//--- create variables
int i=0;
int ng=0;
double tmp=0;
//--- create arrays
int p1[];
int p2[];
//--- initialization
eps=0;
//--- these initializers are not really necessary,
//--- but without them compiler complains about uninitialized locals
ng=0;
//--- Process
if(!CAp::Assert(n==15 || n==21 || n==31 || n==41 || n==51 || n==61,__FUNCTION__+": incorrect N!"))
return;
//--- allocation
ArrayResize(x,n);
ArrayResize(wkronrod,n);
ArrayResize(wgauss,n);
//--- initialization
for(i=0; i<=n-1; i++)
{
x[i]=0;
wkronrod[i]=0;
wgauss[i]=0;
}
eps=MathMax(CMath::m_machineepsilon,1.0E-32);
//--- check
if(n==15)
{
//--- change values
ng=4;
wgauss[0]=0.129484966168869693270611432679082;
wgauss[1]=0.279705391489276667901467771423780;
wgauss[2]=0.381830050505118944950369775488975;
wgauss[3]=0.417959183673469387755102040816327;
x[0]=0.991455371120812639206854697526329;
x[1]=0.949107912342758524526189684047851;
x[2]=0.864864423359769072789712788640926;
x[3]=0.741531185599394439863864773280788;
x[4]=0.586087235467691130294144838258730;
x[5]=0.405845151377397166906606412076961;
x[6]=0.207784955007898467600689403773245;
x[7]=0.000000000000000000000000000000000;
wkronrod[0]=0.022935322010529224963732008058970;
wkronrod[1]=0.063092092629978553290700663189204;
wkronrod[2]=0.104790010322250183839876322541518;
wkronrod[3]=0.140653259715525918745189590510238;
wkronrod[4]=0.169004726639267902826583426598550;
wkronrod[5]=0.190350578064785409913256402421014;
wkronrod[6]=0.204432940075298892414161999234649;
wkronrod[7]=0.209482141084727828012999174891714;
}
//--- check
if(n==21)
{
//--- change values
ng=5;
wgauss[0]=0.066671344308688137593568809893332;
wgauss[1]=0.149451349150580593145776339657697;
wgauss[2]=0.219086362515982043995534934228163;
wgauss[3]=0.269266719309996355091226921569469;
wgauss[4]=0.295524224714752870173892994651338;
x[0]=0.995657163025808080735527280689003;
x[1]=0.973906528517171720077964012084452;
x[2]=0.930157491355708226001207180059508;
x[3]=0.865063366688984510732096688423493;
x[4]=0.780817726586416897063717578345042;
x[5]=0.679409568299024406234327365114874;
x[6]=0.562757134668604683339000099272694;
x[7]=0.433395394129247190799265943165784;
x[8]=0.294392862701460198131126603103866;
x[9]=0.148874338981631210884826001129720;
x[10]=0.000000000000000000000000000000000;
wkronrod[0]=0.011694638867371874278064396062192;
wkronrod[1]=0.032558162307964727478818972459390;
wkronrod[2]=0.054755896574351996031381300244580;
wkronrod[3]=0.075039674810919952767043140916190;
wkronrod[4]=0.093125454583697605535065465083366;
wkronrod[5]=0.109387158802297641899210590325805;
wkronrod[6]=0.123491976262065851077958109831074;
wkronrod[7]=0.134709217311473325928054001771707;
wkronrod[8]=0.142775938577060080797094273138717;
wkronrod[9]=0.147739104901338491374841515972068;
wkronrod[10]=0.149445554002916905664936468389821;
}
//--- check
if(n==31)
{
//--- change values
ng=8;
wgauss[0]=0.030753241996117268354628393577204;
wgauss[1]=0.070366047488108124709267416450667;
wgauss[2]=0.107159220467171935011869546685869;
wgauss[3]=0.139570677926154314447804794511028;
wgauss[4]=0.166269205816993933553200860481209;
wgauss[5]=0.186161000015562211026800561866423;
wgauss[6]=0.198431485327111576456118326443839;
wgauss[7]=0.202578241925561272880620199967519;
x[0]=0.998002298693397060285172840152271;
x[1]=0.987992518020485428489565718586613;
x[2]=0.967739075679139134257347978784337;
x[3]=0.937273392400705904307758947710209;
x[4]=0.897264532344081900882509656454496;
x[5]=0.848206583410427216200648320774217;
x[6]=0.790418501442465932967649294817947;
x[7]=0.724417731360170047416186054613938;
x[8]=0.650996741297416970533735895313275;
x[9]=0.570972172608538847537226737253911;
x[10]=0.485081863640239680693655740232351;
x[11]=0.394151347077563369897207370981045;
x[12]=0.299180007153168812166780024266389;
x[13]=0.201194093997434522300628303394596;
x[14]=0.101142066918717499027074231447392;
x[15]=0.000000000000000000000000000000000;
wkronrod[0]=0.005377479872923348987792051430128;
wkronrod[1]=0.015007947329316122538374763075807;
wkronrod[2]=0.025460847326715320186874001019653;
wkronrod[3]=0.035346360791375846222037948478360;
wkronrod[4]=0.044589751324764876608227299373280;
wkronrod[5]=0.053481524690928087265343147239430;
wkronrod[6]=0.062009567800670640285139230960803;
wkronrod[7]=0.069854121318728258709520077099147;
wkronrod[8]=0.076849680757720378894432777482659;
wkronrod[9]=0.083080502823133021038289247286104;
wkronrod[10]=0.088564443056211770647275443693774;
wkronrod[11]=0.093126598170825321225486872747346;
wkronrod[12]=0.096642726983623678505179907627589;
wkronrod[13]=0.099173598721791959332393173484603;
wkronrod[14]=0.100769845523875595044946662617570;
wkronrod[15]=0.101330007014791549017374792767493;
}
//--- check
if(n==41)
{
//--- change values
ng=10;
wgauss[0]=0.017614007139152118311861962351853;
wgauss[1]=0.040601429800386941331039952274932;
wgauss[2]=0.062672048334109063569506535187042;
wgauss[3]=0.083276741576704748724758143222046;
wgauss[4]=0.101930119817240435036750135480350;
wgauss[5]=0.118194531961518417312377377711382;
wgauss[6]=0.131688638449176626898494499748163;
wgauss[7]=0.142096109318382051329298325067165;
wgauss[8]=0.149172986472603746787828737001969;
wgauss[9]=0.152753387130725850698084331955098;
x[0]=0.998859031588277663838315576545863;
x[1]=0.993128599185094924786122388471320;
x[2]=0.981507877450250259193342994720217;
x[3]=0.963971927277913791267666131197277;
x[4]=0.940822633831754753519982722212443;
x[5]=0.912234428251325905867752441203298;
x[6]=0.878276811252281976077442995113078;
x[7]=0.839116971822218823394529061701521;
x[8]=0.795041428837551198350638833272788;
x[9]=0.746331906460150792614305070355642;
x[10]=0.693237656334751384805490711845932;
x[11]=0.636053680726515025452836696226286;
x[12]=0.575140446819710315342946036586425;
x[13]=0.510867001950827098004364050955251;
x[14]=0.443593175238725103199992213492640;
x[15]=0.373706088715419560672548177024927;
x[16]=0.301627868114913004320555356858592;
x[17]=0.227785851141645078080496195368575;
x[18]=0.152605465240922675505220241022678;
x[19]=0.076526521133497333754640409398838;
x[20]=0.000000000000000000000000000000000;
wkronrod[0]=0.003073583718520531501218293246031;
wkronrod[1]=0.008600269855642942198661787950102;
wkronrod[2]=0.014626169256971252983787960308868;
wkronrod[3]=0.020388373461266523598010231432755;
wkronrod[4]=0.025882133604951158834505067096153;
wkronrod[5]=0.031287306777032798958543119323801;
wkronrod[6]=0.036600169758200798030557240707211;
wkronrod[7]=0.041668873327973686263788305936895;
wkronrod[8]=0.046434821867497674720231880926108;
wkronrod[9]=0.050944573923728691932707670050345;
wkronrod[10]=0.055195105348285994744832372419777;
wkronrod[11]=0.059111400880639572374967220648594;
wkronrod[12]=0.062653237554781168025870122174255;
wkronrod[13]=0.065834597133618422111563556969398;
wkronrod[14]=0.068648672928521619345623411885368;
wkronrod[15]=0.071054423553444068305790361723210;
wkronrod[16]=0.073030690332786667495189417658913;
wkronrod[17]=0.074582875400499188986581418362488;
wkronrod[18]=0.075704497684556674659542775376617;
wkronrod[19]=0.076377867672080736705502835038061;
wkronrod[20]=0.076600711917999656445049901530102;
}
//--- check
if(n==51)
{
//--- change values
ng=13;
wgauss[0]=0.011393798501026287947902964113235;
wgauss[1]=0.026354986615032137261901815295299;
wgauss[2]=0.040939156701306312655623487711646;
wgauss[3]=0.054904695975835191925936891540473;
wgauss[4]=0.068038333812356917207187185656708;
wgauss[5]=0.080140700335001018013234959669111;
wgauss[6]=0.091028261982963649811497220702892;
wgauss[7]=0.100535949067050644202206890392686;
wgauss[8]=0.108519624474263653116093957050117;
wgauss[9]=0.114858259145711648339325545869556;
wgauss[10]=0.119455763535784772228178126512901;
wgauss[11]=0.122242442990310041688959518945852;
wgauss[12]=0.123176053726715451203902873079050;
x[0]=0.999262104992609834193457486540341;
x[1]=0.995556969790498097908784946893902;
x[2]=0.988035794534077247637331014577406;
x[3]=0.976663921459517511498315386479594;
x[4]=0.961614986425842512418130033660167;
x[5]=0.942974571228974339414011169658471;
x[6]=0.920747115281701561746346084546331;
x[7]=0.894991997878275368851042006782805;
x[8]=0.865847065293275595448996969588340;
x[9]=0.833442628760834001421021108693570;
x[10]=0.797873797998500059410410904994307;
x[11]=0.759259263037357630577282865204361;
x[12]=0.717766406813084388186654079773298;
x[13]=0.673566368473468364485120633247622;
x[14]=0.626810099010317412788122681624518;
x[15]=0.577662930241222967723689841612654;
x[16]=0.526325284334719182599623778158010;
x[17]=0.473002731445714960522182115009192;
x[18]=0.417885382193037748851814394594572;
x[19]=0.361172305809387837735821730127641;
x[20]=0.303089538931107830167478909980339;
x[21]=0.243866883720988432045190362797452;
x[22]=0.183718939421048892015969888759528;
x[23]=0.122864692610710396387359818808037;
x[24]=0.061544483005685078886546392366797;
x[25]=0.000000000000000000000000000000000;
wkronrod[0]=0.001987383892330315926507851882843;
wkronrod[1]=0.005561932135356713758040236901066;
wkronrod[2]=0.009473973386174151607207710523655;
wkronrod[3]=0.013236229195571674813656405846976;
wkronrod[4]=0.016847817709128298231516667536336;
wkronrod[5]=0.020435371145882835456568292235939;
wkronrod[6]=0.024009945606953216220092489164881;
wkronrod[7]=0.027475317587851737802948455517811;
wkronrod[8]=0.030792300167387488891109020215229;
wkronrod[9]=0.034002130274329337836748795229551;
wkronrod[10]=0.037116271483415543560330625367620;
wkronrod[11]=0.040083825504032382074839284467076;
wkronrod[12]=0.042872845020170049476895792439495;
wkronrod[13]=0.045502913049921788909870584752660;
wkronrod[14]=0.047982537138836713906392255756915;
wkronrod[15]=0.050277679080715671963325259433440;
wkronrod[16]=0.052362885806407475864366712137873;
wkronrod[17]=0.054251129888545490144543370459876;
wkronrod[18]=0.055950811220412317308240686382747;
wkronrod[19]=0.057437116361567832853582693939506;
wkronrod[20]=0.058689680022394207961974175856788;
wkronrod[21]=0.059720340324174059979099291932562;
wkronrod[22]=0.060539455376045862945360267517565;
wkronrod[23]=0.061128509717053048305859030416293;
wkronrod[24]=0.061471189871425316661544131965264;
wkronrod[25]=0.061580818067832935078759824240055;
}
//--- check
if(n==61)
{
//--- change values
ng=15;
wgauss[0]=0.007968192496166605615465883474674;
wgauss[1]=0.018466468311090959142302131912047;
wgauss[2]=0.028784707883323369349719179611292;
wgauss[3]=0.038799192569627049596801936446348;
wgauss[4]=0.048402672830594052902938140422808;
wgauss[5]=0.057493156217619066481721689402056;
wgauss[6]=0.065974229882180495128128515115962;
wgauss[7]=0.073755974737705206268243850022191;
wgauss[8]=0.080755895229420215354694938460530;
wgauss[9]=0.086899787201082979802387530715126;
wgauss[10]=0.092122522237786128717632707087619;
wgauss[11]=0.096368737174644259639468626351810;
wgauss[12]=0.099593420586795267062780282103569;
wgauss[13]=0.101762389748405504596428952168554;
wgauss[14]=0.102852652893558840341285636705415;
x[0]=0.999484410050490637571325895705811;
x[1]=0.996893484074649540271630050918695;
x[2]=0.991630996870404594858628366109486;
x[3]=0.983668123279747209970032581605663;
x[4]=0.973116322501126268374693868423707;
x[5]=0.960021864968307512216871025581798;
x[6]=0.944374444748559979415831324037439;
x[7]=0.926200047429274325879324277080474;
x[8]=0.905573307699907798546522558925958;
x[9]=0.882560535792052681543116462530226;
x[10]=0.857205233546061098958658510658944;
x[11]=0.829565762382768397442898119732502;
x[12]=0.799727835821839083013668942322683;
x[13]=0.767777432104826194917977340974503;
x[14]=0.733790062453226804726171131369528;
x[15]=0.697850494793315796932292388026640;
x[16]=0.660061064126626961370053668149271;
x[17]=0.620526182989242861140477556431189;
x[18]=0.579345235826361691756024932172540;
x[19]=0.536624148142019899264169793311073;
x[20]=0.492480467861778574993693061207709;
x[21]=0.447033769538089176780609900322854;
x[22]=0.400401254830394392535476211542661;
x[23]=0.352704725530878113471037207089374;
x[24]=0.304073202273625077372677107199257;
x[25]=0.254636926167889846439805129817805;
x[26]=0.204525116682309891438957671002025;
x[27]=0.153869913608583546963794672743256;
x[28]=0.102806937966737030147096751318001;
x[29]=0.051471842555317695833025213166723;
x[30]=0.000000000000000000000000000000000;
wkronrod[0]=0.001389013698677007624551591226760;
wkronrod[1]=0.003890461127099884051267201844516;
wkronrod[2]=0.006630703915931292173319826369750;
wkronrod[3]=0.009273279659517763428441146892024;
wkronrod[4]=0.011823015253496341742232898853251;
wkronrod[5]=0.014369729507045804812451432443580;
wkronrod[6]=0.016920889189053272627572289420322;
wkronrod[7]=0.019414141193942381173408951050128;
wkronrod[8]=0.021828035821609192297167485738339;
wkronrod[9]=0.024191162078080601365686370725232;
wkronrod[10]=0.026509954882333101610601709335075;
wkronrod[11]=0.028754048765041292843978785354334;
wkronrod[12]=0.030907257562387762472884252943092;
wkronrod[13]=0.032981447057483726031814191016854;
wkronrod[14]=0.034979338028060024137499670731468;
wkronrod[15]=0.036882364651821229223911065617136;
wkronrod[16]=0.038678945624727592950348651532281;
wkronrod[17]=0.040374538951535959111995279752468;
wkronrod[18]=0.041969810215164246147147541285970;
wkronrod[19]=0.043452539701356069316831728117073;
wkronrod[20]=0.044814800133162663192355551616723;
wkronrod[21]=0.046059238271006988116271735559374;
wkronrod[22]=0.047185546569299153945261478181099;
wkronrod[23]=0.048185861757087129140779492298305;
wkronrod[24]=0.049055434555029778887528165367238;
wkronrod[25]=0.049795683427074206357811569379942;
wkronrod[26]=0.050405921402782346840893085653585;
wkronrod[27]=0.050881795898749606492297473049805;
wkronrod[28]=0.051221547849258772170656282604944;
wkronrod[29]=0.051426128537459025933862879215781;
wkronrod[30]=0.051494729429451567558340433647099;
}
//--- copy nodes
for(i=n-1; i>=n/2; i--)
x[i]=-x[n-1-i];
//--- copy Kronrod weights
for(i=n-1; i>=n/2; i--)
wkronrod[i]=wkronrod[n-1-i];
//--- copy Gauss weights
for(i=ng-1; i>=0; i--)
{
wgauss[n-2-2*i]=wgauss[i];
wgauss[1+2*i]=wgauss[i];
}
for(i=0; i<=n/2; i++)
wgauss[2*i]=0;
//--- reorder
CTSort::TagSort(x,n,p1,p2);
for(i=0; i<=n-1; i++)
{
//--- swap
tmp=wkronrod[i];
wkronrod[i]=wkronrod[p2[i]];
wkronrod[p2[i]]=tmp;
tmp=wgauss[i];
wgauss[i]=wgauss[p2[i]];
wgauss[p2[i]]=tmp;
}
}
//+------------------------------------------------------------------+
//| Integration report: |
//| * TerminationType = completetion code: |
//| * -5 non-convergence of Gauss-Kronrod nodes |
//| calculation subroutine. |
//| * -1 incorrect parameters were specified |
//| * 1 OK |
//| * Rep.NFEV countains number of function calculations |
//| * Rep.NIntervals contains number of intervals [a,b] |
//| was partitioned into. |
//+------------------------------------------------------------------+
class CAutoGKReport
{
public:
//--- variables
int m_terminationtype;
int m_nfev;
int m_nintervals;
//--- constructor, destructor
CAutoGKReport(void) { ZeroMemory(this); }
~CAutoGKReport(void) {}
//--- copy
void Copy(CAutoGKReport &obj);
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CAutoGKReport::Copy(CAutoGKReport &obj)
{
//--- copy variables
m_terminationtype=obj.m_terminationtype;
m_nfev=obj.m_nfev;
m_nintervals=obj.m_nintervals;
}
//+------------------------------------------------------------------+
//| Integration report: |
//| * TerminationType = completetion code: |
//| * -5 non-convergence of Gauss-Kronrod nodes |
//| calculation subroutine. |
//| * -1 incorrect parameters were specified |
//| * 1 OK |
//| * Rep.NFEV countains number of function calculations |
//| * Rep.NIntervals contains number of intervals [a,b] |
//| was partitioned into. |
//+------------------------------------------------------------------+
class CAutoGKReportShell
{
private:
CAutoGKReport m_innerobj;
public:
//--- constructors, destructor
CAutoGKReportShell(void) {}
CAutoGKReportShell(CAutoGKReport &obj) { m_innerobj.Copy(obj); }
~CAutoGKReportShell(void) {}
//--- methods
int GetTerminationType(void);
void SetTerminationType(const int i);
int GetNFev(void);
void SetNFev(const int i);
int GetNIntervals(void);
void SetNIntervals(const int i);
CAutoGKReport *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable nfev |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetNFev(void)
{
return(m_innerobj.m_nfev);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable nfev |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetNFev(const int i)
{
m_innerobj.m_nfev=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable terminationtype |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetTerminationType(void)
{
return(m_innerobj.m_terminationtype);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable terminationtype |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetTerminationType(const int i)
{
m_innerobj.m_terminationtype=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable nintervals |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetNIntervals(void)
{
return(m_innerobj.m_nintervals);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable nintervals |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetNIntervals(const int i)
{
m_innerobj.m_nintervals=i;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CAutoGKReport *CAutoGKReportShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Auxiliary class for CAutoGK |
//+------------------------------------------------------------------+
struct CAutoGKInternalState
{
//--- variables
double m_a;
double m_b;
double m_eps;
double m_xwidth;
double m_x;
double m_f;
int m_info;
double m_r;
int m_heapsize;
int m_heapwidth;
int m_heapused;
double m_sumerr;
double m_sumabs;
int m_n;
RCommState m_rstate;
//--- arrays
double m_qn[];
double m_wg[];
double m_wk[];
double m_wr[];
//--- matrix
CMatrixDouble m_heap;
//--- constructor, destructor
CAutoGKInternalState(void);
~CAutoGKInternalState(void) {}
//--- copy
void Copy(CAutoGKInternalState &obj);
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CAutoGKInternalState::CAutoGKInternalState(void)
{
m_a=0;
m_b=0;
m_eps=0;
m_xwidth=0;
m_x=0;
m_f=0;
m_info=0;
m_r=0;
m_heapsize=0;
m_heapwidth=0;
m_heapused=0;
m_sumerr=0;
m_sumabs=0;
m_n=0;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CAutoGKInternalState::Copy(CAutoGKInternalState &obj)
{
//--- copy variables
m_a=obj.m_a;
m_b=obj.m_b;
m_eps=obj.m_eps;
m_xwidth=obj.m_xwidth;
m_x=obj.m_x;
m_f=obj.m_f;
m_info=obj.m_info;
m_r=obj.m_r;
m_heapsize=obj.m_heapsize;
m_heapwidth=obj.m_heapwidth;
m_heapused=obj.m_heapused;
m_sumerr=obj.m_sumerr;
m_sumabs=obj.m_sumabs;
m_n=obj.m_n;
m_rstate.Copy(obj.m_rstate);
//--- copy arrays
ArrayCopy(m_qn,obj.m_qn);
ArrayCopy(m_wg,obj.m_wg);
ArrayCopy(m_wk,obj.m_wk);
ArrayCopy(m_wr,obj.m_wr);
//--- copy matrix
m_heap=obj.m_heap;
}
//+------------------------------------------------------------------+
//| This structure stores state of the integration algorithm. |
//| Although this class has public fields, they are not intended for|
//| external use. You should use ALGLIB functions to work with this |
//| class: |
//| * autogksmooth()/AutoGKSmoothW()/... to create objects |
//| * autogkintegrate() to begin integration |
//| * autogkresults() to get results |
//+------------------------------------------------------------------+
class CAutoGKState
{
public:
//--- variables
double m_a;
double m_b;
double m_alpha;
double m_beta;
double m_xwidth;
double m_x;
double m_xminusa;
double m_bminusx;
bool m_needf;
double m_f;
int m_wrappermode;
double m_v;
int m_terminationtype;
int m_nfev;
int m_nintervals;
CAutoGKInternalState m_internalstate;
RCommState m_rstate;
//--- constructor, destructor
CAutoGKState(void);
~CAutoGKState(void) {}
//--- copy
void Copy(CAutoGKState &obj);
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CAutoGKState::CAutoGKState(void)
{
m_a=0;
m_b=0;
m_alpha=0;
m_beta=0;
m_xwidth=0;
m_x=0;
m_xminusa=0;
m_bminusx=0;
m_needf=false;
m_f=0;
m_wrappermode=0;
m_v=0;
m_terminationtype=0;
m_nfev=0;
m_nintervals=0;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CAutoGKState::Copy(CAutoGKState &obj)
{
//--- copy variables
m_a=obj.m_a;
m_b=obj.m_b;
m_alpha=obj.m_alpha;
m_beta=obj.m_beta;
m_xwidth=obj.m_xwidth;
m_x=obj.m_x;
m_xminusa=obj.m_xminusa;
m_bminusx=obj.m_bminusx;
m_needf=obj.m_needf;
m_f=obj.m_f;
m_wrappermode=obj.m_wrappermode;
m_v=obj.m_v;
m_terminationtype=obj.m_terminationtype;
m_nfev=obj.m_nfev;
m_nintervals=obj.m_nintervals;
m_internalstate.Copy(obj.m_internalstate);
m_rstate.Copy(obj.m_rstate);
}
//+------------------------------------------------------------------+
//| This structure stores state of the integration algorithm. |
//| Although this class has public fields, they are not intended for|
//| external use. You should use ALGLIB functions to work with this |
//| class: |
//| * autogksmooth()/AutoGKSmoothW()/... to create objects |
//| * autogkintegrate() to begin integration |
//| * autogkresults() to get results |
//+------------------------------------------------------------------+
class CAutoGKStateShell
{
private:
CAutoGKState m_innerobj;
public:
//--- constructors, destructor
CAutoGKStateShell(void) {}
CAutoGKStateShell(CAutoGKState &obj) { m_innerobj.Copy(obj); }
~CAutoGKStateShell(void) {}
//--- methods
bool GetNeedF(void);
void SetNeedF(const bool b);
double GetX(void);
void SetX(const double d);
double GetXMinusA(void);
void SetXMinusA(const double d);
double GetBMinusX(void);
void SetBMinusX(const double d);
double GetF(void);
void SetF(const double d);
CAutoGKState *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable needf |
//+------------------------------------------------------------------+
bool CAutoGKStateShell::GetNeedF(void)
{
return(m_innerobj.m_needf);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable needf |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetNeedF(const bool b)
{
m_innerobj.m_needf=b;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable x |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetX(void)
{
return(m_innerobj.m_x);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable x |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetX(const double d)
{
m_innerobj.m_x=d;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable xminusa |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetXMinusA(void)
{
return(m_innerobj.m_xminusa);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable xminusa |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetXMinusA(const double d)
{
m_innerobj.m_xminusa=d;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable bminusx |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetBMinusX(void)
{
return(m_innerobj.m_bminusx);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable bminusx |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetBMinusX(const double d)
{
m_innerobj.m_bminusx=d;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable f |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetF(void)
{
return(m_innerobj.m_f);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable f |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetF(const double d)
{
m_innerobj.m_f=d;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CAutoGKState *CAutoGKStateShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Auto Gauss-Kronrod |
//+------------------------------------------------------------------+
class CAutoGK
{
public:
//--- class constant
static const int m_maxsubintervals;
//--- public methods
static void AutoGKSmooth(const double a,const double b,CAutoGKState &state);
static void AutoGKSmoothW(const double a,const double b,const double xwidth,CAutoGKState &state);
static void AutoGKSingular(const double a,const double b,const double alpha,const double beta,CAutoGKState &state);
static void AutoGKResults(CAutoGKState &state,double &v,CAutoGKReport &rep);
static bool AutoGKIteration(CAutoGKState &state);
private:
//--- private methods
static void AutoGKInternalPrepare(const double a,const double b,const double eps,const double xwidth,CAutoGKInternalState &state);
static void MHeapPop(CMatrixDouble &heap,const int heapsize,const int heapwidth);
static void MHeapPush(CMatrixDouble &heap,const int heapsize,const int heapwidth);
static void MHeapResize(CMatrixDouble &heap,int &heapsize,const int newheapsize,const int heapwidth);
static bool AutoGKInternalIteration(CAutoGKInternalState &state);
//--- auxiliary functions for AutoGKInternalIteration
static void Func_Internal_lbl_rcomm(CAutoGKInternalState &state,int i,int j,int ns,int info,double c1,double c2,double intg,double intk,double inta,double v,double ta,double tb,double qeps);
static bool Func_Internal_lbl_5(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_7(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_8(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_11(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_14(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_16(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
static bool Func_Internal_lbl_19(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
//--- auxiliary functions for AutoGKIteration
static void Func_lbl_rcomm(CAutoGKState &state,double s,double tmp,double eps,double a,double b,double x,double t,double alpha,double beta,double v1,double v2);
static bool Func_lbl_5(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
static bool Func_lbl_9(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
static bool Func_lbl_11(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
};
//+------------------------------------------------------------------+
//| Initialize constant |
//+------------------------------------------------------------------+
const int CAutoGK::m_maxsubintervals=10000;
//+------------------------------------------------------------------+
//| Integration of a smooth function F(x) on a finite interval [a,b].|
//| Fast-convergent algorithm based on a Gauss-Kronrod formula is |
//| used. Result is calculated with accuracy close to the machine |
//| precision. |
//| Algorithm works well only with smooth integrands. It may be used |
//| with continuous non-smooth integrands, but with less performance.|
//| It should never be used with integrands which have integrable |
//| singularities at lower or upper limits - algorithm may crash. |
//| Use AutoGKSingular in such cases. |
//| INPUT PARAMETERS: |
//| A, B - interval boundaries (A<B, A=B or A>B) |
//| OUTPUT PARAMETERS |
//| State - structure which stores algorithm state |
//| SEE ALSO |
//| AutoGKSmoothW, AutoGKSingular, AutoGKResults. |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSmooth(const double a,const double b,
CAutoGKState &state)
{
//--- check
if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
return;
//--- function call
AutoGKSmoothW(a,b,0.0,state);
}
//+------------------------------------------------------------------+
//| Integration of a smooth function F(x) on a finite interval [a,b].|
//| This subroutine is same as AutoGKSmooth(), but it guarantees that|
//| interval [a,b] is partitioned into subintervals which have width |
//| at most XWidth. |
//| Subroutine can be used when integrating nearly-constant function |
//| with narrow "bumps" (about XWidth wide). If "bumps" are too |
//| narrow, AutoGKSmooth subroutine can overlook them. |
//| INPUT PARAMETERS: |
//| A, B - interval boundaries (A<B, A=B or A>B) |
//| OUTPUT PARAMETERS |
//| State - structure which stores algorithm state |
//| SEE ALSO |
//| AutoGKSmooth, AutoGKSingular, AutoGKResults. |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSmoothW(const double a,const double b,
const double xwidth,CAutoGKState &state)
{
//--- check
if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(xwidth),__FUNCTION__+": XWidth is not finite!"))
return;
//--- initialization
state.m_wrappermode=0;
state.m_a=a;
state.m_b=b;
state.m_xwidth=xwidth;
state.m_needf=false;
//--- allocation
state.m_rstate.ra.Resize(11);
state.m_rstate.stage=-1;
}
//+------------------------------------------------------------------+
//| Integration on a finite interval [A,B]. |
//| Integrand have integrable singularities at A/B. |
//| F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, |
//| with known alpha/beta (alpha>-1, beta>-1). If alpha/beta are not |
//| known, estimates from below can be used (but these estimates |
//| should be greater than -1 too). |
//| One of alpha/beta variables (or even both alpha/beta) may be |
//| equal to 0, which means than function F(x) is non-singular at |
//| A/B. Anyway (singular at bounds or not), function F(x) is |
//| supposed to be continuous on (A,B). |
//| Fast-convergent algorithm based on a Gauss-Kronrod formula is |
//| used. Result is calculated with accuracy close to the machine |
//| precision. |
//| INPUT PARAMETERS: |
//| A, B - interval boundaries (A<B, A=B or A>B) |
//| Alpha - power-law coefficient of the F(x) at A, |
//| Alpha>-1 |
//| Beta - power-law coefficient of the F(x) at B, |
//| Beta>-1 |
//| OUTPUT PARAMETERS |
//| State - structure which stores algorithm state |
//| SEE ALSO |
//| AutoGKSmooth, AutoGKSmoothW, AutoGKResults. |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSingular(const double a,const double b,
const double alpha,const double beta,
CAutoGKState &state)
{
//--- check
if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(alpha),__FUNCTION__+": Alpha is not finite!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(beta),__FUNCTION__+": Beta is not finite!"))
return;
//--- initialization
state.m_wrappermode=1;
state.m_a=a;
state.m_b=b;
state.m_alpha=alpha;
state.m_beta=beta;
state.m_xwidth=0.0;
state.m_needf=false;
//--- allocation
state.m_rstate.ra.Resize(11);
state.m_rstate.stage=-1;
}
//+------------------------------------------------------------------+
//| Adaptive integration results |
//| Called after AutoGKIteration returned False. |
//| Input parameters: |
//| State - algorithm state (used by AutoGKIteration). |
//| Output parameters: |
//| V - integral(f(x)dx,a,b) |
//| Rep - optimization report (see AutoGKReport |
//| description) |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKResults(CAutoGKState &state,double &v,CAutoGKReport &rep)
{
//--- change values
v=state.m_v;
rep.m_terminationtype=state.m_terminationtype;
rep.m_nfev=state.m_nfev;
rep.m_nintervals=state.m_nintervals;
}
//+------------------------------------------------------------------+
//| Internal AutoGK subroutine |
//| eps<0 - error |
//| eps=0 - automatic eps selection |
//| width<0 - error |
//| width=0 - no width requirements |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKInternalPrepare(const double a,const double b,
const double eps,const double xwidth,
CAutoGKInternalState &state)
{
//--- Save settings
state.m_a=a;
state.m_b=b;
state.m_eps=eps;
state.m_xwidth=xwidth;
//--- Prepare RComm structure
state.m_rstate.ia.Resize(4);
state.m_rstate.ra.Resize(9);
state.m_rstate.stage=-1;
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CAutoGK::MHeapPop(CMatrixDouble &heap,const int heapsize,
const int heapwidth)
{
//--- create variables
int i=0;
int p=0;
double t=0;
int maxcp=0;
//--- check
if(heapsize==1)
return;
//--- initialization
for(i=0; i<=heapwidth-1; i++)
{
t=heap[heapsize-1][i];
heap.Set(heapsize-1,i,heap[0][i]);
heap.Set(0,i,t);
}
p=0;
//--- cycle
while(2*p+1<heapsize-1)
{
maxcp=2*p+1;
//--- check
if(2*p+2<heapsize-1)
{
//--- check
if(heap[2*p+2][0]>heap[2*p+1][0])
maxcp=2*p+2;
}
//--- check
if(heap[p][0]<heap[maxcp][0])
{
for(i=0; i<=heapwidth-1; i++)
{
//--- change values
t=heap[p][i];
heap.Set(p,i,heap[maxcp][i]);
heap.Set(maxcp,i,t);
}
p=maxcp;
}
else
break;
}
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CAutoGK::MHeapPush(CMatrixDouble &heap,const int heapsize,
const int heapwidth)
{
//--- create variables
int i=0;
int p=0;
double t=0;
int parent=0;
//--- check
if(heapsize==0)
return;
p=heapsize;
//--- cycle
while(p!=0)
{
parent=(p-1)/2;
//--- check
if(heap[p][0]>heap[parent][0])
{
for(i=0; i<=heapwidth-1; i++)
{
//--- change values
t=heap[p][i];
heap.Set(p,i,heap[parent][i]);
heap.Set(parent,i,t);
}
p=parent;
}
else
break;
}
}
//+------------------------------------------------------------------+
//| Internal subroutine |
//+------------------------------------------------------------------+
void CAutoGK::MHeapResize(CMatrixDouble &heap,int &heapsize,
const int newheapsize,const int heapwidth)
{
//--- create variables
int i=0;
int i_=0;
//--- create matrix
CMatrixDouble tmp;
//--- allocation
tmp.Resize(heapsize,heapwidth);
//--- copy
for(i=0; i<=heapsize-1; i++)
{
for(i_=0; i_<=heapwidth-1; i_++)
tmp.Set(i,i_,heap[i][i_]);
}
//--- allocation
heap.Resize(newheapsize,heapwidth);
//--- copy
for(i=0; i<=heapsize-1; i++)
{
for(i_=0; i_<=heapwidth-1; i_++)
heap.Set(i,i_,tmp[i][i_]);
}
//--- change value
heapsize=newheapsize;
}
//+------------------------------------------------------------------+
//| Internal AutoGK subroutine |
//+------------------------------------------------------------------+
bool CAutoGK::AutoGKInternalIteration(CAutoGKInternalState &state)
{
//--- create variables
double c1=0;
double c2=0;
int i=0;
int j=0;
double intg=0;
double intk=0;
double inta=0;
double v=0;
double ta=0;
double tb=0;
int ns=0;
double qeps=0;
int info=0;
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
//--- initialization
i=state.m_rstate.ia[0];
j=state.m_rstate.ia[1];
ns=state.m_rstate.ia[2];
info=state.m_rstate.ia[3];
c1=state.m_rstate.ra[0];
c2=state.m_rstate.ra[1];
intg=state.m_rstate.ra[2];
intk=state.m_rstate.ra[3];
inta=state.m_rstate.ra[4];
v=state.m_rstate.ra[5];
ta=state.m_rstate.ra[6];
tb=state.m_rstate.ra[7];
qeps=state.m_rstate.ra[8];
}
else
{
//--- initialization
i=497;
j=-271;
ns=-581;
info=745;
c1=-533;
c2=-77;
intg=678;
intk=-293;
inta=316;
v=647;
ta=-756;
tb=830;
qeps=-871;
}
//--- check
if(state.m_rstate.stage==0)
{
//--- change value
v=state.m_f;
//--- Gauss-Kronrod formula
intk=intk+v*state.m_wk[i];
//--- check
if(i%2==1)
intg=intg+v*state.m_wg[i];
//--- Integral |F(x)|
//--- Use rectangles method
inta=inta+MathAbs(v)*state.m_wr[i];
i=i+1;
//--- function call, return result
return(Func_Internal_lbl_5(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- check
if(state.m_rstate.stage==1)
{
//--- change value
v=state.m_f;
//--- Gauss-Kronrod formula
intk=intk+v*state.m_wk[i];
//--- check
if(i%2==1)
intg=intg+v*state.m_wg[i];
//--- Integral |F(x)|
//--- Use rectangles method
inta=inta+MathAbs(v)*state.m_wr[i];
i=i+1;
//--- function call, return result
return(Func_Internal_lbl_11(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- check
if(state.m_rstate.stage==2)
{
//--- change value
v=state.m_f;
//--- Gauss-Kronrod formula
intk=intk+v*state.m_wk[i];
//--- check
if(i%2==1)
intg=intg+v*state.m_wg[i];
//--- Integral |F(x)|
//--- Use rectangles method
inta=inta+MathAbs(v)*state.m_wr[i];
i=i+1;
//--- function call, return result
return(Func_Internal_lbl_19(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- Routine body
//--- initialize quadratures.
//--- use 15-point Gauss-Kronrod formula.
state.m_n=15;
CGaussKronrodQ::GKQGenerateGaussLegendre(state.m_n,info,state.m_qn,state.m_wk,state.m_wg);
//--- check
if(info<0)
{
//--- change values
state.m_info=-5;
state.m_r=0;
//--- return result
return(false);
}
//--- allocation
ArrayResize(state.m_wr,state.m_n);
for(i=0; i<=state.m_n-1; i++)
{
//--- check
if(i==0)
{
state.m_wr[i]=0.5*MathAbs(state.m_qn[1]-state.m_qn[0]);
continue;
}
//--- check
if(i==state.m_n-1)
{
state.m_wr[state.m_n-1]=0.5*MathAbs(state.m_qn[state.m_n-1]-state.m_qn[state.m_n-2]);
continue;
}
state.m_wr[i]=0.5*MathAbs(state.m_qn[i-1]-state.m_qn[i+1]);
}
//--- special case
if(state.m_a==state.m_b)
{
//--- change values
state.m_info=1;
state.m_r=0;
//--- return result
return(false);
}
//--- test parameters
if(state.m_eps<0.0 || state.m_xwidth<0.0)
{
//--- change values
state.m_info=-1;
state.m_r=0;
//--- return result
return(false);
}
state.m_info=1;
//--- check
if(state.m_eps==0.0)
state.m_eps=100000*CMath::m_machineepsilon;
//--- First,prepare heap
//--- * column 0 - absolute error
//--- * column 1 - integral of a F(x) (calculated using Kronrod extension nodes)
//--- * column 2 - integral of a |F(x)| (calculated using modified rect. method)
//--- * column 3 - left boundary of a subinterval
//--- * column 4 - right boundary of a subinterval
if(state.m_xwidth!=0.0)
{
//--- maximum subinterval should be no more than XWidth.
//--- so we create Ceil((B-A)/XWidth)+1 small subintervals
ns=(int)MathCeil(MathAbs(state.m_b-state.m_a)/state.m_xwidth)+1;
state.m_heapsize=ns;
state.m_heapused=ns;
state.m_heapwidth=5;
state.m_heap.Resize(state.m_heapsize,state.m_heapwidth);
state.m_sumerr=0;
state.m_sumabs=0;
j=0;
//--- function call, return result
return(Func_Internal_lbl_8(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- no maximum width requirements
//--- start from one big subinterval
state.m_heapwidth=5;
state.m_heapsize=1;
state.m_heapused=1;
state.m_heap.Resize(state.m_heapsize,state.m_heapwidth);
c1=0.5*(state.m_b-state.m_a);
c2=0.5*(state.m_b+state.m_a);
intg=0;
intk=0;
inta=0;
i=0;
//--- function call, return result
return(Func_Internal_lbl_5(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
void CAutoGK::Func_Internal_lbl_rcomm(CAutoGKInternalState &state,
int i,int j,int ns,int info,
double c1,double c2,double intg,
double intk,double inta,double v,
double ta,double tb,double qeps)
{
//--- save
state.m_rstate.ia.Set(0,i);
state.m_rstate.ia.Set(1,j);
state.m_rstate.ia.Set(2,ns);
state.m_rstate.ia.Set(3,info);
state.m_rstate.ra.Set(0,c1);
state.m_rstate.ra.Set(1,c2);
state.m_rstate.ra.Set(2,intg);
state.m_rstate.ra.Set(3,intk);
state.m_rstate.ra.Set(4,inta);
state.m_rstate.ra.Set(5,v);
state.m_rstate.ra.Set(6,ta);
state.m_rstate.ra.Set(7,tb);
state.m_rstate.ra.Set(8,qeps);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_5(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- check
if(i>state.m_n-1)
return(Func_Internal_lbl_7(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
//--- obtain F
state.m_x=c1*state.m_qn[i]+c2;
state.m_rstate.stage=0;
//--- Saving state
Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_7(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- calculation
intk=intk*(state.m_b-state.m_a)*0.5;
intg=intg*(state.m_b-state.m_a)*0.5;
inta=inta*(state.m_b-state.m_a)*0.5;
state.m_heap.Set(0,0,MathAbs(intg-intk));
state.m_heap.Set(0,1,intk);
state.m_heap.Set(0,2,inta);
state.m_heap.Set(0,3,state.m_a);
state.m_heap.Set(0,4,state.m_b);
state.m_sumerr=state.m_heap[0][0];
state.m_sumabs=MathAbs(inta);
//--- method iterations
return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_8(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- check
if(j>ns-1)
{
//--- method iterations
return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- calculation
ta=state.m_a+j*(state.m_b-state.m_a)/ns;
tb=state.m_a+(j+1)*(state.m_b-state.m_a)/ns;
c1=0.5*(tb-ta);
c2=0.5*(tb+ta);
intg=0;
intk=0;
inta=0;
i=0;
//--- function call, return result
return(Func_Internal_lbl_11(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_11(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- check
if(i>state.m_n-1)
{
//--- calculation
intk=intk*(tb-ta)*0.5;
intg=intg*(tb-ta)*0.5;
inta=inta*(tb-ta)*0.5;
state.m_heap.Set(j,0,MathAbs(intg-intk));
state.m_heap.Set(j,1,intk);
state.m_heap.Set(j,2,inta);
state.m_heap.Set(j,3,ta);
state.m_heap.Set(j,4,tb);
state.m_sumerr=state.m_sumerr+state.m_heap[j][0];
state.m_sumabs=state.m_sumabs+MathAbs(inta);
j=j+1;
//--- function call, return result
return(Func_Internal_lbl_8(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- obtain F
state.m_x=c1*state.m_qn[i]+c2;
state.m_rstate.stage=1;
//--- Saving state
Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_14(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- additional memory if needed
if(state.m_heapused==state.m_heapsize)
MHeapResize(state.m_heap,state.m_heapsize,4*state.m_heapsize,state.m_heapwidth);
//--- TODO: every 20 iterations recalculate errors/sums
if(state.m_sumerr<=state.m_eps*state.m_sumabs || state.m_heapused>=m_maxsubintervals)
{
state.m_r=0;
for(j=0; j<=state.m_heapused-1; j++)
state.m_r=state.m_r+state.m_heap[j][1];
//--- return result
return(false);
}
//--- Exclude interval with maximum absolute error
MHeapPop(state.m_heap,state.m_heapused,state.m_heapwidth);
state.m_sumerr=state.m_sumerr-state.m_heap[state.m_heapused-1][0];
state.m_sumabs=state.m_sumabs-state.m_heap[state.m_heapused-1][2];
//--- Divide interval,create subintervals
ta=state.m_heap[state.m_heapused-1][3];
tb=state.m_heap[state.m_heapused-1][4];
state.m_heap.Set(state.m_heapused-1,3,ta);
state.m_heap.Set(state.m_heapused-1,4,0.5*(ta+tb));
state.m_heap.Set(state.m_heapused,3,0.5*(ta+tb));
state.m_heap.Set(state.m_heapused,4,tb);
j=state.m_heapused-1;
//--- function call, return result
return(Func_Internal_lbl_16(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_16(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- check
if(j>state.m_heapused)
{
//--- function call
MHeapPush(state.m_heap,state.m_heapused-1,state.m_heapwidth);
//--- function call
MHeapPush(state.m_heap,state.m_heapused,state.m_heapwidth);
state.m_heapused=state.m_heapused+1;
//--- method iterations
return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- calculation
c1=0.5*(state.m_heap[j][4]-state.m_heap[j][3]);
c2=0.5*(state.m_heap[j][4]+state.m_heap[j][3]);
intg=0;
intk=0;
inta=0;
i=0;
//--- function call, return result
return(Func_Internal_lbl_19(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to |
//| get rid of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_19(CAutoGKInternalState &state,
int &i,int &j,int &ns,int &info,
double &c1,double &c2,double &intg,
double &intk,double &inta,double &v,
double &ta,double &tb,double &qeps)
{
//--- check
if(i>state.m_n-1)
{
//--- calculation
intk=intk*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
intg=intg*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
inta=inta*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
state.m_heap.Set(j,0,MathAbs(intg-intk));
state.m_heap.Set(j,1,intk);
state.m_heap.Set(j,2,inta);
state.m_sumerr=state.m_sumerr+state.m_heap[j][0];
state.m_sumabs=state.m_sumabs+state.m_heap[j][2];
j=j+1;
//--- function call, return result
return(Func_Internal_lbl_16(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
}
//--- F(x)
state.m_x=c1*state.m_qn[i]+c2;
state.m_rstate.stage=2;
//--- Saving state
Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Iterative method |
//+------------------------------------------------------------------+
bool CAutoGK::AutoGKIteration(CAutoGKState &state)
{
//--- create variables
double s=0;
double tmp=0;
double eps=0;
double a=0;
double b=0;
double x=0;
double t=0;
double alpha=0;
double beta=0;
double v1=0;
double v2=0;
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
//--- initialization
s=state.m_rstate.ra[0];
tmp=state.m_rstate.ra[1];
eps=state.m_rstate.ra[2];
a=state.m_rstate.ra[3];
b=state.m_rstate.ra[4];
x=state.m_rstate.ra[5];
t=state.m_rstate.ra[6];
alpha=state.m_rstate.ra[7];
beta=state.m_rstate.ra[8];
v1=state.m_rstate.ra[9];
v2=state.m_rstate.ra[10];
}
else
{
//--- initialization
s=-983;
tmp=-989;
eps=-834;
a=900;
b=-287;
x=364;
t=214;
alpha=-338;
beta=-686;
v1=912;
v2=585;
}
//--- check
if(state.m_rstate.stage==0)
{
//--- change values
state.m_needf=false;
state.m_nfev=state.m_nfev+1;
state.m_internalstate.m_f=state.m_f;
//--- function call, return result
return(Func_lbl_5(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//--- check
if(state.m_rstate.stage==1)
{
//--- change value
state.m_needf=false;
//--- check
if(alpha!=0.0)
state.m_internalstate.m_f=state.m_f*MathPow(x,-(alpha/(1+alpha)))/(1+alpha);
else
state.m_internalstate.m_f=state.m_f;
state.m_nfev=state.m_nfev+1;
//--- function call, return result
return(Func_lbl_9(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//--- check
if(state.m_rstate.stage==2)
{
//--- change value
state.m_needf=false;
//--- check
if(beta!=0.0)
state.m_internalstate.m_f=state.m_f*MathPow(x,-(beta/(1+beta)))/(1+beta);
else
state.m_internalstate.m_f=state.m_f;
state.m_nfev=state.m_nfev+1;
//--- function call, return result
return(Func_lbl_11(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//--- Routine body
eps=0;
a=state.m_a;
b=state.m_b;
alpha=state.m_alpha;
beta=state.m_beta;
state.m_terminationtype=-1;
state.m_nfev=0;
state.m_nintervals=0;
//--- smooth function at a finite interval
if(state.m_wrappermode!=0)
{
//--- function with power-law singularities at the ends of a finite interval
if(state.m_wrappermode!=1)
return(false);
//--- test coefficients
if(alpha<=-1.0 || beta<=-1.0)
{
//--- change values
state.m_terminationtype=-1;
state.m_v=0;
//--- return result
return(false);
}
//--- special cases
if(a==b)
{
//--- change values
state.m_terminationtype=1;
state.m_v=0;
//--- return result
return(false);
}
//--- reduction to general form
if(a<b)
s=1;
else
{
//--- change values
s=-1;
tmp=a;
a=b;
b=tmp;
tmp=alpha;
alpha=beta;
beta=tmp;
}
//--- change values
alpha=MathMin(alpha,0);
beta=MathMin(beta,0);
//--- first,integrate left half of [a,b]:
//--- integral(f(x)dx,a,(b+a)/2)=
//--- =1/(1+alpha) * integral(t^(-alpha/(1+alpha))*f(a+t^(1/(1+alpha)))dt,0,(0.5*(b-a))^(1+alpha))
AutoGKInternalPrepare(0,MathPow(0.5*(b-a),1+alpha),eps,state.m_xwidth,state.m_internalstate);
//--- function call, return result
return(Func_lbl_9(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//--- special case
if(a==b)
{
//--- change values
state.m_terminationtype=1;
state.m_v=0;
//--- return result
return(false);
}
//--- general case
AutoGKInternalPrepare(a,b,eps,state.m_xwidth,state.m_internalstate);
//--- function call, return result
return(Func_lbl_5(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid |
//| of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
void CAutoGK::Func_lbl_rcomm(CAutoGKState &state,double s,double tmp,
double eps,double a,double b,double x,
double t,double alpha,double beta,
double v1,double v2)
{
//--- save
state.m_rstate.ra.Set(0,s);
state.m_rstate.ra.Set(1,tmp);
state.m_rstate.ra.Set(2,eps);
state.m_rstate.ra.Set(3,a);
state.m_rstate.ra.Set(4,b);
state.m_rstate.ra.Set(5,x);
state.m_rstate.ra.Set(6,t);
state.m_rstate.ra.Set(7,alpha);
state.m_rstate.ra.Set(8,beta);
state.m_rstate.ra.Set(9,v1);
state.m_rstate.ra.Set(10,v2);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid |
//| of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_5(CAutoGKState &state,double &s,double &tmp,
double &eps,double &a,double &b,double &x,
double &t,double &alpha,double &beta,
double &v1,double &v2)
{
//--- check
if(!AutoGKInternalIteration(state.m_internalstate))
{
//--- change values
state.m_v=state.m_internalstate.m_r;
state.m_terminationtype=state.m_internalstate.m_info;
state.m_nintervals=state.m_internalstate.m_heapused;
//--- return result
return(false);
}
//--- change values
x=state.m_internalstate.m_x;
state.m_x=x;
state.m_xminusa=x-a;
state.m_bminusx=b-x;
state.m_needf=true;
state.m_rstate.stage=0;
//--- Saving state
Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid |
//| of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_9(CAutoGKState &state,double &s,double &tmp,
double &eps,double &a,double &b,double &x,
double &t,double &alpha,double &beta,
double &v1,double &v2)
{
//--- check
if(!AutoGKInternalIteration(state.m_internalstate))
{
v1=state.m_internalstate.m_r;
state.m_nintervals=state.m_nintervals+state.m_internalstate.m_heapused;
//--- then,integrate right half of [a,b]:
//--- integral(f(x)dx,(b+a)/2,b)=
//--- =1/(1+beta) * integral(t^(-beta/(1+beta))*f(b-t^(1/(1+beta)))dt,0,(0.5*(b-a))^(1+beta))
AutoGKInternalPrepare(0,MathPow(0.5*(b-a),1+beta),eps,state.m_xwidth,state.m_internalstate);
//--- function call, return result
return(Func_lbl_11(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
}
//--- Fill State.X,State.XMinusA,State.BMinusX.
//--- Latter two are filled correctly even if B<A.
x=state.m_internalstate.m_x;
t=MathPow(x,1/(1+alpha));
state.m_x=a+t;
//--- check
if(s>0.0)
{
state.m_xminusa=t;
state.m_bminusx=b-(a+t);
}
else
{
state.m_xminusa=a+t-b;
state.m_bminusx=-t;
}
//--- change values
state.m_needf=true;
state.m_rstate.stage=1;
//--- Saving state
Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid |
//| of the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_11(CAutoGKState &state,double &s,double &tmp,
double &eps,double &a,double &b,double &x,
double &t,double &alpha,double &beta,
double &v1,double &v2)
{
//--- check
if(!AutoGKInternalIteration(state.m_internalstate))
{
//--- change values
v2=state.m_internalstate.m_r;
state.m_nintervals=state.m_nintervals+state.m_internalstate.m_heapused;
//--- final result
state.m_v=s*(v1+v2);
state.m_terminationtype=1;
//--- return result
return(false);
}
//--- Fill State.X,State.XMinusA,State.BMinusX.
//--- Latter two are filled correctly (X-A,B-X) even if B<A.
x=state.m_internalstate.m_x;
t=MathPow(x,1/(1+beta));
state.m_x=b-t;
//--- check
if(s>0.0)
{
state.m_xminusa=b-t-a;
state.m_bminusx=t;
}
else
{
state.m_xminusa=-t;
state.m_bminusx=a-(b-t);
}
//--- change values
state.m_needf=true;
state.m_rstate.stage=2;
//--- Saving state
Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
Reference in a new issueCopy permalink