//+------------------------------------------------------------------+
//|                                               fasttransforms.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 "alglibinternal.mqh"
//+------------------------------------------------------------------+
//| Fast Fourier Transform                                           |
//+------------------------------------------------------------------+
class CFastFourierTransform
  {
public:
   static void       FFTC1D(complex &a[],const int n);
   static void       FFTC1DInv(complex &a[],const int n);
   static void       FFTR1D(double &a[],const int n,complex &f[]);
   static void       FFTR1DInv(complex &f[],const int n,double &a[]);
   static void       FFTR1DInternalEven(double &a[],const int n,double &buf[],CFtPlan &plan);
   static void       FFTR1DInvInternalEven(double &a[],const int n,double &buf[],CFtPlan &plan);
   static void       FFTC1D(CRowComplex &a,const int n);
   static void       FFTC1DInv(CRowComplex &a,const int n);
   static void       FFTR1D(CRowDouble &a,const int n,CRowComplex &f);
   static void       FFTR1DInv(CRowComplex &f,const int n,CRowDouble &a);
   static void       FFTR1DInternalEven(CRowDouble &a,const int n,CRowDouble &buf,CFtPlan &plan);
   static void       FFTR1DInvInternalEven(CRowDouble &a,const int n,CRowDouble &buf,CFtPlan &plan);
  };
//+------------------------------------------------------------------+
//| 1-dimensional complex FFT.                                       |
//| Array size N may be arbitrary number (composite or prime).       |
//| Composite N's are handled with cache-oblivious variation of a    |
//| Cooley-Tukey algorithm. Small prime-factors are transformed using|
//| hard coded codelets (similar to FFTW codelets, but without       |
//| low-level optimization), large prime-factors are handled with    |
//| Bluestein's algorithm.                                           |
//| Fastests transforms are for smooth N's (prime factors are 2, 3,  |
//| 5 only), most fast for powers of 2. When N have prime factors    |
//| larger than these, but orders of magnitude smaller than N,       |
//| computations will be about 4 times slower than for nearby highly |
//| composite N's. When N itself is prime, speed will be 6 times     |
//| lower.                                                           |
//| Algorithm has O(N*logN) complexity for any N (composite or       |
//| prime).                                                          |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..N-1] - complex function to be transformed   |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     A   -   DFT of a input array, array[0..N-1]                  |
//|             A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N),    |
//|             k = 0..N-1)                                          |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTC1D(complex &a[],const int n)
  {
//--- create a variable
   int i=0;
//--- create array
   CRowDouble buf;
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Special case: N=1,FFT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
      return;
//--- convert input array to the more convinient format
   buf.Resize(2*n);
   for(i=0; i<n; i++)
     {
      buf.Set(2*i+0,a[i].real);
      buf.Set(2*i+1,a[i].imag);
     }
//--- Generate plan and execute it.
//--- Plan is a combination of a successive factorizations of N and
//--- precomputed data. It is much like a FFTW plan,but is not stored
//--- between subroutine calls and is much simpler.
   CFtBase::FtComplexFFTPlan(n,1,plan);
   CFtBase::FtApplyPlan(plan,buf,0,1);
//--- result
   for(i=0; i<n; i++)
     {
      a[i].real=buf[2*i+0];
      a[i].imag=buf[2*i+1];
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTC1D(CRowComplex &a,const int n)
  {
//--- create a variable
   int i=0;
//--- create array
   CRowDouble buf;
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Special case: N=1,FFT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
      return;
//--- convert input array to the more convinient format
   buf.Resize(2*n);
   for(i=0; i<n; i++)
     {
      buf.Set(2*i+0,a[i].real);
      buf.Set(2*i+1,a[i].imag);
     }
//--- Generate plan and execute it.
//--- Plan is a combination of a successive factorizations of N and
//--- precomputed data. It is much like a FFTW plan,but is not stored
//--- between subroutine calls and is much simpler.
   CFtBase::FtComplexFFTPlan(n,1,plan);
   CFtBase::FtApplyPlan(plan,buf,0,1);
//--- result
   for(i=0; i<n; i++)
     {
      a.SetRe(i,buf[2*i+0]);
      a.SetIm(i,buf[2*i+1]);
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional complex inverse FFT.                               |
//| Array size N may be arbitrary number (composite or prime).       |
//| Algorithm has O(N*logN) complexity for any N (composite or prime)|
//| See FFTC1D() description for more information about algorithm    |
//| performance.                                                     |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..N-1] - complex array to be transformed      |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     A   -   inverse DFT of a input array, array[0..N-1]          |
//|             A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N),  |
//|             k = 0..N-1)                                          |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTC1DInv(complex &a[],const int n)
  {
//--- create a variable
   int i=0;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Inverse DFT can be expressed in terms of the DFT as
//---     invfft(x)=fft(x')'/N
//--- here x' means conj(x).
   for(i=0; i<=n-1; i++)
      a[i].imag=-a[i].imag;
//--- function call
   FFTC1D(a,n);
//--- change values
   for(i=0; i<=n-1; i++)
     {
      a[i].real=a[i].real/n;
      a[i].imag=-(a[i].imag/n);
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTC1DInv(CRowComplex &a,const int n)
  {
//--- create a variable
   int i=0;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Inverse DFT can be expressed in terms of the DFT as
//---     invfft(x)=fft(x')'/N
//--- here x' means conj(x).
   for(i=0; i<n; i++)
      a.SetIm(i,-a[i].imag);
//--- function call
   FFTC1D(a,n);
//--- change values
   for(i=0; i<=n-1; i++)
     {
      a.SetRe(i,a[i].real/n);
      a.SetIm(i,-(a[i].imag/n));
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional real FFT.                                          |
//| Algorithm has O(N*logN) complexity for any N (composite or       |
//| prime).                                                          |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..N-1] - real function to be transformed      |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     F   -   DFT of a input array, array[0..N-1]                  |
//|             F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N),           |
//|             k = 0..N-1)                                          |
//| NOTE:                                                            |
//|     F[] satisfies symmetry property F[k] = conj(F[N-k]), so just |
//| one half of array is usually needed. But for convinience         |
//| subroutine returns full complex array (with frequencies above    |
//| N/2), so its result may be used by other FFT-related subroutines.|
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1D(double &a[],const int n,complex &f[])
  {
//--- create variables
   int     i=0;
   int     n2=0;
   int     idx=0;
   complex hn=0;
   complex hmnc=0;
   complex v=0;
   int     i_=0;
//--- create array
   CRowDouble buf;
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Special cases:
//--- * N=1,FFT is just identity transform.
//--- * N=2,FFT is simple too
//--- After this block we assume that N is strictly greater than 2
   if(n==1)
     {
      //--- allocation
      ArrayResize(f,1);
      f[0]=a[0];
      //--- exit the function
      return;
     }
//--- check
   if(n==2)
     {
      //--- allocation
      ArrayResize(f,2);
      f[0].real=a[0]+a[1];
      f[0].imag=0;
      f[1].real=a[0]-a[1];
      f[1].imag=0;
      //--- exit the function
      return;
     }
//--- Choose between odd-size and even-size FFTs
   if(n%2==0)
     {
      //--- even-size real FFT,use reduction to the complex task
      n2=n/2;
      //--- allocation
      buf=a;
      buf.Resize(n);
      //--- function call
      CFtBase::FtComplexFFTPlan(n2,1,plan);
      CFtBase::FtApplyPlan(plan,buf,0,1);
      //--- allocation
      ArrayResize(f,n);
      //--- calculation
      for(i=0; i<=n2; i++)
        {
         idx=2*(i%n2);
         hn.real=buf[idx+0];
         hn.imag=buf[idx+1];
         idx=2*((n2-i)%n2);
         hmnc.real=buf[idx+0];
         hmnc.imag=-buf[idx+1];
         v.real=-MathSin(-(2*M_PI*i/n));
         v.imag=MathCos(-(2*M_PI*i/n));
         f[i]=hn+hmnc-v*(hn-hmnc);
         f[i].real=0.5*f[i].real;
         f[i].imag=0.5*f[i].imag;
        }
      //--- copy
      for(i=n2+1; i<=n-1; i++)
         f[i]=CMath::Conj(f[n-i]);
     }
   else
     {
      //--- use complex FFT
      ArrayResize(f,n);
      //--- copy
      for(i=0; i<=n-1; i++)
         f[i]=a[i];
      //--- function call
      FFTC1D(f,n);
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1D(CRowDouble &a,const int n,CRowComplex &f)
  {
//--- create variables
   int     i=0;
   int     n2=0;
   int     idx=0;
   complex hn=0;
   complex hmnc=0;
   complex v=0;
   int     i_=0;
//--- create array
   CRowDouble buf;
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
      return;
//--- Special cases:
//--- * N=1,FFT is just identity transform.
//--- * N=2,FFT is simple too
//--- After this block we assume that N is strictly greater than 2
   if(n==1)
     {
      //--- allocation
      f.Resize(1);
      f.Set(0,a[0]);
      //--- exit the function
      return;
     }
//--- check
   if(n==2)
     {
      //--- allocation
      f.Resize(2);
      f.Set(0,a[0]+a[1]);
      f.Set(1,a[0]-a[1]);
      //--- exit the function
      return;
     }
//--- Choose between odd-size and even-size FFTs
   if(n%2==0)
     {
      //--- even-size real FFT,use reduction to the complex task
      n2=n/2;
      //--- allocation
      buf=a;
      buf.Resize(n);
      //--- function call
      CFtBase::FtComplexFFTPlan(n2,1,plan);
      CFtBase::FtApplyPlan(plan,buf,0,1);
      //--- allocation
      f.Resize(n);
      //--- calculation
      for(i=0; i<=n2; i++)
        {
         idx=2*(i%n2);
         hn.real=buf[idx+0];
         hn.imag=buf[idx+1];
         idx=2*((n2-i)%n2);
         hmnc.real=buf[idx+0];
         hmnc.imag=-buf[idx+1];
         v.real=-MathSin(-(2*M_PI*i/n));
         v.imag=MathCos(-(2*M_PI*i/n));
         f.Set(i,hn+hmnc-v*(hn-hmnc));
         f.SetRe(i,0.5*f[i].real);
         f.SetIm(i,0.5*f[i].imag);
        }
      //--- copy
      for(i=n2+1; i<n; i++)
         f.Set(i,CMath::Conj(f[n-i]));
     }
   else
     {
      //--- use complex FFT
      f.Resize(n);
      //--- copy
      for(i=0; i<n; i++)
         f.Set(i,a[i]);
      //--- function call
      FFTC1D(f,n);
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional real inverse FFT.                                  |
//| Algorithm has O(N*logN) complexity for any N (composite or       |
//| prime).                                                          |
//| INPUT PARAMETERS                                                 |
//|     F   -   array[0..floor(N/2)] - frequencies from forward real |
//|             FFT                                                  |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     A   -   inverse DFT of a input array, array[0..N-1]          |
//| NOTE:                                                            |
//|     F[] should satisfy symmetry property F[k] = conj(F[N-k]),    |
//| so just one half of frequencies array is needed - elements from 0|
//| to floor(N/2). F[0] is ALWAYS real. If N is even F[floor(N/2)] is|
//| real too. If N is odd, then F[floor(N/2)] has no special         |
//| properties.                                                      |
//| Relying on properties noted above, FFTR1DInv subroutine uses only|
//| elements from 0th to floor(N/2)-th. It ignores imaginary part of |
//| F[0], and in case N is even it ignores imaginary part of         |
//| F[floor(N/2)] too.                                               |
//| When you call this function using full arguments list -          |
//| "FFTR1DInv(F,N,A)"                                               |
//| - you can pass either either frequencies array with N elements or|
//| reduced array with roughly N/2 elements - subroutine will        |
//| successfully transform both.                                     |
//| If you call this function using reduced arguments list -         |
//| "FFTR1DInv(F,A)" - you must pass FULL array with N elements      |
//| (although higher N/2 are still not used) because array size is   |
//| used to automatically determine FFT length                       |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInv(complex &f[],const int n,double &a[])
  {
//--- create a variable
   int i=0;
//--- create arrays
   double  h[];
   complex fh[];
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(f)>=(int)MathFloor((double)n/2.0)+1,__FUNCTION__+": Length(F)<Floor(N/2)+1!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(f[0].real),__FUNCTION__+": F contains infinite or NAN values!"))
      return;
   for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
     {
      //--- check
      if(!CAp::Assert(CMath::IsFinite(f[i].real) && CMath::IsFinite(f[i].imag),__FUNCTION__+": F contains infinite or NAN values!"))
         return;
     }
//--- check
   if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].real),__FUNCTION__+": F contains infinite or NAN values!"))
      return;
//--- check
   if(n%2!=0)
     {
      //--- check
      if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].imag),__FUNCTION__+": F contains infinite or NAN values!"))
         return;
     }
//--- Special case: N=1,FFT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
     {
      //--- allocation
      ArrayResize(a,1);
      a[0]=f[0].real;
      //--- exit the function
      return;
     }
//--- inverse real FFT is reduced to the inverse real FHT,
//--- which is reduced to the forward real FHT,
//--- which is reduced to the forward real FFT.
//--- Don't worry,it is really compact and efficient reduction :)
   ArrayResize(h,n);
   ArrayResize(a,n);
//--- change values
   h[0]=f[0].real;
   for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
     {
      h[i]=f[i].real-f[i].imag;
      h[n-i]=f[i].real+f[i].imag;
     }
//--- check
   if(n%2==0)
      h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real;
   else
     {
      h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real-f[(int)MathFloor((double)n/2.0)].imag;
      h[(int)MathFloor((double)n/2.0)+1]=f[(int)MathFloor((double)n/2.0)].real+f[(int)MathFloor((double)n/2.0)].imag;
     }
//--- function call
   FFTR1D(h,n,fh);
//--- change values
   for(i=0; i<=n-1; i++)
      a[i]=(fh[i].real-fh[i].imag)/n;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInv(CRowComplex &f,const int n,CRowDouble &a)
  {
//--- create a variable
   int i=0;
//--- create arrays
   double  h[];
   complex fh[];
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(f)>=(int)MathFloor((double)n/2.0)+1,__FUNCTION__+": Length(F)<Floor(N/2)+1!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(f[0].real),__FUNCTION__+": F contains infinite or NAN values!"))
      return;
   for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
     {
      //--- check
      if(!CAp::Assert(CMath::IsFinite(f[i].real) && CMath::IsFinite(f[i].imag),__FUNCTION__+": F contains infinite or NAN values!"))
         return;
     }
//--- check
   if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].real),__FUNCTION__+": F contains infinite or NAN values!"))
      return;
//--- check
   if(n%2!=0)
     {
      //--- check
      if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].imag),__FUNCTION__+": F contains infinite or NAN values!"))
         return;
     }
