//+------------------------------------------------------------------+ //| ap.mqh | //| Copyright 2003-2012 Sergey Bochkanov (ALGLIB project) | //| Copyright 2012-2017, MetaQuotes Software Corp. | //| 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 #include "complex.mqh" #include "matrix.mqh" #include "bitconvert.mqh" //+------------------------------------------------------------------+ //| Reverse communication structure | //+------------------------------------------------------------------+ struct RCommState { public: int stage; int ia[]; bool ba[]; double ra[]; complex ca[]; RCommState(void) { stage=-1; } ~RCommState(void) { }; void Copy(RCommState &obj); }; //+------------------------------------------------------------------+ //| Create a copy | //+------------------------------------------------------------------+ void RCommState::Copy(RCommState &obj) { //--- copy a variable stage=obj.stage; //--- copy arrays ArrayCopy(ia,obj.ia); ArrayCopy(ba,obj.ba); ArrayCopy(ra,obj.ra); ArrayCopy(ca,obj.ca); } //+------------------------------------------------------------------+ //| Internal functions | //+------------------------------------------------------------------+ class CAp { public: //--- variable that determines whether an exception happened static bool exception_happened; //--- constructor, destructor CAp(void); ~CAp(void); //--- len static int Len(const int &a[]); static int Len(const bool &a[]); static int Len(const double &a[]); static int Len(const complex &a[]); //--- rows count static int Rows(const CMatrixInt &a); static int Rows(const CMatrixDouble &a); static int Rows(const CMatrixComplex &a); //--- cols count static int Cols(const CMatrixInt &a); static int Cols(const CMatrixDouble &a); static int Cols(const CMatrixComplex &a); //--- swap static void Swap(int &a,int &b); static void Swap(double &a,double &b); static void Swap(complex &a,complex &b); static void Swap(bool &a[],bool &b[]); static void Swap(int &a[],int &b[]); static void Swap(double &a[],double &b[]); static void Swap(complex &a[],complex &b[]); static void Swap(CMatrixInt &a,CMatrixInt &b); static void Swap(CMatrixDouble &a,CMatrixDouble &b); static void Swap(CMatrixComplex &a,CMatrixComplex &b); //--- check assertions static bool Assert(const bool cond); static bool Assert(const bool cond,const string s); //--- determination of accuracy static int ThresHoldToDPS(const double threshold); //--- join string static string StringJoin(const string sep,const string &a[]); //--- convert to string static string Format(const complex &a,const int dps); static string Format(const bool &a[]); static string Format(const int &a[]); static string Format(const double &a[],const int dps); static string Format(const complex &a[],const int dps); static string FormatB(const CMatrixInt &a); static string Format(const CMatrixInt &a); static string Format(const CMatrixDouble &a,const int dps); static string Format(const CMatrixComplex &a,const int dps); //--- work with matrix static bool IsSymmetric(const CMatrixDouble &a); static bool IsHermitian(const CMatrixComplex &a); static bool ForceSymmetric(CMatrixDouble &a); static bool ForceHermitian(CMatrixComplex &a); }; //+------------------------------------------------------------------+ //| Initialize variable | //+------------------------------------------------------------------+ bool CAp::exception_happened=false; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CAp::CAp(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CAp::~CAp(void) { } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ static int CAp::Len(const int &a[]) { return(ArraySize(a)); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ static int CAp::Len(const bool &a[]) { return(ArraySize(a)); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ static int CAp::Len(const double &a[]) { return(ArraySize(a)); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ static int CAp::Len(const complex &a[]) { return(ArraySize(a)); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ static int CAp::Rows(const CMatrixInt &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ static int CAp::Rows(const CMatrixDouble &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ static int CAp::Rows(const CMatrixComplex &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get cols count | //+------------------------------------------------------------------+ static int CAp::Cols(const CMatrixInt &a) { //--- check if(a.Size()==0) return(0); //--- return result return(a[0].Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ static int CAp::Cols(const CMatrixDouble &a) { //--- check if(a.Size()==0) return(0); //--- return result return(a[0].Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ static int CAp::Cols(const CMatrixComplex &a) { //--- check if(a.Size()==0) return(0); //--- return result return(a[0].Size()); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(int &a,int &b) { int t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(double &a,double &b) { double t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(complex &a,complex &b) { complex t(a.re,a.im); a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(bool &a[],bool &b[]) { //--- calculation int na=ArraySize(a); int nb=ArraySize(b); //--- create array bool t[]; ArrayResizeAL(t,na); //--- swap ArrayCopy(t,a); ArrayResizeAL(a,nb); ArrayCopy(a,b); ArrayResizeAL(b,na); ArrayCopy(b,t); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(int &a[],int &b[]) { //--- calculation int na=ArraySize(a); int nb=ArraySize(b); //--- create array int t[]; ArrayResizeAL(t,na); //--- swap ArrayCopy(t,a); ArrayResizeAL(a,nb); ArrayCopy(a,b); ArrayResizeAL(b,na); ArrayCopy(b,t); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(double &a[],double &b[]) { //--- calculation int na=ArraySize(a); int nb=ArraySize(b); //--- create array double t[]; ArrayResizeAL(t,na); //--- swap ArrayCopy(t,a); ArrayResizeAL(a,nb); ArrayCopy(a,b); ArrayResizeAL(b,na); ArrayCopy(b,t); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(complex &a[],complex &b[]) { //--- calculation int na=ArraySize(a); int nb=ArraySize(b); //--- create array complex t[]; ArrayResizeAL(t,na); //--- swap ArrayCopy(t,a); ArrayResizeAL(a,nb); ArrayCopy(a,b); ArrayResizeAL(b,na); ArrayCopy(b,t); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(CMatrixInt &a,CMatrixInt &b) { //--- create matrix CMatrixInt t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(CMatrixDouble &a,CMatrixDouble &b) { //--- create matrix CMatrixDouble t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ static void CAp::Swap(CMatrixComplex &a,CMatrixComplex &b) { //--- create matrix CMatrixComplex t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Check assertions | //+------------------------------------------------------------------+ static bool CAp::Assert(const bool cond) { return(Assert(cond,"ALGLIB: assertion failed")); } //+------------------------------------------------------------------+ //| Check assertions | //+------------------------------------------------------------------+ static bool CAp::Assert(const bool cond,const string s) { //--- check if(cond==0) { Print(__FUNCTION__+" "+s); exception_happened=true; return(false); } //--- the assertion is true return(true); } //+------------------------------------------------------------------+ //| returns dps (digits-of-precision) value corresponding to | //| threshold. | //| dps(0.9) = dps(0.5) = dps(0.1) = 0 | //| dps(0.09) = dps(0.05) = dps(0.01) = 1 | //| and so on | //+------------------------------------------------------------------+ static int CAp::ThresHoldToDPS(const double threshold) { //--- initialization int res=0; double t=1.0; for(res=0;t/10>threshold*(1+1E-10);res++) t/=10; //--- return result return(res); } //+------------------------------------------------------------------+ //| Concatenation | //+------------------------------------------------------------------+ static string CAp::StringJoin(const string sep,const string &a[]) { int size=ArraySize(a); //--- check if(size==0) { Print(__FUNCTION__+": array size error"); return(NULL); } //--- concatenation string res=""; for(int i=0;i=0) fmt="f"; else fmt="e"; //--- get sign of the imaginary part string sign; if(a.im>=0) sign="+"; else sign="-"; //--- converting int d=(int)MathAbs(dps); string fmtx=StringFormat(".%d"+fmt,d); string fmty=StringFormat(".%d"+fmt,d); //--- get result string res=StringFormat("%"+fmtx,a.re)+sign+ StringFormat("%"+fmty,MathAbs(a.im))+"i"; StringReplace(res,",","."); //--- return result return(res); } //+------------------------------------------------------------------+ //| Prints formatted array | //+------------------------------------------------------------------+ static string CAp::Format(const bool &a[]) { int size=ArraySize(a); //--- check if(size==0) { Print(__FUNCTION__+": array size error"); return(NULL); } //--- converting string result[]; ArrayResizeAL(result,size); for(int i=0;i=0) sfmt="f"; else sfmt="e"; //--- converting int d=(int)MathAbs(dps); string fmt=StringFormat(".%d"+sfmt,d); for(int i=0;i=0) fmt="f"; else fmt="e"; //--- converting int d=(int)MathAbs(dps); string fmtx=StringFormat(".%d"+fmt,d); string fmty=StringFormat(".%d"+fmt,d); string sign; for(int i=0;i=0) sign="+"; else sign="-"; //--- fill result result[i]=StringFormat("%"+fmtx,a[i].re)+sign+ StringFormat("%"+fmty,MathAbs(a[i].im))+"i"; StringReplace(result[i],",","."); } //--- return result return("{"+StringJoin(",",result)+"}"); } //+------------------------------------------------------------------+ //| Prints formatted matrix | //+------------------------------------------------------------------+ static string CAp::FormatB(const CMatrixInt &a) { int m=a.Size(); //--- check if(m==0) { Print(__FUNCTION__+": array size error"); return(NULL); } int n=a[0].Size(); //--- check if(n==0) { Print(__FUNCTION__+": array size error"); return(NULL); } //--- prepare arrays bool line[]; string result[]; ArrayResizeAL(line,n); ArrayResizeAL(result,m); //--- converting for(int i=0;iK*N) if(!CAp::Assert(n>0,__FUNCTION__+": N<=0!")) return(-1); //--- check if(!CAp::Assert(n=RNDBaseMax-1!")) return(-1); //--- initialization mx=m_HQRndMax-1-(m_HQRndMax-1)%n; do result=HQRndIntegerBase(state)-1; while(result>=mx); //--- get result result=result%n; //--- return result return(result); } //+------------------------------------------------------------------+ //| Random number generator: normal numbers | //| This function generates one random number from normal | //| distribution. | //| Its performance is equal to that of HQRNDNormal2() | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ static double CHighQualityRand::HQRndNormal(CHighQualityRandState &state) { //--- create variables double v1=0; double v2=0; //--- function call HQRndNormal2(state,v1,v2); //--- return result return(v1); } //+------------------------------------------------------------------+ //| Random number generator: random X and Y such that X^2+Y^2=1 | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ static void CHighQualityRand::HQRndUnit2(CHighQualityRandState &state, double &x,double &y) { //--- create variables double v=0; double mx=0; double mn=0; //--- initialization x=0; y=0; //--- function call do HQRndNormal2(state,x,y); while(!(x!=0.0 || y!=0.0)); //--- change values mx=MathMax(MathAbs(x),MathAbs(y)); mn=MathMin(MathAbs(x),MathAbs(y)); v=mx*MathSqrt(1+CMath::Sqr(mn/mx)); //--- get result x=x/v; y=y/v; } //+------------------------------------------------------------------+ //| Random number generator: normal numbers | //| This function generates two independent random numbers from | //| normal distribution. Its performance is equal to that of | //| HQRNDNormal() | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ static void CHighQualityRand::HQRndNormal2(CHighQualityRandState &state, double &x1,double &x2) { //--- create variables double u=0; double v=0; double s=0; //--- initialization x1=0; x2=0; //--- cycle while(true) { u=2*HQRndUniformR(state)-1; v=2*HQRndUniformR(state)-1; s=CMath::Sqr(u)+CMath::Sqr(v); //--- check if(s>0.0 && s<1.0) { //--- two Sqrt's instead of one to //--- avoid overflow when S is too small s=MathSqrt(-(2*MathLog(s)))/MathSqrt(s); x1=u*s; x2=v*s; //--- exit the function return; } } } //+------------------------------------------------------------------+ //| Random number generator: exponential distribution | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ static double CHighQualityRand::HQRndExponential(CHighQualityRandState &state, const double lambdav) { //--- check if(!CAp::Assert(lambdav>0.0,__FUNCTION__+": LambdaV<=0!")) return(EMPTY_VALUE); //--- return result return(-(MathLog(HQRndUniformR(state))/lambdav)); } //+------------------------------------------------------------------+ //| L'Ecuyer, Efficient and portable combined random number | //| generators | //+------------------------------------------------------------------+ static int CHighQualityRand::HQRndIntegerBase(CHighQualityRandState &state) { //--- create variables int result=0; int k=0; //--- check if(!CAp::Assert(state.m_magicv==m_HQRndMagic,__FUNCTION__+": State is not correctly initialized!")) return(-1); //--- initialization k=state.m_s1/53668; state.m_s1=40014*(state.m_s1-k*53668)-k*12211; //--- check if(state.m_s1<0) state.m_s1=state.m_s1+2147483563; //--- change values k=state.m_s2/52774; state.m_s2=40692*(state.m_s2-k*52774)-k*3791; //--- check if(state.m_s2<0) state.m_s2=state.m_s2+2147483399; //--- Result result=state.m_s1-state.m_s2; //--- check if(result<1) result=result+2147483562; //--- return result return(result); } //+------------------------------------------------------------------+ //| Math functions | //+------------------------------------------------------------------+ class CMath { public: //--- class variables static bool m_first_call; static double m_last; static CHighQualityRandState m_state; //--- machine constants static const double m_machineepsilon; static const double m_maxrealnumber; static const double m_minrealnumber; //--- constructor, destructor CMath(void); ~CMath(void); //--- methods static bool IsFinite(const double d); static double RandomReal(void); static int RandomInteger(const int n); static double Sqr(const double x) { return(x*x); } static double AbsComplex(const complex &z); static double AbsComplex(const double r); static complex Conj(const complex &z); static complex Csqr(const complex &z); }; //+------------------------------------------------------------------+ //| Initialize class constants | //+------------------------------------------------------------------+ const double CMath::m_machineepsilon=5E-16; const double CMath::m_maxrealnumber=1E300; const double CMath::m_minrealnumber=1E-300; bool CMath::m_first_call=true; double CMath::m_last=0.0; CHighQualityRandState CMath::m_state; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CMath::CMath(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CMath::~CMath(void) { } //+------------------------------------------------------------------+ //| Check on +-inf | //+------------------------------------------------------------------+ static bool CMath::IsFinite(const double d) { //--- return result return(MathIsValidNumber(d)); } //+------------------------------------------------------------------+ //| Random real value [0,1) | //+------------------------------------------------------------------+ static double CMath::RandomReal(void) { //--- create variable double result; //--- check if(m_first_call) { CHighQualityRand::HQRndSeed(1+MathRand(),1+MathRand(),m_state); m_first_call=false; } //--- get value result=CHighQualityRand::HQRndUniformR(m_state); //--- check if(result==m_last) { m_first_call=true; return(RandomReal()); } //--- change value m_last=result; //--- return result return(CHighQualityRand::HQRndUniformR(m_state)); } //+------------------------------------------------------------------+ //| Random integer value | //+------------------------------------------------------------------+ static int CMath::RandomInteger(const int n) { //--- check if(m_first_call) { CHighQualityRand::HQRndSeed(1+MathRand(),1+MathRand(),m_state); m_first_call=false; } //--- check and return result if(n>=CHighQualityRand::m_HQRndM1-1) return(CHighQualityRand::HQRndUniformI(m_state,CHighQualityRand::m_HQRndM1-2)); else return(CHighQualityRand::HQRndUniformI(m_state,n)); } //+------------------------------------------------------------------+ //| The absolute value of a complex number | //+------------------------------------------------------------------+ static double CMath::AbsComplex(const complex &z) { //--- initialization double w=0.0; double v=0.0; double xabs=MathAbs(z.re); double yabs=MathAbs(z.im); //--- check if(xabs>yabs) w=xabs; else w=yabs; //--- check if(xabsyabs) w=xabs; else w=yabs; //--- check if(xabs63) return('?'); //--- return result return(_sixbits2char_tbl[v]); } //+------------------------------------------------------------------+ //| This function converts character to six-bit value (from 0 to 63).| //| This function is inverse of ae_sixbits2char() | //| If c is not correct character, this function returns -1. | //+------------------------------------------------------------------+ static int CSerializer::Char2SixBits(const char c) { //--- check if(c>=0 && c<127) return(_char2sixbits_tbl[c]); //--- return result return(-1); } //+------------------------------------------------------------------+ //| This function converts three bytes (24 bits) to four six-bit | //| values (24 bits again). | //| src array | //| src_offs offset of three-bytes chunk | //| dst array for ints | //| dst_offs offset of four-ints chunk | //+------------------------------------------------------------------+ static void CSerializer::ThreeBytes2FourSixBits(uchar &src[],const int src_offs, int &dst[],const int dst_offs) { //--- get bits dst[dst_offs+0]=src[src_offs+0] & 0x3F; dst[dst_offs+1]=(src[src_offs+0]>>6) | ((src[src_offs+1]&0x0F)<<2); dst[dst_offs+2]=(src[src_offs+1]>>4) | ((src[src_offs+2]&0x03)<<4); dst[dst_offs+3]=src[src_offs+2]>>2; } //+------------------------------------------------------------------+ //| This function converts four six-bit values (24 bits) to three | //| bytes (24 bits again). | //| src pointer to four ints | //| src_offs offset of the chunk | //| dst pointer to three bytes | //| dst_offs offset of the chunk | //+------------------------------------------------------------------+ static void CSerializer::FourSixBits2ThreeBytes(int &src[],const int src_offs, uchar &dst[],const int dst_offs) { //--- get bytes dst[dst_offs+0]=(uchar)(src[src_offs+0] | ((src[src_offs+1]&0x03)<<6)); dst[dst_offs+1]=(uchar)((src[src_offs+1]>>2) | ((src[src_offs+2]&0x0F)<<4)); dst[dst_offs+2]=(uchar)((src[src_offs+2]>>4) | (src[src_offs+3]<<2)); } //+------------------------------------------------------------------+ //| This function serializes boolean value into buffer | //| v boolean value to be serialized | //| buf buffer, at least 11 characters wide | //| offs offset in the buffer | //| after return(from this function, offs points to the char's past | //| the value being read. | //+------------------------------------------------------------------+ static void CSerializer::Bool2Str(const bool v,char &buf[],int &offs) { //--- create variables char c; int i; //--- check if(v) c='1'; else c='0'; //--- copy c for(i=0;i=m_ser_entry_length) { Print(__FUNCTION__+" "+emsg); //--- return result return(-1); } sixbits[sixbitsread]=d; sixbitsread++; offs++; } //--- check if(sixbitsread==0) { Print(__FUNCTION__+" "+emsg); //--- return result return(-1); } for(i=sixbitsread;i<12;i++) sixbits[i]=0; //--- function call FourSixBits2ThreeBytes(sixbits,0,bytes,0); //--- function call FourSixBits2ThreeBytes(sixbits,4,bytes,3); //--- function call FourSixBits2ThreeBytes(sixbits,8,bytes,6); //--- check if((bytes[sizeof(int)-1]&0x80)!=0) c=(uchar)0xFF; else c=(uchar)0x00; for(i=sizeof(int);i<8;i++) //--- check if(bytes[i]!=c) { Print(__FUNCTION__+" "+emsg3264); //--- return result return(-1); } //--- copy for(i=0;i=m_ser_entry_length) { Print(__FUNCTION__+"emsg"); //--- return result return(EMPTY_VALUE); } sixbits[sixbitsread]=d; sixbitsread++; offs++; } //--- check if(sixbitsread!=m_ser_entry_length) { Print(__FUNCTION__+"emsg"); //--- return result return(EMPTY_VALUE); } sixbits[m_ser_entry_length]=0; //--- function call FourSixBits2ThreeBytes(sixbits,0,bytes,0); //--- function call FourSixBits2ThreeBytes(sixbits,4,bytes,3); //--- function call FourSixBits2ThreeBytes(sixbits,8,bytes,6); //--- copy for(i=0;i