//+------------------------------------------------------------------+ //| 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 (AB) | //| 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 (AB) | //| 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 (AB) | //| 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+1heap[2*p+1][0]) maxcp=2*p+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(a0.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 B0.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); } //+------------------------------------------------------------------+