//--- Special case: N=1,FFT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
     {
      //--- allocation
      a.Resize(1);
      a.Set(0,f[0].real);
      //--- exit the function
      return;
     }
//--- inverse real FFT is reduced to the inverse real FHT,
//--- which is reduced to the forward real FHT,
//--- which is reduced to the forward real FFT.
//--- Don't worry,it is really compact and efficient reduction :)
   ArrayResize(h,n);
   a.Resize(n);
//--- change values
   h[0]=f[0].real;
   for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
     {
      h[i]=f[i].real-f[i].imag;
      h[n-i]=f[i].real+f[i].imag;
     }
//--- check
   if(n%2==0)
      h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real;
   else
     {
      h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real-f[(int)MathFloor((double)n/2.0)].imag;
      h[(int)MathFloor((double)n/2.0)+1]=f[(int)MathFloor((double)n/2.0)].real+f[(int)MathFloor((double)n/2.0)].imag;
     }
//--- function call
   FFTR1D(h,n,fh);
//--- change values
   for(i=0; i<n; i++)
      a.Set(i,(fh[i].real-fh[i].imag)/n);
  }
//+------------------------------------------------------------------+
//| Internal subroutine. Never call it directly!                     |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInternalEven(double &a[],const int n,
                                               double &buf[],CFtPlan &plan)
  {
   CRowDouble A=a;
   CRowDouble Buf;
   FFTR1DInternalEven(A,n,Buf,plan);
   A.ToArray(a);
   Buf.ToArray(buf);
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInternalEven(CRowDouble &a,const int n,
                                               CRowDouble &buf,CFtPlan &plan)
  {
//--- create variables
   double  x=0;
   double  y=0;
   int     i=0;
   int     n2=0;
   int     idx=0;
   complex hn=0;
   complex hmnc=0;
   complex v=0;
   int     i_=0;
//--- check
   if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
      return;
//--- Special cases:
//--- * N=2
//--- After this block we assume that N is strictly greater than 2
   if(n==2)
     {
      //--- change values
      x=a[0]+a[1];
      y=a[0]-a[1];
      a.Set(0,x);
      a.Set(1,y);
      //--- exit the function
      return;
     }
//--- even-size real FFT,use reduction to the complex task
   n2=n/2;
//--- copy
   buf=a;
   buf.Resize(n);
//--- function call
   CFtBase::FtApplyPlan(plan,buf,0,1);
//--- change value
   a.Set(0,buf[0]+buf[1]);
   for(i=1; i<=n2-1; i++)
     {
      //--- calculation
      idx=2*(i%n2);
      hn.real=buf[idx+0];
      hn.imag=buf[idx+1];
      idx=2*(n2-i);
      hmnc.real=buf[idx+0];
      hmnc.imag=-buf[idx+1];
      v.real=-MathSin(-(2*M_PI*i/n));
      v.imag=MathCos(-(2*M_PI*i/n));
      v=hn+hmnc-v*(hn-hmnc);
      a.Set(2*i+0,0.5*v.real);
      a.Set(2*i+1,0.5*v.imag);
     }
//--- change value
   a.Set(1,buf[0]-buf[1]);
  }
//+------------------------------------------------------------------+
//| Internal subroutine. Never call it directly!                     |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInvInternalEven(double &a[],const int n,
                                                  double &buf[],CFtPlan &plan)
  {
//--- create variables
   double x=0;
   double y=0;
   double t=0;
   int    i=0;
   int    n2=0;
//--- check
   if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
      return;
//--- Special cases:
//--- * N=2
//--- After this block we assume that N is strictly greater than 2
   if(n==2)
     {
      //--- change values
      x=0.5*(a[0]+a[1]);
      y=0.5*(a[0]-a[1]);
      a[0]=x;
      a[1]=y;
      //--- exit the function
      return;
     }
//--- inverse real FFT is reduced to the inverse real FHT,
//--- which is reduced to the forward real FHT,
//--- which is reduced to the forward real FFT.
//--- Don't worry,it is really compact and efficient reduction :)
   n2=n/2;
   buf[0]=a[0];
//--- calculation
   for(i=1; i<=n2-1; i++)
     {
      x=a[2*i+0];
      y=a[2*i+1];
      buf[i]=x-y;
      buf[n-i]=x+y;
     }
   buf[n2]=a[1];
//--- function call
   FFTR1DInternalEven(buf,n,a,plan);
//--- change values
   a[0]=buf[0]/n;
   t=1.0/(double)n;
//--- calculation
   for(i=1; i<=n2-1; i++)
     {
      x=buf[2*i+0];
      y=buf[2*i+1];
      a[i]=t*(x-y);
      a[n-i]=t*(x+y);
     }
//--- change value
   a[n2]=buf[1]/n;
  }
