| //+------------------------------------------------------------------+
|
| //| 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);
|
| }
|
| //+------------------------------------------------------------------+
|
