//+------------------------------------------------------------------+ //| ap.mqh | //| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) | //| Copyright 2012-2026, MetaQuotes Ltd. | //| 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 "matrix.mqh" #include "bitconvert.mqh" //--- #define IsPosInf(value) (AL_POSINF==value) #define IsNegInf(value) (AL_NEGINF==value) //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ double AL_NaN =CInfOrNaN::NaN(); //---AL_NaN double AL_POSINF=CInfOrNaN::PositiveInfinity(); //---+infinity 1.#INF double AL_NEGINF=CInfOrNaN::NegativeInfinity(); //----infinity -1.#INF //+------------------------------------------------------------------+ //| Reverse communication structure | //+------------------------------------------------------------------+ struct RCommState { public: int stage; CRowInt ia; bool ba[]; CRowDouble ra; CRowComplex ca; RCommState(void) { stage=-1; } ~RCommState(void) {} void Copy(const RCommState &obj); //--- overloading void operator=(const RCommState &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| Create a copy | //+------------------------------------------------------------------+ void RCommState::Copy(const RCommState &obj) { //--- copy a variable stage=obj.stage; //--- copy arrays ia=obj.ia; ArrayCopy(ba,obj.ba); ra=obj.ra; ca=obj.ca; } //+------------------------------------------------------------------+ //| Internal functions | //+------------------------------------------------------------------+ enum TRACE_MODE { TRACE_NONE,TRACE_FILE }; //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ class CAp { public: //--- variable that determines whether an exception happened static bool exception_happened; //--- len template static int Len(const T &a[]); static int Len(const CRowInt &a); static int Len(const CRowDouble &a); static int Len(const CRowComplex &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 template static void Swap(T &a,T &b); template static void Swap(T &a[],T &b[]); static void Swap(vector &a,vector &b); static void Swap(vector &a,vector &b); static void Swap(vector &a,vector &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); static bool SetErrorFlag(bool &flag,bool cond,string xdesc); static void SetErrorFlagDiff(bool &flag,double val,double refval,double tol,double 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 vector &a,const int dps); static string Format(const complex &a[],const int dps); static string Format(const vector &a,const int dps); static string FormatB(CMatrixInt &a); static string Format(CMatrixInt &a); static string Format(CMatrixDouble &a,const int dps); static string Format(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); //--- trace static void TraceFile(string tags,string filename); static void TraceDisable(void); static bool IsTraceEnabled(string tag); static void Trace(string s); private: static TRACE_MODE m_TraceMode; static string m_TraceTags; static string m_TraceFilename; }; //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ TRACE_MODE CAp::m_TraceMode=TRACE_NONE; string CAp::m_TraceTags=""; string CAp::m_TraceFilename=""; //+------------------------------------------------------------------+ //| Initialize variable | //+------------------------------------------------------------------+ bool CAp::exception_happened=false; //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ template int CAp::Len(const T &a[]) { return(ArraySize(a)); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ int CAp::Len(const CRowInt &a) { return((int)a.Size()); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ int CAp::Len(const CRowDouble &a) { return((int)a.Size()); } //+------------------------------------------------------------------+ //| Get array lenght | //+------------------------------------------------------------------+ int CAp::Len(const CRowComplex &a) { return((int)a.Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ int CAp::Rows(const CMatrixInt &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ int CAp::Rows(const CMatrixDouble &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ int CAp::Rows(const CMatrixComplex &a) { return(a.Size()); } //+------------------------------------------------------------------+ //| Get cols count | //+------------------------------------------------------------------+ int CAp::Cols(const CMatrixInt &a) { //--- check if(a.Size()==0) return(0); //--- return result return(a[0].Size()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ int CAp::Cols(const CMatrixDouble &a) { return(a.Cols()); } //+------------------------------------------------------------------+ //| Get rows count | //+------------------------------------------------------------------+ int CAp::Cols(const CMatrixComplex &a) { return(a.Cols()); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ template void CAp::Swap(T &a,T &b) { T t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ template void CAp::Swap(T &a[],T &b[]) { ArraySwap(a,b); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(vector &a,vector &b) { a.Swap(b); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(vector &a,vector &b) { a.Swap(b); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(vector &a,vector &b) { a.Swap(b); } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(CMatrixInt &a,CMatrixInt &b) { //--- create matrix CMatrixInt t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(CMatrixDouble &a,CMatrixDouble &b) { //--- create matrix CMatrixDouble t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Swap | //+------------------------------------------------------------------+ void CAp::Swap(CMatrixComplex &a,CMatrixComplex &b) { //--- create matrix CMatrixComplex t; //--- swap t=a; a=b; b=t; } //+------------------------------------------------------------------+ //| Check assertions | //+------------------------------------------------------------------+ bool CAp::Assert(const bool cond) { return(Assert(cond,"ALGLIB: assertion failed")); } //+------------------------------------------------------------------+ //| Check assertions | //+------------------------------------------------------------------+ bool CAp::Assert(const bool cond,const string s) { //--- check if(!cond) { Print(__FUNCTION__+" " + s); exception_happened=true; return(false); } //--- the assertion is true return(true); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ bool CAp::SetErrorFlag(bool &flag,bool cond,string xdesc) { if(cond) { flag=true; Print(xdesc); } //--- return(flag); } //+------------------------------------------------------------------+ //| Internally calls SetErrorFlag() with condition: | //| Abs(Val-RefVal)>Tol*Max(Abs(RefVal),S) | //| This function is used to test relative error in Val against | //| RefVal, with relative error being replaced by absolute when scale| //| of RefVal is less than S. | //| This function returns value of COND. | //+------------------------------------------------------------------+ void CAp::SetErrorFlagDiff(bool &flag,double val,double refval, double tol,double s) { CAp::SetErrorFlag(flag,MathAbs(val-refval)>(tol*MathMax(MathAbs(refval),s)),"apserv.ap:162"); } //+------------------------------------------------------------------+ //| 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 | //+------------------------------------------------------------------+ 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 | //+------------------------------------------------------------------+ 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.imag>=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.real) + sign + StringFormat("%" + fmty,MathAbs(a.imag)) + "i"; StringReplace(res,",","."); //--- return result return(res); } //+------------------------------------------------------------------+ //| Prints formatted array | //+------------------------------------------------------------------+ 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 &a,const int dps) { ulong size=a.Size(); //--- check if(size==0) { Print(__FUNCTION__+": vector size error"); return(NULL); } string result[]; ArrayResizeAL(result,(int)size); //--- definition of output style string sfmt; if(dps>=0) sfmt="f"; else sfmt="e"; //--- converting int d=(int)MathAbs(dps); string fmt=StringFormat(".%d" + sfmt,d); for(ulong 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].real) + sign + StringFormat("%" + fmty,MathAbs(a[i].imag)) + "i"; StringReplace(result[i],",","."); } //--- return result return("{" + StringJoin(",",result) + "}"); } //+------------------------------------------------------------------+ //| Prints formatted matrix | //+------------------------------------------------------------------+ string CAp::FormatB(CMatrixInt &a) { int m=a.Rows(); //--- check if(m==0) { Print(__FUNCTION__+": array size error"); return(NULL); } int n=a.Cols(); //--- 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; i v=a[i]; result[i]=Format(v,dps); } //--- return result return("{" + StringJoin(",",result) + "}"); } //+------------------------------------------------------------------+ //| Prints formatted matrix | //+------------------------------------------------------------------+ string CAp::Format(CMatrixComplex &a,const int dps) { ulong m=a.Rows(); //--- check if(m==0) { Print(__FUNCTION__+": matrix size error"); return(NULL); } ulong n=a.Cols(); //--- check if(n==0) { Print(__FUNCTION__+": matrix size error"); return(NULL); } //--- prepare arrays string result[]; ArrayResizeAL(result,(int)m); //--- converting for(ulong i=0; i v=a[i]; result[i]=Format(v,dps); } //--- return result return("{" + StringJoin(",",result) + "}"); } //+------------------------------------------------------------------+ //| checks that matrix is symmetric. | //| max|A-A^T| is calculated; if it is within 1.0E-14 of max|A|, | //| matrix is considered symmetric | //+------------------------------------------------------------------+ bool CAp::IsSymmetric(const CMatrixDouble &a) { ulong n=a.Rows(); //--- check if(n!=a.Cols()) return(false); //--- check if(n==0) return(true); //--- check for symmetry double mx=0; double err=0; for(ulong i=0; i m=a.Transpose(); a=a.TriL()+m.TriU(1); //--- return result return(true); } //+------------------------------------------------------------------+ //| Forces Hermiticity by copying upper half of A to the lower one | //+------------------------------------------------------------------+ bool CAp::ForceHermitian(CMatrixComplex &a) { ulong n=a.Rows(); //--- check if(n!=a.Cols()) return(false); //--- check if(n==0) return(true); //--- change matrix complex c; for(ulong i=0; i= 0) return(true); //--- contains tag (followed by dot, which means match with child) if(StringFind(m_TraceTags,("," + tag + ".")) >= 0) return(true); //---nothing return(false); } //+------------------------------------------------------------------+ //| Trace | //+------------------------------------------------------------------+ void CAp::Trace(string s) { if(m_TraceMode==TRACE_FILE) { int handle=FileOpen(m_TraceFilename,FILE_READ|FILE_WRITE|FILE_TXT|FILE_ANSI); if(handle<=0 || !FileSeek(handle,0,SEEK_END)) return; FileWriteString(handle,s); FileClose(handle); } } //+------------------------------------------------------------------+ //| Portable high quality random number generator state. | //| Initialized with HQRNDRandomize() or HQRNDSeed(). | //| Fields: | //| S1, S2 - seed values | //| V - precomputed value | //| MagicV - 'magic' value used to determine whether State| //| structure was correctly initialized. | //+------------------------------------------------------------------+ class CHighQualityRandState { public: //--- variables int m_s1; int m_s2; double m_v; int m_magicv; //--- constructor, destructor CHighQualityRandState(void) { ZeroMemory(this); } ~CHighQualityRandState(void) {} //--- void Copy(const CHighQualityRandState &obj); //--- overloading void operator=(const CHighQualityRandState &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CHighQualityRandState::Copy(const CHighQualityRandState &obj) { m_s1=obj.m_s1; m_s2=obj.m_s2; m_magicv=obj.m_magicv; m_v=obj.m_v; } //+------------------------------------------------------------------+ //| Portable high quality random number generator state. | //| Initialized with HQRNDRandomize() or HQRNDSeed(). | //| Fields: | //| S1, S2 - seed values | //| V - precomputed value | //| MagicV - 'magic' value used to determine whether State| //| structure was correctly initialized. | //+------------------------------------------------------------------+ class CHighQualityRandStateShell { private: CHighQualityRandState m_innerobj; public: //--- constructors, destructor CHighQualityRandStateShell(void) {} CHighQualityRandStateShell(CHighQualityRandState &obj); ~CHighQualityRandStateShell(void) {} //--- method CHighQualityRandState *GetInnerObj(void) { return(GetPointer(m_innerobj)); } }; //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ CHighQualityRandStateShell::CHighQualityRandStateShell(CHighQualityRandState &obj) { //--- copy m_innerobj.m_s1=obj.m_s1; m_innerobj.m_s2=obj.m_s2; m_innerobj.m_v=obj.m_v; m_innerobj.m_magicv=obj.m_magicv; } //+------------------------------------------------------------------+ //| Portable high quality random number generator | //+------------------------------------------------------------------+ class CHighQualityRand { public: //--- static class members static const int m_HQRndMax; static const int m_HQRndM1; static const int m_HQRndM2; static const int m_HQRndMagic; //--- public methods static void HQRndRandomize(CHighQualityRandState &state); static void HQRndSeed(const int s1,const int s2,CHighQualityRandState &state); static double HQRndUniformR(CHighQualityRandState &state); static int HQRndUniformI(CHighQualityRandState &state,const int n); static double HQRndNormal(CHighQualityRandState &state); static void HQRndNormalV(CHighQualityRandState &state,int n,CRowDouble &x); static void HQRndNormalM(CHighQualityRandState &state,int m,int n,CMatrixDouble &x); static void HQRndUnit2(CHighQualityRandState &state,double &x,double &y); static void HQRndNormal2(CHighQualityRandState &state,double &x1,double &x2); static double HQRndExponential(CHighQualityRandState &state,const double lambdav); static double HQRndDiscrete(CHighQualityRandState &state,int n,CRowDouble &x); template static T HQRndDiscrete(CHighQualityRandState &state,int n,vector &x); static double HQRndContinuous(CHighQualityRandState &state,int n,CRowDouble &x); template static T HQRndContinuous(CHighQualityRandState &state,int n,vector &x); private: static int HQRndIntegerBase(CHighQualityRandState &state); }; //+------------------------------------------------------------------+ //| Initialize constants | //+------------------------------------------------------------------+ const int CHighQualityRand::m_HQRndMax=2147483561; const int CHighQualityRand::m_HQRndM1 =2147483563; const int CHighQualityRand::m_HQRndM2 =2147483399; const int CHighQualityRand::m_HQRndMagic=1634357784; //+------------------------------------------------------------------+ //| HQRNDState initialization with random values which come from | //| standard RNG. | //+------------------------------------------------------------------+ void CHighQualityRand::HQRndRandomize(CHighQualityRandState &state) { //--- get result HQRndSeed(CMath::RandomInteger(m_HQRndM1),CMath::RandomInteger(m_HQRndM2),state); } //+------------------------------------------------------------------+ //| HQRNDState initialization with seed values | //+------------------------------------------------------------------+ void CHighQualityRand::HQRndSeed(const int s1,const int s2, CHighQualityRandState &state) { //--- calculation parameters state.m_s1=s1%(m_HQRndM1-1)+1; state.m_s2=s2%(m_HQRndM2-1)+1; state.m_v=1.0/(double)m_HQRndMax; state.m_magicv=m_HQRndMagic; } //+------------------------------------------------------------------+ //| This function generates random real number in [0,1). | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ double CHighQualityRand::HQRndUniformR(CHighQualityRandState &state) { return(state.m_v*(HQRndIntegerBase(state)-1)); } //+------------------------------------------------------------------+ //| This function generates random integer number in [0, N) | //| 1. N must be less than HQRNDMax-1. | //| 2. State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed() | //+------------------------------------------------------------------+ int CHighQualityRand::HQRndUniformI(CHighQualityRandState &state,const int n) { //--- create variables int result=0; int maxcnt=0; int mx=0; int a=0; int b=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0!")) return(INT_MAX); //--- maxcnt=m_HQRndMax+1; //---Two branches: one for N<=MaxCnt, another for N>MaxCnt. if(n>maxcnt) { //---N>=MaxCnt. //---We have two options here: //---a) N is exactly divisible by MaxCnt //---b) N is not divisible by MaxCnt //---In both cases we reduce problem on interval spanning [0,N) //---to several subproblems on intervals spanning [0,MaxCnt). if(n%maxcnt==0) { //---N is exactly divisible by MaxCnt. //---[0,N) range is dividided into N/MaxCnt bins, //---each of them having length equal to MaxCnt. //---We generate: //---* random bin number B //---* random offset within bin A //---Both random numbers are generated by recursively //---calling HQRNDUniformI(). //---Result is equal to A+MaxCnt*B. //--- check if(!CAp::Assert(n/maxcnt<=maxcnt,__FUNCTION__+": N is too large")) return(INT_MAX); a=HQRndUniformI(state,maxcnt); b=HQRndUniformI(state,n/maxcnt); result=a+maxcnt*b; } else { //---N is NOT exactly divisible by MaxCnt. //---[0,N) range is dividided into Ceil(N/MaxCnt) bins, //---each of them having length equal to MaxCnt. //---We generate: //---* random bin number B in [0, Ceil(N/MaxCnt)-1] //---* random offset within bin A //---* if both of what is below is true //--- 1) bin number B is that of the last bin //--- 2) A >= N mod MaxCnt //--- then we repeat generation of A/B. //--- This stage is essential in order to avoid bias in the result. //---* otherwise, we return A*MaxCnt+N //--- check if(!CAp::Assert(n/maxcnt+1<=maxcnt,__FUNCTION__+": N is too large")) return(INT_MAX); result=-1; do { a=HQRndUniformI(state,maxcnt); b=HQRndUniformI(state,n/maxcnt+1); if(b==n/maxcnt && a>=n%maxcnt) continue; result=a+maxcnt*b; } while(result<0); } } else { //---N<=MaxCnt //---Code below is a bit complicated because we can not simply //---return "HQRNDIntegerBase() mod N" - it will be skewed for //---large N's in [0.1*HQRNDMax...HQRNDMax]. mx=maxcnt-maxcnt%n; do { result=HQRndIntegerBase(state); } while(result>=mx); 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(). | //+------------------------------------------------------------------+ 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: vector with random entries (normal | //| distribution) | //| This function generates N random numbers from normal | //| distribution. | //| State structure must be initialized with HQRNDRandomize() or | //| HQRNDSeed(). | //+------------------------------------------------------------------+ void CHighQualityRand::HQRndNormalV(CHighQualityRandState &state, int n,CRowDouble &x) { //--- create variables int n2=0; double v1=0; double v2=0; //--- function call if(n<0) { n=(int)x.Size(); } else { if((int)x.Size()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(). | //+------------------------------------------------------------------+ 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)); } //+------------------------------------------------------------------+ //| This function generates random number from discrete distribution| //| given by finite sample X. | //| INPUT PARAMETERS | //| State - high quality random number generator, must be | //| initialized with HQRNDRandomize() or HQRNDSeed(). | //| X - finite sample | //| N - number of elements to use, N>=1 | //| RESULT | //| this function returns one of the X[i] for random i=0..N-1 | //+------------------------------------------------------------------+ double CHighQualityRand::HQRndDiscrete(CHighQualityRandState &state, int n,CRowDouble &x) { double result=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0")) return (EMPTY_VALUE); if(!CAp::Assert(n<=(int)x.Size(),__FUNCTION__+": Length(X)=1 | //| RESULT | //| this function returns one of the X[i] for random i=0..N-1 | //+------------------------------------------------------------------+ template T CHighQualityRand::HQRndDiscrete(CHighQualityRandState &state, int n,vector &x) { T result=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0")) return (EMPTY_VALUE); if(!CAp::Assert(n<=(int)x.Size(),__FUNCTION__+": Length(X)=1 | //| RESULT | //| this function returns random number from continuous | //| distribution which tries to approximate X as mush as possible. | //| min(X)<=Result<=max(X). | //+------------------------------------------------------------------+ double CHighQualityRand::HQRndContinuous(CHighQualityRandState &state, int n,CRowDouble &x) { //--- create variables double result=0; double mx=0; double mn=0; int i=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0")) return (EMPTY_VALUE); if(!CAp::Assert(n<=(int)x.Size(),__FUNCTION__+": Length(X)=mn,__FUNCTION__+": X is not sorted by ascending")) return(EMPTY_VALUE); if(mx!=mn) result=(mx-mn)*HQRndUniformR(state)+mn; else result=mn; //--- return result return(result); } //+------------------------------------------------------------------+ //| This function generates random number from continuous | //| distribution given by finite sample X. | //| INPUT PARAMETERS | //| State - high quality random number generator, must be | //| initialized with HQRNDRandomize() or HQRNDSeed(). | //| X - finite sample, array[N] (can be larger, in this | //| case only leading N elements are used). THIS ARRAY | //| MUST BE SORTED BY ASCENDING. | //| N - number of elements to use, N>=1 | //| RESULT | //| this function returns random number from continuous | //| distribution which tries to approximate X as mush as possible. | //| min(X)<=Result<=max(X). | //+------------------------------------------------------------------+ template T CHighQualityRand::HQRndContinuous(CHighQualityRandState &state, int n,vector &x) { //--- create variables T result=0; double mx=0; double mn=0; int i=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0")) return (EMPTY_VALUE); if(!CAp::Assert(n<=(int)x.Size(),__FUNCTION__+": Length(X)=mn,__FUNCTION__+": X is not sorted by ascending")) return(EMPTY_VALUE); if(mx!=mn) result=(T)((mx-mn)*HQRndUniformR(state)+mn); else result=mn; //--- return result return(result); } //+------------------------------------------------------------------+ //| L'Ecuyer, Efficient and portable combined random number | //| generators | //+------------------------------------------------------------------+ 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; //--- 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; //+------------------------------------------------------------------+ //| Check on +-inf | //+------------------------------------------------------------------+ bool CMath::IsFinite(const double d) { return(MathIsValidNumber(d)); } //+------------------------------------------------------------------+ //| Random real value [0,1) | //+------------------------------------------------------------------+ double CMath::RandomReal(void) { 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 | //+------------------------------------------------------------------+ 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 | //+------------------------------------------------------------------+ double CMath::AbsComplex(const complex z) { //--- initialization double xabs=MathAbs(z.real); double yabs=MathAbs(z.imag); //--- check if(xabs==0) return(yabs); //--- check if(yabs==0) return(xabs); //--- check //--- calculation double t=MathSqrt(MathPow(xabs,2)+MathPow(yabs,2)); //--- return result return(t); } //+------------------------------------------------------------------+ //| The absolute value of a complex number | //+------------------------------------------------------------------+ double CMath::AbsComplex(const double r) { //--- initialization complex z=r; double w=0.0; double v=0.0; double xabs=MathAbs(z.real); double yabs=MathAbs(z.imag); //--- check if(xabs>yabs) w=xabs; else w=yabs; //--- check if(xabs0?1:0); m_entries_needed+=1+n; } //+------------------------------------------------------------------+ //| Switching state on TO_STRING | //+------------------------------------------------------------------+ void CSerializer::SStart_Str(void) { //--- get size int allocsize=Get_Alloc_Size(); //--- clear input/output buffers which may hold pointers to unneeded memory //--- NOTE: it also helps us to avoid errors when data are written to incorrect location ClearBuffers(); //--- check and change m_mode if(m_mode!=ALLOC) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- exit the function return; } m_mode=TO_STRING; //--- other preparations ArrayResizeAL(m_out_str,allocsize); m_entries_saved=0; m_bytes_written=0; } //+------------------------------------------------------------------+ //| Switching state on FROM_STRING | //+------------------------------------------------------------------+ void CSerializer::UStart_Str(const string s) { //--- check and change m_mode if(m_mode!=DEFAULT) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- exit the function return; } m_mode=FROM_STRING; StringToCharArray(s,m_in_str); m_bytes_read=0; } //+------------------------------------------------------------------+ //| Reset parameters | //+------------------------------------------------------------------+ void CSerializer::Reset(void) { m_mode=DEFAULT; m_entries_needed=0; m_bytes_asked=0; ClearBuffers(); } //+------------------------------------------------------------------+ //| Serialize bool | //+------------------------------------------------------------------+ void CSerializer::Serialize_Bool(const bool v) { //--- check if(m_mode!=TO_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- exit the function return; } //--- function call Bool2Str(v,m_out_str,m_bytes_written); m_entries_saved++; //--- check if(m_entries_saved%m_ser_entries_per_row!=0) { m_out_str[m_bytes_written]=' '; m_bytes_written++; } else { m_out_str[m_bytes_written+0]='\r'; m_out_str[m_bytes_written+1]='\n'; m_bytes_written+=2; } } //+------------------------------------------------------------------+ //| Serialize int | //+------------------------------------------------------------------+ void CSerializer::Serialize_Int(const int v) { //--- check if(m_mode!=TO_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- exit the function return; } //--- function call Int2Str(v,m_out_str,m_bytes_written); m_entries_saved++; //--- check if(m_entries_saved%m_ser_entries_per_row!=0) { m_out_str[m_bytes_written]=' '; m_bytes_written++; } else { m_out_str[m_bytes_written+0]='\r'; m_out_str[m_bytes_written+1]='\n'; m_bytes_written+=2; } } //+------------------------------------------------------------------+ //| Serialize double | //+------------------------------------------------------------------+ void CSerializer::Serialize_Double(const double v) { //--- check if(m_mode!=TO_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- exit the function return; } //--- function call Double2Str(v,m_out_str,m_bytes_written); m_entries_saved++; //--- check if(m_entries_saved%m_ser_entries_per_row!=0) { m_out_str[m_bytes_written]=' '; m_bytes_written++; } else { m_out_str[m_bytes_written+0]='\r'; m_out_str[m_bytes_written+1]='\n'; m_bytes_written+=2; } } //+------------------------------------------------------------------+ //| Unserialize bool | //+------------------------------------------------------------------+ bool CSerializer::Unserialize_Bool(void) { //--- check if(m_mode!=FROM_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- return result return(false); } //--- return result return(Str2Bool(m_in_str,m_bytes_read)); } //+------------------------------------------------------------------+ //| Unserialize int | //+------------------------------------------------------------------+ int CSerializer::Unserialize_Int(void) { //--- check if(m_mode!=FROM_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- return result return(-1); } //--- return result return(Str2Int(m_in_str,m_bytes_read)); } //+------------------------------------------------------------------+ //| Unserialize double | //+------------------------------------------------------------------+ double CSerializer::Unserialize_Double(void) { //--- check if(m_mode!=FROM_STRING) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- return result return(EMPTY_VALUE); } //--- return result return(Str2Double(m_in_str,m_bytes_read)); } //+------------------------------------------------------------------+ //| Stop | //+------------------------------------------------------------------+ void CSerializer::Stop(void) { if(m_mode==TO_STRING) { m_out_str[m_bytes_written]='.'; m_bytes_written++; return; } if(m_mode==FROM_STRING) { //--- because input string may be from pre-3.11 serializer, //--- which does not include trailing dot, we do not test //--- string for presence of "." symbol. Anyway, because string //--- is not stream, we do not have to read ALL trailing symbols. return; } Print(__FUNCTION__": internal error during unserialization"); } //+------------------------------------------------------------------+ //| Get string | //+------------------------------------------------------------------+ string CSerializer::Get_String(void) { return(GetSelectionString(m_out_str,0,m_bytes_written)); } //+------------------------------------------------------------------+ //| Get alloc size | //+------------------------------------------------------------------+ int CSerializer::Get_Alloc_Size(void) { //--- create variables int rows; int lastrowsize; int result; //--- check and change m_mode if(m_mode!=ALLOC) { Print(__FUNCTION__+": internal error during (un)serialization"); //--- return result return(-1); } //--- if no entries needes (degenerate case) if(m_entries_needed==0) { m_bytes_asked=1; //--- return result return(m_bytes_asked); } //--- non-degenerate case rows=m_entries_needed/m_ser_entries_per_row; lastrowsize=m_ser_entries_per_row; //--- check if(m_entries_needed%m_ser_entries_per_row!=0) { lastrowsize=m_entries_needed%m_ser_entries_per_row; rows++; } //--- calculate result size result=((rows-1)*m_ser_entries_per_row+lastrowsize)*m_ser_entry_length; result+=(rows-1)*(m_ser_entries_per_row-1)+(lastrowsize-1); result+=rows*2; //--- save result m_bytes_asked=result; //--- return result return(result); } //+------------------------------------------------------------------+ //| This function converts six-bit value (from 0 to 63) to character | //| (only digits, lowercase and uppercase letters, minus and | //| underscore are used). If v is negative or greater than 63, this | //| function returns '?'. | //+------------------------------------------------------------------+ char CSerializer::SixBits2Char(const int v) { //--- check if(v<0 || v>63) 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. | //+------------------------------------------------------------------+ 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 | //+------------------------------------------------------------------+ 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 | //+------------------------------------------------------------------+ 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. | //+------------------------------------------------------------------+ void CSerializer::Bool2Str(const bool v,char &buf[],int &offs) { char c; //--- check if(v) c='1'; else c='0'; //--- copy c for(int 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(int i=0; i