//+------------------------------------------------------------------+
//| Internal subroutine. Never call it directly!                     |
//+------------------------------------------------------------------+
void CFastFourierTransform::FFTR1DInvInternalEven(CRowDouble &a,const int n,
                                                  CRowDouble &buf,CFtPlan &plan)
  {
//--- create variables
   double x=0;
   double y=0;
   double t=0;
   int    i=0;
   int    n2=0;
//--- check
   if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
      return;
//--- Special cases:
//--- * N=2
//--- After this block we assume that N is strictly greater than 2
   if(n==2)
     {
      //--- change values
      x=0.5*(a[0]+a[1]);
      y=0.5*(a[0]-a[1]);
      a.Set(0,x);
      a.Set(1,y);
      //--- exit the function
      return;
     }
//--- inverse real FFT is reduced to the inverse real FHT,
//--- which is reduced to the forward real FHT,
//--- which is reduced to the forward real FFT.
//--- Don't worry,it is really compact and efficient reduction :)
   n2=n/2;
   buf.Set(0,a[0]);
//--- calculation
   for(i=1; i<=n2-1; i++)
     {
      x=a[2*i+0];
      y=a[2*i+1];
      buf.Set(i,x-y);
      buf.Set(n-i,x+y);
     }
   buf.Set(n2,a[1]);
//--- function call
   FFTR1DInternalEven(buf,n,a,plan);
//--- change values
   a.Set(0,buf[0]/n);
   t=1.0/(double)n;
//--- calculation
   for(i=1; i<=n2-1; i++)
     {
      x=buf[2*i+0];
      y=buf[2*i+1];
      a.Set(i,t*(x-y));
      a.Set(n-i,t*(x+y));
     }
//--- change value
   a.Set(n2,buf[1]/n);
  }
//+------------------------------------------------------------------+
//| Convolution class                                                |
//+------------------------------------------------------------------+
class CConv
  {
public:
   static void       ConvC1D(complex &a[],const int m,complex &b[],const int n,complex &r[]);
   static void       ConvC1DInv(complex &a[],const int m,complex &b[],const int n,complex &r[]);
   static void       ConvC1DCircular(complex &s[],const int m,complex &r[],const int n,complex &c[]);
   static void       ConvC1DCircularInv(complex &a[],const int m,complex &b[],const int n,complex &r[]);
   static void       ConvR1D(double &a[],const int m,double &b[],const int n,double &r[]);
   static void       ConvR1DInv(double &a[],const int m,double &b[],const int n,double &r[]);
   static void       ConvR1DCircular(double &s[],const int m,double &r[],const int n,double &c[]);
   static void       ConvR1DCircularInv(double &a[],const int m,double &b[],const int n,double &r[]);
   static void       ConvC1DX(complex &a[],const int m,complex &b[],const int n,const bool circular,int alg,int q,complex &r[]);
   static void       ConvR1DX(double &a[],const int m,double &b[],const int n,const bool circular,int alg,int q,double &r[]);
  };
//+------------------------------------------------------------------+
//| 1-dimensional complex convolution.                               |
//| For given A/B returns conv(A,B) (non-circular). Subroutine can   |
//| automatically choose between three implementations:              |
//| straightforward O(M*N) formula for very small N (or M),          |
//| significantly larger than min(M,N), but O(M*N) algorithm is too  |
//| slow, and general FFT-based formula for cases where two previois |
//| algorithms are too slow.                                         |
//| Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.     |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - complex function to be transformed   |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - complex function to be transformed   |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..N+M-2].                   |
//| NOTE:                                                            |
//|     It is assumed that A is zero at T<0, B is zero too. If one or|
//| both functions have non-zero values at negative T's, you can     |
//| still use this subroutine - just shift its result                |
//| correspondingly.                                                 |
//+------------------------------------------------------------------+
void CConv::ConvC1D(complex &a[],const int m,complex &b[],const int n,complex &r[])
  {
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer that B.
   if(m<n)
     {
      //--- function call
      ConvC1D(b,n,a,m,r);
      return;
     }
//--- function call
   ConvC1DX(a,m,b,n,false,-1,0,r);
  }
//+------------------------------------------------------------------+
//| 1-dimensional complex non-circular deconvolution (inverse of     |
//| ConvC1D()).                                                      |
//| Algorithm has M*log(M)) complexity for any M (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - convolved signal, A = conv(R, B)     |
//|     M   -   convolved signal length                              |
//|     B   -   array[0..N-1] - response                             |
//|     N   -   response length, N<=M                                |
//| OUTPUT PARAMETERS                                                |
//|     R   -   deconvolved signal. array[0..M-N].                   |
//| NOTE:                                                            |
//|     deconvolution is unstable process and may result in division |
//| by zero (if your response function is degenerate, i.e. has zero  |
//| Fourier coefficient).                                            |
//| NOTE:                                                            |
//|     It is assumed that A is zero at T<0, B is zero too. If one   |
//| or both functions have non-zero values at negative T's, you can  |
//| still use this subroutine - just shift its result correspondingly|
//+------------------------------------------------------------------+
void CConv::ConvC1DInv(complex &a[],const int m,complex &b[],
                       const int n,complex &r[])
  {
//--- create variables
   int     i=0;
   int     p=0;
   complex c1=0;
   complex c2=0;
   complex c3=0;
   double  t=0;
//--- create arrays
   double buf[];
   double buf2[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert((n>0 && m>0) && n<=m,__FUNCTION__+": incorrect N or M!"))
      return;
//--- function call
   p=CFtBase::FtBaseFindSmooth(m);
   CFtBase::FtComplexFFTPlan(p,1,plan);
//--- allocation
   ArrayResize(buf,2*p);
//--- copy
   for(i=0; i<=m-1; i++)
     {
      buf[2*i+0]=a[i].real;
      buf[2*i+1]=a[i].imag;
     }
//--- make zero
   for(i=m; i<=p-1; i++)
     {
      buf[2*i+0]=0;
      buf[2*i+1]=0;
     }
//--- allocation
   ArrayResize(buf2,2*p);
//--- copy
   for(i=0; i<=n-1; i++)
     {
      buf2[2*i+0]=b[i].real;
      buf2[2*i+1]=b[i].imag;
     }
//--- make zero
   for(i=n; i<=p-1; i++)
     {
      buf2[2*i+0]=0;
      buf2[2*i+1]=0;
     }
//--- function call
   CFtBase::FtApplyPlan(plan,buf,0,1);
   CFtBase::FtApplyPlan(plan,buf2,0,1);
//--- calculation
   for(i=0; i<=p-1; i++)
     {
      c1.real=buf[2*i+0];
      c1.imag=buf[2*i+1];
      c2.real=buf2[2*i+0];
      c2.imag=buf2[2*i+1];
      c3=c1/c2;
      buf[2*i+0]=c3.real;
      buf[2*i+1]=-c3.imag;
     }
//--- function call
   CFtBase::FtApplyPlan(plan,buf,0,1);
   t=1.0/(double)p;
//--- allocation
   ArrayResize(r,m-n+1);
//--- change values
   for(i=0; i<=m-n; i++)
     {
      r[i].real=t*buf[2*i+0];
      r[i].imag=-(t*buf[2*i+1]);
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional circular complex convolution.                      |
//| For given S/R returns conv(S,R) (circular). Algorithm has        |
//| linearithmic complexity for any M/N.                             |
//| IMPORTANT:  normal convolution is commutative, i.e. it is        |
//| symmetric - conv(A,B)=conv(B,A). Cyclic convolution IS NOT. One  |
//| function - S - is a signal, periodic function, and another - R - |
//| is a response, non-periodic function with limited length.        |
//| INPUT PARAMETERS                                                 |
//|     S   -   array[0..M-1] - complex periodic signal              |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - complex non-periodic response        |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..M-1].                     |
//| NOTE:                                                            |
//|     It is assumed that B is zero at T<0. If it has non-zero      |
//| values at negative T's, you can still use this subroutine - just |
//| shift its result correspondingly.                                |
//+------------------------------------------------------------------+
void CConv::ConvC1DCircular(complex &s[],const int m,complex &r[],
                            const int n,complex &c[])
  {
//--- create variables
   int i1=0;
   int i2=0;
   int j2=0;
   int i_=0;
   int i1_=0;
//--- create array
   complex buf[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(buf,m);
      //--- initialization
      for(i1=0; i1<=m-1; i1++)
         buf[i1]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         //--- change values
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            buf[i_]=buf[i_]+r[i_+i1_];
         i1=i1+m;
        }
      //--- function call
      ConvC1DCircular(s,m,buf,m,c);
      //--- exit the function
      return;
     }
//--- function call
   ConvC1DX(s,m,r,n,true,-1,0,c);
  }
//+------------------------------------------------------------------+
//| 1-dimensional circular complex deconvolution (inverse of         |
//| ConvC1DCircular()).                                              |
//| Algorithm has M*log(M)) complexity for any M (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - convolved periodic signal,           |
//|             A = conv(R, B)                                       |
//|     M   -   convolved signal length                              |
//|     B   -   array[0..N-1] - non-periodic response                |
//|     N   -   response length                                      |
//| OUTPUT PARAMETERS                                                |
//|     R   -   deconvolved signal. array[0..M-1].                   |
//| NOTE:                                                            |
//|     deconvolution is unstable process and may result in division |
//| by zero (if your response function is degenerate, i.e. has zero  |
//| Fourier coefficient).                                            |
//| NOTE:                                                            |
//|     It is assumed that B is zero at T<0. If it has non-zero      |
//| values at negative T's, you can still use this subroutine - just |
//| shift its result correspondingly.                                |
//+------------------------------------------------------------------+
void CConv::ConvC1DCircularInv(complex &a[],const int m,complex &b[],
                               const int n,complex &r[])
  {
//--- create variables
   int     i=0;
   int     i1=0;
   int     i2=0;
   int     j2=0;
   complex c1=0;
   complex c2=0;
   complex c3=0;
   double  t=0;
   int     i_=0;
   int     i1_=0;
//--- create arrays
   double  buf[];
   double  buf2[];
   complex cbuf[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(cbuf,m);
      //--- initialization
      for(i=0; i<=m-1; i++)
         cbuf[i]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         //--- change values
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            cbuf[i_]=cbuf[i_]+b[i_+i1_];
         i1=i1+m;
        }
      //--- function call
      ConvC1DCircularInv(a,m,cbuf,m,r);
      //--- exit the function
      return;
     }
//--- Task is normalized
   CFtBase::FtComplexFFTPlan(m,1,plan);
//--- allocation
   ArrayResize(buf,2*m);
//--- copy
   for(i=0; i<=m-1; i++)
     {
      buf[2*i+0]=a[i].real;
      buf[2*i+1]=a[i].imag;
     }
//--- allocation
   ArrayResize(buf2,2*m);
//--- copy
   for(i=0; i<=n-1; i++)
     {
      buf2[2*i+0]=b[i].real;
      buf2[2*i+1]=b[i].imag;
     }
//--- make zero
   for(i=n; i<=m-1; i++)
     {
      buf2[2*i+0]=0;
      buf2[2*i+1]=0;
     }
//--- function call
   CFtBase::FtApplyPlan(plan,buf,0,1);
   CFtBase::FtApplyPlan(plan,buf2,0,1);
//--- calculation
   for(i=0; i<=m-1; i++)
     {
      c1.real=buf[2*i+0];
      c1.imag=buf[2*i+1];
      c2.real=buf2[2*i+0];
      c2.imag=buf2[2*i+1];
      c3=c1/c2;
      buf[2*i+0]=c3.real;
      buf[2*i+1]=-c3.imag;
     }
//--- function call
   CFtBase::FtApplyPlan(plan,buf,0,1);
   t=1.0/(double)m;
//--- allocation
   ArrayResize(r,m);
//--- change values
   for(i=0; i<=m-1; i++)
     {
      r[i].real=t*buf[2*i+0];
      r[i].imag=-(t*buf[2*i+1]);
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional real convolution.                                  |
//| Analogous to ConvC1D(), see ConvC1D() comments for more details. |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - real function to be transformed      |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - real function to be transformed      |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..N+M-2].                   |
//| NOTE:                                                            |
//|     It is assumed that A is zero at T<0, B is zero too. If one   |
//| or both functions have non-zero values at negative T's, you can  |
//| still use this subroutine - just shift its result correspondingly|
//+------------------------------------------------------------------+
void CConv::ConvR1D(double &a[],const int m,double &b[],
                    const int n,double &r[])
  {
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer that B.
   if(m<n)
     {
      //--- function call
      ConvR1D(b,n,a,m,r);
      return;
     }
//--- function call
   ConvR1DX(a,m,b,n,false,-1,0,r);
  }
//+------------------------------------------------------------------+
//| 1-dimensional real deconvolution (inverse of ConvC1D()).         |
//| Algorithm has M*log(M)) complexity for any M (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - convolved signal, A = conv(R, B)     |
//|     M   -   convolved signal length                              |
//|     B   -   array[0..N-1] - response                             |
//|     N   -   response length, N<=M                                |
//| OUTPUT PARAMETERS                                                |
//|     R   -   deconvolved signal. array[0..M-N].                   |
//| NOTE:                                                            |
//|     deconvolution is unstable process and may result in division |
//| by zero (if your response function is degenerate, i.e. has zero  |
//| Fourier coefficient).                                            |
//| NOTE:                                                            |
//|     It is assumed that A is zero at T<0, B is zero too. If one or|
//| both functions have non-zero values at negative T's, you can     |
//| still use this subroutine - just shift its result correspondingly|
//+------------------------------------------------------------------+
void CConv::ConvR1DInv(double &a[],const int m,double &b[],
                       const int n,double &r[])
  {
//--- create variables
   int     i=0;
   int     p=0;
   complex c1=0;
   complex c2=0;
   complex c3=0;
   int     i_=0;
//--- create arrays
   double buf[];
   double buf2[];
   double buf3[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert((n>0 && m>0) && n<=m,__FUNCTION__+": incorrect N or M!"))
      return;
//--- function call
   p=CFtBase::FtBaseFindSmoothEven(m);
//--- allocation
   ArrayResize(buf,p);
//--- copy
   for(i_=0; i_<=m-1; i_++)
      buf[i_]=a[i_];
//--- make zero
   for(i=m; i<=p-1; i++)
      buf[i]=0;
//--- allocation
   ArrayResize(buf2,p);
   for(i_=0; i_<=n-1; i_++)
      buf2[i_]=b[i_];
//--- make zero
   for(i=n; i<=p-1; i++)
      buf2[i]=0;
//--- allocation
   ArrayResize(buf3,p);
//--- function call
   CFtBase::FtComplexFFTPlan(p/2,1,plan);
//--- function call
   CFastFourierTransform::FFTR1DInternalEven(buf,p,buf3,plan);
//--- function call
   CFastFourierTransform::FFTR1DInternalEven(buf2,p,buf3,plan);
//--- change values
   buf[0]=buf[0]/buf2[0];
   buf[1]=buf[1]/buf2[1];
//--- calculation
   for(i=1; i<=p/2-1; i++)
     {
      c1.real=buf[2*i+0];
      c1.imag=buf[2*i+1];
      c2.real=buf2[2*i+0];
      c2.imag=buf2[2*i+1];
      c3=c1/c2;
      buf[2*i+0]=c3.real;
      buf[2*i+1]=c3.imag;
     }
//--- function call
   CFastFourierTransform::FFTR1DInvInternalEven(buf,p,buf3,plan);
//--- allocation
   ArrayResize(r,m-n+1);
//--- copy
   for(i_=0; i_<=m-n; i_++)
      r[i_]=buf[i_];
  }
//+------------------------------------------------------------------+
//| 1-dimensional circular real convolution.                         |
//| Analogous to ConvC1DCircular(), see ConvC1DCircular() comments   |
//| for more details.                                                |
//| INPUT PARAMETERS                                                 |
//|     S   -   array[0..M-1] - real signal                          |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - real response                        |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..M-1].                     |
//| NOTE:                                                            |
//|     It is assumed that B is zero at T<0. If it has non-zero      |
//| values at negative T's, you can still use this subroutine - just |
//| shift its result correspondingly.                                |
//+------------------------------------------------------------------+
void CConv::ConvR1DCircular(double &s[],const int m,double &r[],
                            const int n,double &c[])
  {
//--- create variables
   int i1=0;
   int i2=0;
   int j2=0;
   int i_=0;
   int i1_=0;
//--- create array
   double buf[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(buf,m);
      //--- initialization
      for(i1=0; i1<=m-1; i1++)
         buf[i1]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            buf[i_]=buf[i_]+r[i_+i1_];
         i1=i1+m;
        }
      ConvR1DCircular(s,m,buf,m,c);
      //--- exit the function
      return;
     }
//--- reduce to usual convolution
   ConvR1DX(s,m,r,n,true,-1,0,c);
  }
//+------------------------------------------------------------------+
//| 1-dimensional complex deconvolution (inverse of ConvC1D()).      |
//| Algorithm has M*log(M)) complexity for any M (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - convolved signal, A = conv(R, B)     |
//|     M   -   convolved signal length                              |
//|     B   -   array[0..N-1] - response                             |
//|     N   -   response length                                      |
//| OUTPUT PARAMETERS                                                |
//|     R   -   deconvolved signal. array[0..M-N].                   |
//| NOTE:                                                            |
//|     deconvolution is unstable process and may result in division |
//| by zero (if your response function is degenerate, i.e. has zero  |
//| Fourier coefficient).                                            |
//| NOTE:                                                            |
//|     It is assumed that B is zero at T<0. If it has non-zero      |
//| values at negative T's, you can still use this subroutine - just |
//| shift its result correspondingly.                                |
//+------------------------------------------------------------------+
void CConv::ConvR1DCircularInv(double &a[],const int m,double &b[],
                               const int n,double &r[])
  {
//--- create variables
   int     i=0;
   int     i1=0;
   int     i2=0;
   int     j2=0;
   complex c1=0;
   complex c2=0;
   complex c3=0;
   int     i_=0;
   int     i1_=0;
//--- create arrays
   double  buf[];
   double  buf2[];
   double  buf3[];
   complex cbuf[];
   complex cbuf2[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(buf,m);
      //--- initialization
      for(i=0; i<=m-1; i++)
         buf[i]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            buf[i_]=buf[i_]+b[i_+i1_];
         i1=i1+m;
        }
      //--- function call
      ConvR1DCircularInv(a,m,buf,m,r);
      //--- exit the function
      return;
     }
//--- Task is normalized
   if(m%2==0)
     {
      //--- size is even,use fast even-size FFT
      ArrayResize(buf,m);
      for(i_=0; i_<=m-1; i_++)
         buf[i_]=a[i_];
      //--- allocation
      ArrayResize(buf2,m);
      //--- copy
      for(i_=0; i_<=n-1; i_++)
         buf2[i_]=b[i_];
      //--- make zero
      for(i=n; i<=m-1; i++)
         buf2[i]=0;
      //--- allocation
      ArrayResize(buf3,m);
      //--- function call
      CFtBase::FtComplexFFTPlan(m/2,1,plan);
      //--- function call
      CFastFourierTransform::FFTR1DInternalEven(buf,m,buf3,plan);
      //--- function call
      CFastFourierTransform::FFTR1DInternalEven(buf2,m,buf3,plan);
      //--- change values
      buf[0]=buf[0]/buf2[0];
      buf[1]=buf[1]/buf2[1];
      //--- calculation
      for(i=1; i<=m/2-1; i++)
        {
         c1.real=buf[2*i+0];
         c1.imag=buf[2*i+1];
         c2.real=buf2[2*i+0];
         c2.imag=buf2[2*i+1];
         c3=c1/c2;
         buf[2*i+0]=c3.real;
         buf[2*i+1]=c3.imag;
        }
      //--- function call
      CFastFourierTransform::FFTR1DInvInternalEven(buf,m,buf3,plan);
      //--- allocation
      ArrayResize(r,m);
      //--- copy
      for(i_=0; i_<=m-1; i_++)
         r[i_]=buf[i_];
     }
   else
     {
      //--- odd-size,use general real FFT
      CFastFourierTransform::FFTR1D(a,m,cbuf);
      //--- allocation
      ArrayResize(buf2,m);
      //--- copy
      for(i_=0; i_<=n-1; i_++)
         buf2[i_]=b[i_];
      //--- initialization
      for(i=n; i<=m-1; i++)
         buf2[i]=0;
      //--- function call
      CFastFourierTransform::FFTR1D(buf2,m,cbuf2);
      //--- calculation
      for(i=0; i<=(int)MathFloor((double)m/2.0); i++)
         cbuf[i]=cbuf[i]/cbuf2[i];
      //--- function call
      CFastFourierTransform::FFTR1DInv(cbuf,m,r);
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional complex convolution.                               |
//| Extended subroutine which allows to choose convolution algorithm.|
//| Intended for internal use, ALGLIB users should call              |
//| ConvC1D()/ConvC1DCircular().                                     |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - complex function to be transformed   |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - complex function to be transformed   |
//|     N   -   problem size, N<=M                                   |
//|     Alg -   algorithm type:                                      |
//|             *-2     auto-select Q for overlap-add                |
//|             *-1     auto-select algorithm and parameters         |
//|             * 0     straightforward formula for small N's        |
//|             * 1     general FFT-based code                       |
//|             * 2     overlap-add with length Q                    |
//|     Q   -   length for overlap-add                               |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..N+M-1].                   |
//+------------------------------------------------------------------+
void CConv::ConvC1DX(complex &a[],const int m,complex &b[],const int n,
                     const bool circular,int alg,int q,complex &r[])
  {
//--- create variables
   int     i=0;
   int     j=0;
   int     p=0;
   int     ptotal=0;
   int     i1=0;
   int     i2=0;
   int     j1=0;
   int     j2=0;
   complex v=0;
   double  ax=0;
   double  ay=0;
   double  bx=0;
   double  by=0;
   double  t=0;
   double  tx=0;
   double  ty=0;
   double  flopcand=0;
   double  flopbest=0;
   int     algbest=0;
   int     i_=0;
   int     i1_=0;
//--- create arrays
   complex bbuf[];
   double  buf[];
   double  buf2[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- check
   if(!CAp::Assert(n<=m,__FUNCTION__+": N<M assumption is false!"))
      return;
//--- Auto-select
   if(alg==-1 || alg==-2)
     {
      //--- Initial candidate: straightforward implementation.
      //--- If we want to use auto-fitted overlap-add,
      //--- flop count is initialized by large real number - to force
      //--- another algorithm selection
      algbest=0;
      //--- check
      if(alg==-1)
         flopbest=2*m*n;
      else
         flopbest=CMath::m_maxrealnumber;
      //--- Another candidate - generic FFT code
      if(alg==-1)
        {
         //--- check
         if(circular && CFtBase::FtBaseIsSmooth(m))
           {
            //--- special code for circular convolution of a sequence with a smooth length
            flopcand=3*CFtBase::FtBaseGetFlopEstimate(m)+6*m;
            //--- check
            if(flopcand<flopbest)
              {
               algbest=1;
               flopbest=flopcand;
              }
           }
         else
           {
            //--- general cyclic/non-cyclic convolution
            p=CFtBase::FtBaseFindSmooth(m+n-1);
            flopcand=3*CFtBase::FtBaseGetFlopEstimate(p)+6*p;
            //--- check
            if(flopcand<flopbest)
              {
               algbest=1;
               flopbest=flopcand;
              }
           }
        }
      //--- Another candidate - overlap-add
      q=1;
      ptotal=1;
      while(ptotal<n)
         ptotal=ptotal*2;
      //--- calculation
      while(ptotal<=m+n-1)
        {
         //--- change values
         p=ptotal-n+1;
         flopcand=(int)MathCeil((double)m/(double)p)*(2*CFtBase::FtBaseGetFlopEstimate(ptotal)+8*ptotal);
         //--- check
         if(flopcand<flopbest)
           {
            flopbest=flopcand;
            algbest=2;
            q=p;
           }
         ptotal=ptotal*2;
        }
      //--- change value
      alg=algbest;
      //--- function call
      ConvC1DX(a,m,b,n,circular,alg,q,r);
      //--- exit the function
      return;
     }
//--- straightforward formula for
//--- circular and non-circular convolutions.
//--- Very simple code,no further comments needed.
   if(alg==0)
     {
      //--- Special case: N=1
      if(n==1)
        {
         //--- allocation
         ArrayResize(r,m);
         v=b[0];
         for(i_=0; i_<=m-1; i_++)
            r[i_]=v*a[i_];
         //--- exit the function
         return;
        }
      //--- use straightforward formula
      if(circular)
        {
         //--- circular convolution
         ArrayResize(r,m);
         //--- change values
         v=b[0];
         for(i_=0; i_<=m-1; i_++)
            r[i_]=v*a[i_];
         //--- calculation
         for(i=1; i<=n-1; i++)
           {
            //--- change values
            v=b[i];
            i1=0;
            i2=i-1;
            j1=m-i;
            j2=m-1;
            i1_=j1-i1;
            //--- calculation
            for(i_=i1; i_<=i2; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
            //--- change values
            i1=i;
            i2=m-1;
            j1=0;
            j2=m-i-1;
            i1_=j1-i1;
            //--- calculation
            for(i_=i1; i_<=i2; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
           }
        }
      else
        {
         //--- non-circular convolution
         ArrayResize(r,m+n-1);
         //--- initialization
         for(i=0; i<=m+n-2; i++)
            r[i]=0;
         //--- calculation
         for(i=0; i<=n-1; i++)
           {
            v=b[i];
            i1_=-i;
            for(i_=i; i_<=i+m-1; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
           }
        }
      //--- exit the function
      return;
     }
//--- general FFT-based code for
//--- circular and non-circular convolutions.
//--- First,if convolution is circular,we test whether M is smooth or not.
//--- If it is smooth,we just use M-length FFT to calculate convolution.
//--- If it is not,we calculate non-circular convolution and wrap it arount.
//--- IF convolution is non-circular,we use zero-padding + FFT.
   if(alg==1)
     {
      //--- check
      if(circular && CFtBase::FtBaseIsSmooth(m))
        {
         //--- special code for circular convolution with smooth M
         CFtBase::FtComplexFFTPlan(m,1,plan);
         //--- allocation
         ArrayResize(buf,2*m);
         //--- copy
         for(i=0; i<=m-1; i++)
           {
            buf[2*i+0]=a[i].real;
            buf[2*i+1]=a[i].imag;
           }
         //--- allocation
         ArrayResize(buf2,2*m);
         //--- copy
         for(i=0; i<=n-1; i++)
           {
            buf2[2*i+0]=b[i].real;
            buf2[2*i+1]=b[i].imag;
           }
         //--- make zero
         for(i=n; i<=m-1; i++)
           {
            buf2[2*i+0]=0;
            buf2[2*i+1]=0;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         CFtBase::FtApplyPlan(plan,buf2,0,1);
         //--- calculation
         for(i=0; i<=m-1; i++)
           {
            ax=buf[2*i+0];
            ay=buf[2*i+1];
            bx=buf2[2*i+0];
            by=buf2[2*i+1];
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*i+0]=tx;
            buf[2*i+1]=-ty;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         t=1.0/(double)m;
         //--- allocation
         ArrayResize(r,m);
         //--- change values
         for(i=0; i<=m-1; i++)
           {
            r[i].real=t*buf[2*i+0];
            r[i].imag=-(t*buf[2*i+1]);
           }
        }
      else
        {
         //--- M is non-smooth,general code (circular/non-circular):
         //--- * first part is the same for circular and non-circular
         //---   convolutions. zero padding,FFTs,inverse FFTs
         //--- * second part differs:
         //---   * for non-circular convolution we just copy array
         //---   * for circular convolution we add array tail to its head
         p=CFtBase::FtBaseFindSmooth(m+n-1);
         CFtBase::FtComplexFFTPlan(p,1,plan);
         //--- allocation
         ArrayResize(buf,2*p);
         //--- copy
         for(i=0; i<=m-1; i++)
           {
            buf[2*i+0]=a[i].real;
            buf[2*i+1]=a[i].imag;
           }
         //--- make zero
         for(i=m; i<=p-1; i++)
           {
            buf[2*i+0]=0;
            buf[2*i+1]=0;
           }
         //--- allocation
         ArrayResize(buf2,2*p);
         //--- copy
         for(i=0; i<=n-1; i++)
           {
            buf2[2*i+0]=b[i].real;
            buf2[2*i+1]=b[i].imag;
           }
         //--- make zero
         for(i=n; i<=p-1; i++)
           {
            buf2[2*i+0]=0;
            buf2[2*i+1]=0;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         CFtBase::FtApplyPlan(plan,buf2,0,1);
         //--- calculation
         for(i=0; i<=p-1; i++)
           {
            ax=buf[2*i+0];
            ay=buf[2*i+1];
            bx=buf2[2*i+0];
            by=buf2[2*i+1];
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*i+0]=tx;
            buf[2*i+1]=-ty;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         t=1.0/(double)p;
         //--- check
         if(circular)
           {
            //--- circular,add tail to head
            ArrayResize(r,m);
            //--- copy
            for(i=0; i<=m-1; i++)
              {
               r[i].real=t*buf[2*i+0];
               r[i].imag=-(t*buf[2*i+1]);
              }
            //--- change values
            for(i=m; i<=m+n-2; i++)
              {
               r[i-m].real=r[i-m].real+t*buf[2*i+0];
               r[i-m].imag=r[i-m].imag-t*buf[2*i+1];
              }
           }
         else
           {
            //--- non-circular,just copy
            ArrayResize(r,m+n-1);
            for(i=0; i<=m+n-2; i++)
              {
               r[i].real=t*buf[2*i+0];
               r[i].imag=-(t*buf[2*i+1]);
              }
           }
        }
      //--- exit the function
      return;
     }
//--- overlap-add method for
//--- circular and non-circular convolutions.
//--- First part of code (separate FFTs of input blocks) is the same
//--- for all types of convolution. Second part (overlapping outputs)
//--- differs for different types of convolution. We just copy output
//--- when convolution is non-circular. We wrap it around,if it is
//--- circular.
   if(alg==2)
     {
      //--- allocation
      ArrayResize(buf,2*(q+n-1));
      //--- prepare R
      if(circular)
        {
         //--- allocation
         ArrayResize(r,m);
         for(i=0; i<=m-1; i++)
            r[i]=0;
        }
      else
        {
         //--- allocation
         ArrayResize(r,m+n-1);
         for(i=0; i<=m+n-2; i++)
            r[i]=0;
        }
      //--- pre-calculated FFT(B)
      ArrayResize(bbuf,q+n-1);
      for(i_=0; i_<=n-1; i_++)
         bbuf[i_]=b[i_];
      for(j=n; j<=q+n-2; j++)
         bbuf[j]=0;
      //--- function call
      CFastFourierTransform::FFTC1D(bbuf,q+n-1);
      //--- prepare FFT plan for chunks of A
      CFtBase::FtComplexFFTPlan(q+n-1,1,plan);
      //--- main overlap-add cycle
      i=0;
      //--- calculation
      while(i<=m-1)
        {
         p=MathMin(q,m-i);
         //--- copy
         for(j=0; j<=p-1; j++)
           {
            buf[2*j+0]=a[i+j].real;
            buf[2*j+1]=a[i+j].imag;
           }
         //--- make zero
         for(j=p; j<=q+n-2; j++)
           {
            buf[2*j+0]=0;
            buf[2*j+1]=0;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         //--- calculation
         for(j=0; j<=q+n-2; j++)
           {
            ax=buf[2*j+0];
            ay=buf[2*j+1];
            bx=bbuf[j].real;
            by=bbuf[j].imag;
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*j+0]=tx;
            buf[2*j+1]=-ty;
           }
         //--- function call
         CFtBase::FtApplyPlan(plan,buf,0,1);
         t=1.0/(double)(q+n-1);
         //--- check
         if(circular)
           {
            j1=MathMin(i+p+n-2,m-1)-i;
            j2=j1+1;
           }
         else
           {
            j1=p+n-2;
            j2=j1+1;
           }
         //--- change values
         for(j=0; j<=j1; j++)
           {
            r[i+j].real=r[i+j].real+buf[2*j+0]*t;
            r[i+j].imag=r[i+j].imag-buf[2*j+1]*t;
           }
         //--- change values
         for(j=j2; j<=p+n-2; j++)
           {
            r[j-j2].real=r[j-j2].real+buf[2*j+0]*t;
            r[j-j2].imag=r[j-j2].imag-buf[2*j+1]*t;
           }
         i=i+p;
        }
      //--- exit the function
      return;
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional real convolution.                                  |
//| Extended subroutine which allows to choose convolution algorithm.|
//| Intended for internal use, ALGLIB users should call ConvR1D().   |
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..M-1] - complex function to be transformed   |
//|     M   -   problem size                                         |
//|     B   -   array[0..N-1] - complex function to be transformed   |
//|     N   -   problem size, N<=M                                   |
//|     Alg -   algorithm type:                                      |
//|             *-2     auto-select Q for overlap-add                |
//|             *-1     auto-select algorithm and parameters         |
//|             * 0     straightforward formula for small N's        |
//|             * 1     general FFT-based code                       |
//|             * 2     overlap-add with length Q                    |
//|     Q   -   length for overlap-add                               |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..N+M-1].                   |
//+------------------------------------------------------------------+
void CConv::ConvR1DX(double &a[],const int m,double &b[],const int n,
                     const bool circular,int alg,int q,double &r[])
  {
//--- create variables
   double v=0;
   int    i=0;
   int    j=0;
   int    p=0;
   int    ptotal=0;
   int    i1=0;
   int    i2=0;
   int    j1=0;
   int    j2=0;
   double ax=0;
   double ay=0;
   double bx=0;
   double by=0;
   double tx=0;
   double ty=0;
   double flopcand=0;
   double flopbest=0;
   int    algbest=0;
   int    i_=0;
   int    i1_=0;
//--- create arrays
   double buf[];
   double buf2[];
   double buf3[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- check
   if(!CAp::Assert(n<=m,__FUNCTION__+": N<M assumption is false!"))
      return;
//--- handle special cases
   if(MathMin(m,n)<=2)
      alg=0;
//--- Auto-select
   if(alg<0)
     {
      //--- Initial candidate: straightforward implementation.
      //--- If we want to use auto-fitted overlap-add,
      //--- flop count is initialized by large real number - to force
      //--- another algorithm selection
      algbest=0;
      //--- check
      if(alg==-1)
         flopbest=0.15*m*n;
      else
         flopbest=CMath::m_maxrealnumber;
      //--- Another candidate - generic FFT code
      if(alg==-1)
        {
         //--- check
         if((circular && CFtBase::FtBaseIsSmooth(m)) && m%2==0)
           {
            //--- special code for circular convolution of a sequence with a smooth length
            flopcand=3*CFtBase::FtBaseGetFlopEstimate(m/2)+(double)(6*m)/2.0;
            //--- check
            if(flopcand<flopbest)
              {
               algbest=1;
               flopbest=flopcand;
              }
           }
         else
           {
            //--- general cyclic/non-cyclic convolution
            p=CFtBase::FtBaseFindSmoothEven(m+n-1);
            flopcand=3*CFtBase::FtBaseGetFlopEstimate(p/2)+(double)(6*p)/2.0;
            //--- check
            if(flopcand<flopbest)
              {
               algbest=1;
               flopbest=flopcand;
              }
           }
        }
      //--- Another candidate - overlap-add
      q=1;
      ptotal=1;
      while(ptotal<n)
         ptotal=ptotal*2;
      //--- calculation
      while(ptotal<=m+n-1)
        {
         //--- change values
         p=ptotal-n+1;
         flopcand=(int)MathCeil((double)m/(double)p)*(2*CFtBase::FtBaseGetFlopEstimate(ptotal/2)+1*(ptotal/2));
         //--- check
         if(flopcand<flopbest)
           {
            flopbest=flopcand;
            algbest=2;
            q=p;
           }
         ptotal=ptotal*2;
        }
      alg=algbest;
      //--- function call
      ConvR1DX(a,m,b,n,circular,alg,q,r);
      //--- exit the function
      return;
     }
//--- straightforward formula for
//--- circular and non-circular convolutions.
//--- Very simple code,no further comments needed.
   if(alg==0)
     {
      //--- Special case: N=1
      if(n==1)
        {
         //--- allocation
         ArrayResize(r,m);
         v=b[0];
         for(i_=0; i_<=m-1; i_++)
            r[i_]=v*a[i_];
         //--- exit the function
         return;
        }
      //--- use straightforward formula
      if(circular)
        {
         //--- circular convolution
         ArrayResize(r,m);
         v=b[0];
         for(i_=0; i_<=m-1; i_++)
            r[i_]=v*a[i_];
         //--- calculation
         for(i=1; i<=n-1; i++)
           {
            //--- change values
            v=b[i];
            i1=0;
            i2=i-1;
            j1=m-i;
            j2=m-1;
            i1_=j1-i1;
            //--- calculation
            for(i_=i1; i_<=i2; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
            //--- change values
            i1=i;
            i2=m-1;
            j1=0;
            j2=m-i-1;
            i1_=j1-i1;
            //--- calculation
            for(i_=i1; i_<=i2; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
           }
        }
      else
        {
         //--- non-circular convolution
         ArrayResize(r,m+n-1);
         for(i=0; i<=m+n-2; i++)
            r[i]=0;
         //--- calculation
         for(i=0; i<=n-1; i++)
           {
            v=b[i];
            i1_=-i;
            for(i_=i; i_<=i+m-1; i_++)
               r[i_]=r[i_]+v*a[i_+i1_];
           }
        }
      //--- exit the function
      return;
     }
//--- general FFT-based code for
//--- circular and non-circular convolutions.
//--- First,if convolution is circular,we test whether M is smooth or not.
//--- If it is smooth,we just use M-length FFT to calculate convolution.
//--- If it is not,we calculate non-circular convolution and wrap it arount.
//--- If convolution is non-circular,we use zero-padding + FFT.
//--- We assume that M+N-1>2 - we should call small case code otherwise
   if(alg==1)
     {
      //--- check
      if(!CAp::Assert(m+n-1>2,__FUNCTION__+": internal error!"))
         return;
      //--- check
      if((circular && CFtBase::FtBaseIsSmooth(m)) && m%2==0)
        {
         //--- special code for circular convolution with smooth even M
         ArrayResize(buf,m);
         for(i_=0; i_<=m-1; i_++)
            buf[i_]=a[i_];
         //--- allocation
         ArrayResize(buf2,m);
         for(i_=0; i_<=n-1; i_++)
            buf2[i_]=b[i_];
         for(i=n; i<=m-1; i++)
            buf2[i]=0;
         //--- allocation
         ArrayResize(buf3,m);
         //--- function call
         CFtBase::FtComplexFFTPlan(m/2,1,plan);
         //--- function call
         CFastFourierTransform::FFTR1DInternalEven(buf,m,buf3,plan);
         //--- function call
         CFastFourierTransform::FFTR1DInternalEven(buf2,m,buf3,plan);
         //--- change values
         buf[0]=buf[0]*buf2[0];
         buf[1]=buf[1]*buf2[1];
         //--- calculation
         for(i=1; i<=m/2-1; i++)
           {
            ax=buf[2*i+0];
            ay=buf[2*i+1];
            bx=buf2[2*i+0];
            by=buf2[2*i+1];
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*i+0]=tx;
            buf[2*i+1]=ty;
           }
         //--- function call
         CFastFourierTransform::FFTR1DInvInternalEven(buf,m,buf3,plan);
         //--- allocation
         ArrayResize(r,m);
         //--- copy
         for(i_=0; i_<=m-1; i_++)
            r[i_]=buf[i_];
        }
      else
        {
         //--- M is non-smooth or non-even,general code (circular/non-circular):
         //--- * first part is the same for circular and non-circular
         //---   convolutions. zero padding,FFTs,inverse FFTs
         //--- * second part differs:
         //---   * for non-circular convolution we just copy array
         //---   * for circular convolution we add array tail to its head
         p=CFtBase::FtBaseFindSmoothEven(m+n-1);
         //--- allocation
         ArrayResize(buf,p);
         for(i_=0; i_<=m-1; i_++)
            buf[i_]=a[i_];
         for(i=m; i<=p-1; i++)
            buf[i]=0;
         //--- allocation
         ArrayResize(buf2,p);
         for(i_=0; i_<=n-1; i_++)
            buf2[i_]=b[i_];
         for(i=n; i<=p-1; i++)
            buf2[i]=0;
         //--- allocation
         ArrayResize(buf3,p);
         //--- function call
         CFtBase::FtComplexFFTPlan(p/2,1,plan);
         //--- function call
         CFastFourierTransform::FFTR1DInternalEven(buf,p,buf3,plan);
         //--- function call
         CFastFourierTransform::FFTR1DInternalEven(buf2,p,buf3,plan);
         //--- change values
         buf[0]=buf[0]*buf2[0];
         buf[1]=buf[1]*buf2[1];
         //--- calculation
         for(i=1; i<=p/2-1; i++)
           {
            ax=buf[2*i+0];
            ay=buf[2*i+1];
            bx=buf2[2*i+0];
            by=buf2[2*i+1];
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*i+0]=tx;
            buf[2*i+1]=ty;
           }
         //--- function call
         CFastFourierTransform::FFTR1DInvInternalEven(buf,p,buf3,plan);
         //--- check
         if(circular)
           {
            //--- circular,add tail to head
            ArrayResize(r,m);
            for(i_=0; i_<=m-1; i_++)
               r[i_]=buf[i_];
            //--- check
            if(n>=2)
              {
               i1_=m;
               for(i_=0; i_<=n-2; i_++)
                  r[i_]=r[i_]+buf[i_+i1_];
              }
           }
         else
           {
            //--- non-circular,just copy
            ArrayResize(r,m+n-1);
            for(i_=0; i_<=m+n-2; i_++)
               r[i_]=buf[i_];
           }
        }
      //--- exit the function
      return;
     }
//--- overlap-add method
   if(alg==2)
     {
      //--- check
      if(!CAp::Assert((q+n-1)%2==0,__FUNCTION__+": internal error!"))
         return;
      //--- allocation
      ArrayResize(buf,q+n-1);
      ArrayResize(buf2,q+n-1);
      ArrayResize(buf3,q+n-1);
      //--- function call
      CFtBase::FtComplexFFTPlan((q+n-1)/2,1,plan);
      //--- prepare R
      if(circular)
        {
         //--- allocation
         ArrayResize(r,m);
         for(i=0; i<=m-1; i++)
            r[i]=0;
        }
      else
        {
         //--- allocation
         ArrayResize(r,m+n-1);
         for(i=0; i<=m+n-2; i++)
            r[i]=0;
        }
      //--- pre-calculated FFT(B)
      for(i_=0; i_<=n-1; i_++)
         buf2[i_]=b[i_];
      for(j=n; j<=q+n-2; j++)
         buf2[j]=0;
      //--- function call
      CFastFourierTransform::FFTR1DInternalEven(buf2,q+n-1,buf3,plan);
      //--- main overlap-add cycle
      i=0;
      //--- calculation
      while(i<=m-1)
        {
         p=MathMin(q,m-i);
         i1_=i;
         //--- copy
         for(i_=0; i_<=p-1; i_++)
            buf[i_]=a[i_+i1_];
         //--- make zero
         for(j=p; j<=q+n-2; j++)
            buf[j]=0;
         //--- function call
         CFastFourierTransform::FFTR1DInternalEven(buf,q+n-1,buf3,plan);
         //--- change values
         buf[0]=buf[0]*buf2[0];
         buf[1]=buf[1]*buf2[1];
         //--- calculation
         for(j=1; j<=(q+n-1)/2-1; j++)
           {
            ax=buf[2*j+0];
            ay=buf[2*j+1];
            bx=buf2[2*j+0];
            by=buf2[2*j+1];
            tx=ax*bx-ay*by;
            ty=ax*by+ay*bx;
            buf[2*j+0]=tx;
            buf[2*j+1]=ty;
           }
         //--- function call
         CFastFourierTransform::FFTR1DInvInternalEven(buf,q+n-1,buf3,plan);
         //--- check
         if(circular)
           {
            j1=MathMin(i+p+n-2,m-1)-i;
            j2=j1+1;
           }
         else
           {
            j1=p+n-2;
            j2=j1+1;
           }
         //--- change values
         i1_=-i;
         for(i_=i; i_<=i+j1; i_++)
            r[i_]=r[i_]+buf[i_+i1_];
         //--- check
         if(p+n-2>=j2)
           {
            i1_=j2;
            for(i_=0; i_<=p+n-2-j2; i_++)
               r[i_]=r[i_]+buf[i_+i1_];
           }
         i=i+p;
        }
      //--- exit the function
      return;
     }
  }
//+------------------------------------------------------------------+
//| Cross-correlation class                                          |
//+------------------------------------------------------------------+
class CCorr
  {
public:
   static void       CorrC1D(complex &signal[],const int n,complex &pattern[],const int m,complex &r[]);
   static void       CorrC1DCircular(complex &signal[],const int m,complex &pattern[],const int n,complex &c[]);
   static void       CorrR1D(double &signal[],const int n,double &pattern[],const int m,double &r[]);
   static void       CorrR1DCircular(double &signal[],const int m,double &pattern[],const int n,double &c[]);
  };
//+------------------------------------------------------------------+
//| 1-dimensional complex cross-correlation.                         |
//| For given Pattern/Signal returns corr(Pattern,Signal)            |
//| (non-circular).                                                  |
//| Correlation is calculated using reduction to convolution.        |
//| Algorithm with max(N,N)*log(max(N,N)) complexity is used (see    |
//| ConvC1D() for more info about performance).                      |
//| IMPORTANT:                                                       |
//|     for historical reasons subroutine accepts its parameters in  |
//|     reversed order: CorrC1D(Signal, Pattern) = Pattern x Signal  |
//|     (using traditional definition of cross-correlation, denoting |
//|     cross-correlation as "x").                                   |
//| INPUT PARAMETERS                                                 |
//|     Signal  -   array[0..N-1] - complex function to be           |
//|                 transformed, signal containing pattern           |
//|     N       -   problem size                                     |
//|     Pattern -   array[0..M-1] - complex function to be           |
//|                 transformed, pattern to search withing signal    |
//|     M       -   problem size                                     |
//| OUTPUT PARAMETERS                                                |
//|     R       -   cross-correlation, array[0..N+M-2]:              |
//|                 * positive lags are stored in R[0..N-1],         |
//|                   R[i] = sum(conj(pattern[j])*signal[i+j]        |
//|                 * negative lags are stored in R[N..N+M-2],       |
//|                   R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j] |
//| NOTE:                                                            |
//|     It is assumed that pattern domain is [0..M-1]. If Pattern is |
//| non-zero on [-K..M-1], you can still use this subroutine, just   |
//| shift result by K.                                               |
//+------------------------------------------------------------------+
void CCorr::CorrC1D(complex &signal[],const int n,complex &pattern[],
                    const int m,complex &r[])
  {
//--- create variables
   int i=0;
   int i_=0;
   int i1_=0;
//--- create arrays
   complex p[];
   complex b[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- allocation
   ArrayResize(p,m);
   for(i=0; i<=m-1; i++)
      p[m-1-i]=CMath::Conj(pattern[i]);
//--- function call
   CConv::ConvC1D(p,m,signal,n,b);
//--- allocation
   ArrayResize(r,m+n-1);
   i1_=m-1;
   for(i_=0; i_<=n-1; i_++)
      r[i_]=b[i_+i1_];
//--- check
   if(m+n-2>=n)
     {
      i1_=-n;
      for(i_=n; i_<=m+n-2; i_++)
         r[i_]=b[i_+i1_];
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional circular complex cross-correlation.                |
//| For given Pattern/Signal returns corr(Pattern,Signal) (circular).|
//| Algorithm has linearithmic complexity for any M/N.               |
//| IMPORTANT:                                                       |
//|     for historical reasons subroutine accepts its parameters in  |
//|     reversed order: CorrC1DCircular(Signal, Pattern) = Pattern x |
//|     Signal (using traditional definition of cross-correlation,   |
//|     denoting cross-correlation as "x").                          |
//| INPUT PARAMETERS                                                 |
//|     Signal  -   array[0..N-1] - complex function to be           |
//|                 transformed, periodic signal containing pattern  |
//|     N       -   problem size                                     |
//|     Pattern -   array[0..M-1] - complex function to be           |
//|                 transformed, non-periodic pattern to search      |
//|                 withing signal                                   |
//|     M       -   problem size                                     |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..M-1].                     |
//+------------------------------------------------------------------+
void CCorr::CorrC1DCircular(complex &signal[],const int m,complex &pattern[],
                            const int n,complex &c[])
  {
//--- create variables
   int i1=0;
   int i2=0;
   int i=0;
   int j2=0;
   int i_=0;
   int i1_=0;
//--- create arrays
   complex p[];
   complex b[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(b,m);
      //--- initialization
      for(i1=0; i1<=m-1; i1++)
         b[i1]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         //--- change values
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            b[i_]=b[i_]+pattern[i_+i1_];
         i1=i1+m;
        }
      //--- function call
      CorrC1DCircular(signal,m,b,m,c);
      //--- exit the function
      return;
     }
//--- Task is normalized
   ArrayResize(p,n);
   for(i=0; i<=n-1; i++)
      p[n-1-i]=CMath::Conj(pattern[i]);
//--- function call
   CConv::ConvC1DCircular(signal,m,p,n,b);
//--- allocation
   ArrayResize(c,m);
   i1_=n-1;
//--- copy
   for(i_=0; i_<=m-n; i_++)
      c[i_]=b[i_+i1_];
//--- check
   if(m-n+1<=m-1)
     {
      i1_=-(m-n+1);
      for(i_=m-n+1; i_<=m-1; i_++)
         c[i_]=b[i_+i1_];
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional real cross-correlation.                            |
//| For given Pattern/Signal returns corr(Pattern,Signal)            |
//| (non-circular).                                                  |
//| Correlation is calculated using reduction to convolution.        |
//| Algorithm with max(N,N)*log(max(N,N)) complexity is used (see    |
//| ConvC1D() for more info about performance).                      |
//| IMPORTANT:                                                       |
//|     for  historical reasons subroutine accepts its parameters in |
//|     reversed order: CorrR1D(Signal, Pattern) = Pattern x Signal  |
//|     (using  traditional definition of cross-correlation, denoting|
//|     cross-correlation as "x").                                   |
//| INPUT PARAMETERS                                                 |
//|     Signal  -   array[0..N-1] - real function to be transformed, |
//|                 signal containing pattern                        |
//|     N       -   problem size                                     |
//|     Pattern -   array[0..M-1] - real function to be transformed, |
//|                 pattern to search withing signal                 |
//|     M       -   problem size                                     |
//| OUTPUT PARAMETERS                                                |
//|     R       -   cross-correlation, array[0..N+M-2]:              |
//|                 * positive lags are stored in R[0..N-1],         |
//|                   R[i] = sum(pattern[j]*signal[i+j]              |
//|                 * negative lags are stored in R[N..N+M-2],       |
//|                   R[N+M-1-i] = sum(pattern[j]*signal[-i+j]       |
//| NOTE:                                                            |
//|     It is assumed that pattern domain is [0..M-1]. If Pattern is |
//| non-zero on [-K..M-1],  you can still use this subroutine, just  |
//| shift result by K.                                               |
//+------------------------------------------------------------------+
void CCorr::CorrR1D(double &signal[],const int n,double &pattern[],
                    const int m,double &r[])
  {
//--- create variables
   int i=0;
   int i_=0;
   int i1_=0;
//--- create arrays
   double p[];
   double b[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- allocation
   ArrayResize(p,m);
//--- copy
   for(i=0; i<=m-1; i++)
      p[m-1-i]=pattern[i];
//--- function call
   CConv::ConvR1D(p,m,signal,n,b);
//--- allocation
   ArrayResize(r,m+n-1);
   i1_=m-1;
//--- copy
   for(i_=0; i_<=n-1; i_++)
      r[i_]=b[i_+i1_];
//--- check
   if(m+n-2>=n)
     {
      i1_=-n;
      for(i_=n; i_<=m+n-2; i_++)
         r[i_]=b[i_+i1_];
     }
  }
//+------------------------------------------------------------------+
//| 1-dimensional circular real cross-correlation.                   |
//| For given Pattern/Signal returns corr(Pattern,Signal) (circular).|
//| Algorithm has linearithmic complexity for any M/N.               |
//| IMPORTANT:                                                       |
//|     for  historical reasons subroutine accepts its parameters in |
//|     reversed order: CorrR1DCircular(Signal, Pattern) = Pattern x |
//|     Signal (using traditional definition of cross-correlation,   |
//|     denoting cross-correlation as "x").                          |
//| INPUT PARAMETERS                                                 |
//|     Signal  -   array[0..N-1] - real function to be transformed, |
//|                 periodic signal containing pattern               |
//|     N       -   problem size                                     |
//|     Pattern -   array[0..M-1] - real function to be transformed, |
//|                 non-periodic pattern to search withing signal    |
//|     M       -   problem size                                     |
//| OUTPUT PARAMETERS                                                |
//|     R   -   convolution: A*B. array[0..M-1].                     |
//+------------------------------------------------------------------+
void CCorr::CorrR1DCircular(double &signal[],const int m,double &pattern[],
                            const int n,double &c[])
  {
//--- create variables
   int i1=0;
   int i2=0;
   int i=0;
   int j2=0;
   int i_=0;
   int i1_=0;
//--- create arrays
   double p[];
   double b[];
//--- check
   if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
      return;
//--- normalize task: make M>=N,
//--- so A will be longer (at least - not shorter) that B.
   if(m<n)
     {
      //--- allocation
      ArrayResize(b,m);
      //--- initialization
      for(i1=0; i1<=m-1; i1++)
         b[i1]=0;
      i1=0;
      //--- calculation
      while(i1<n)
        {
         //--- change values
         i2=MathMin(i1+m-1,n-1);
         j2=i2-i1;
         i1_=i1;
         for(i_=0; i_<=j2; i_++)
            b[i_]=b[i_]+pattern[i_+i1_];
         i1=i1+m;
        }
      //--- function call
      CorrR1DCircular(signal,m,b,m,c);
      //--- exit the function
      return;
     }
//--- Task is normalized
   ArrayResize(p,n);
   for(i=0; i<=n-1; i++)
      p[n-1-i]=pattern[i];
//--- function call
   CConv::ConvR1DCircular(signal,m,p,n,b);
//--- allocation
   ArrayResize(c,m);
   i1_=n-1;
//--- copy
   for(i_=0; i_<=m-n; i_++)
      c[i_]=b[i_+i1_];
//--- check
   if(m-n+1<=m-1)
     {
      i1_=-(m-n+1);
      for(i_=m-n+1; i_<=m-1; i_++)
         c[i_]=b[i_+i1_];
     }
  }
//+------------------------------------------------------------------+
//| Fast Hartley Transform                                           |
//+------------------------------------------------------------------+
class CFastHartleyTransform
  {
public:
   static void       FHTR1D(double &a[],const int n);
   static void       FHTR1DInv(double &a[],const int n);
  };
//+------------------------------------------------------------------+
//| 1-dimensional Fast Hartley Transform.                            |
//| Algorithm has O(N*logN) complexity for any N (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A  -  array[0..N-1] - real function to be transformed        |
//|     N  -  problem size                                           |
//| OUTPUT PARAMETERS                                                |
//|     A  -  FHT of a input array, array[0..N-1],                   |
//|           A_out[k]=sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)),|
//|           j=0..N-1)                                              |
//+------------------------------------------------------------------+
void CFastHartleyTransform::FHTR1D(double &a[],const int n)
  {
//--- create a variable
   int i=0;
//--- create array
   complex fa[];
//--- object of class
   CFtPlan plan;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- Special case: N=1,FHT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
      return;
//--- Reduce FHt to real FFT
   CFastFourierTransform::FFTR1D(a,n,fa);
//--- change values
   for(i=0; i<=n-1; i++)
      a[i]=fa[i].real-fa[i].imag;
  }
//+------------------------------------------------------------------+
//| 1-dimensional inverse FHT.                                       |
//| Algorithm has O(N*logN) complexity for any N (composite or prime)|
//| INPUT PARAMETERS                                                 |
//|     A   -   array[0..N-1] - complex array to be transformed      |
//|     N   -   problem size                                         |
//| OUTPUT PARAMETERS                                                |
//|     A   -   inverse FHT of a input array, array[0..N-1]          |
//+------------------------------------------------------------------+
void CFastHartleyTransform::FHTR1DInv(double &a[],const int n)
  {
//--- create a variable
   int i=0;
//--- check
   if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
      return;
//--- Special case: N=1,iFHT is just identity transform.
//--- After this block we assume that N is strictly greater than 1.
   if(n==1)
      return;
//--- Inverse FHT can be expressed in terms of the FHT as
//---     invfht(x)=fht(x)/N
   FHTR1D(a,n);
//--- change values
   for(i=0; i<=n-1; i++)
      a[i]=a[i]/n;
  }
//+------------------------------------------------------------------+