//-----------------------------------------------------------------------------------
//                                                                     TSAnalysis.mqh
//                                                                      2012, victorg
//                                                                http://www.mql5.com
//-----------------------------------------------------------------------------------
#property copyright "2012, victorg"
#property link      "http://www.mql5.com"
#property version   "2.00"

#include <Object.mqh>

#import "shell32.dll"
int ShellExecuteW(int hwnd,string lpOperation,string lpFile,string lpParameters,
                  string lpDirectory,int nShowCmd);
#import
#import "kernel32.dll"
int DeleteFileW(string lpFileName);
int MoveFileW(string lpExistingFileName,string lpNewFileName);
#import
//-----------------------------------------------------------------------------------
class TSAnalysis:public CObject
  {
protected:
  double    TS[];                    // Time series
  double    TSort[];                 // Sorted time series
  double    TSCenter[];              // Centered time series ( TS[] - mean )
  int       NumTS;                   // Number of time series data points
  double    MinTS;                   // Minimum time series value
  double    MaxTS;                   // Maximum time series value
  double    Median;                  // Median
  double    Mean;                    // Mean (average)
  double    Var;                     // Variance
  double    uVar;                    // Unbiased variance
  double    StDev;                   // Standard deviation
  double    uStDev;                  // Unbiaced standard deviation
  double    Skew;                    // Skewness
  double    Kurt;                    // Kurtosis
  double    ExKurt;                  // Excess Kurtosis
  double    JBTest;                  // Jarque-Bera test
  double    JBpVal;                  // JB test p-value
  double    AJBTest;                 // Adjusted Jarque-Bera test
  double    AJBpVal;                 // AJB test p-values
  double    maxOut;                  // Sequence Plot. Border of outliers
  double    minOut;                  // Sequence Plot. Border of outliers
  double    XHist[];                 // Histogram. X-axis
  double    YHist[];                 // Histogram. Y-axis
  double    Xnpp[];                  // Normal Probability Plot. X-axis
  int       NLags;                   // Number of lags for ACF and PACF Plot
  double    ACF[];                   // Autocorrelation function (correlogram)
  double    UPLim;                   // ACF. Upper limit (5% significance level)
  double    LOLim;                   // ACF. Lower limit (5% significance level)
  double    CBup[];                  // ACF. Upper limit (confidence bands)
  double    CBlo[];                  // ACF. Lower limit (confidence bands)
  double    Spect[];                 // ACF Spectrum. Y-axis
  double    PACF[];                  // Partial autocorrelation function
  int       IP;                      // Autoregressive model order
  double    ARSp[];                  // AR Spectrum. Y-axis
  
public:
            TSAnalysis();
  void      Calc(double &ts[]);      // Calculation

protected:
  double    ndtri(double y);         // Inverse of Normal distribution function
  void      LevinsonRecursion(const double &R[],double &A[],double &K[]);
  void      fht(double &f[], ulong ldn); // Fast Hartley Transform
  virtual void show();               // Show result
  };
//-----------------------------------------------------------------------------------
// Constructor
//-----------------------------------------------------------------------------------
void TSAnalysis::TSAnalysis()
  {

  }
//-----------------------------------------------------------------------------------
// Calculate and draw
//-----------------------------------------------------------------------------------
void TSAnalysis::Calc(double &ts[])
  {
  int i,k,m,n,p;
  double sum2,sum3,sum4,a,b,c,v,delta;
  double cor[],ar[],tdat[];
  
  NumTS=ArraySize(ts);                          // Number of time series data points
  if(NumTS<8)                                   // Number of data points is too small
    {
    Print("TSAnalysis: Error. Number of TS data points is too small!");
    return;
    }
  ArrayResize(TS,NumTS);
  ArrayCopy(TS,ts);                                  // Time series
  ArrayResize(TSort,NumTS);
  ArrayCopy(TSort,ts);
  ArraySort(TSort);                                  // Sorted time series
  MinTS=TSort[0];                                    // Minimum time series value
  MaxTS=TSort[NumTS-1];                              // Maximum time series value
  
  i=(NumTS-1)/2;
  Median=TSort[i];                                   // Median
  if((NumTS&0x01)==0)Median=(Median+TSort[i+1])/2.0; // Median
  
  Mean=0; sum2=0; sum3=0; sum4=0;
  for(i=0;i<NumTS;i++)
    {
    n=i+1; delta=TS[i]-Mean;
    a=delta/n; Mean+=a;                              // Mean (average)
    sum4+=a*(a*a*delta*i*(n*(n-3.0)+3.0)+6.0*a*sum2-4.0*sum3); // sum of fourth degree
    b=TS[i]-Mean;
    sum3+=a*(b*delta*(n-2.0)-3.0*sum2);                  // sum of third degree
    sum2+=delta*b;                                   // sum of second degree
    }
  if(sum2<1.e-250)                                   // variance is too small
    {
    Print("TSAnalysis: Error. The variance is too small or zero!");
    return;
    }
  ArrayResize(TSCenter,NumTS);
  for(i=0;i<NumTS;i++)TSCenter[i]=TS[i]-Mean;        // Centered time series
  Var=sum2/NumTS;                                    // Variance
  uVar=sum2/(NumTS-1);                               // Unbiased variance
  StDev=MathSqrt(Var);                               // Standard deviation
  uStDev=MathSqrt(uVar);                             // Unbiased standard deviation
  Skew=MathSqrt(NumTS)*sum3/sum2/MathSqrt(sum2);     // Skewness
  Kurt=NumTS*sum4/sum2/sum2;                         // Kurtosis
  ExKurt=Kurt-3;                                     // Excess kurtosis
  JBTest=(NumTS/6.0)*(Skew*Skew+ExKurt*ExKurt/4);    // Jarque-Bera test
  JBpVal=MathExp(-JBTest/2.0);                       // JB test p-value
  a=6*(NumTS-2.0)/(NumTS+1.0)/(NumTS+3.0);
  b=3*(NumTS-1.0)/(NumTS+1.0);
  AJBTest=Skew*Skew/a+(Kurt-b)*(Kurt-b)/             // Adjusted Jarque-Bera test
          (24.0*NumTS*(NumTS-2.0)*(NumTS-3.0)/(NumTS+1.0)/
          (NumTS+1.0)/(NumTS+3.0)/(NumTS+5.0));
  AJBpVal=MathExp(-AJBTest/2.0);                     // AJB test p-value
  
  // Time Series Plot. Y=TS[],line1=maxOut,line2=Mean,line3=minOut
  delta=(1.55+0.8*MathLog10(NumTS/10.0)*MathSqrt(Kurt-1))*StDev;
  maxOut=Mean+delta;                          // Time Series Plot. Border of outliers
  minOut=Mean-delta;                          // Time Series Plot. Border of outliers
  
  // Histogram. X=XHist[],Y=YHist[]
  n=(int)MathRound((Kurt+1.5)*MathPow(NumTS,0.4)/6.0);
  if((n&0x01)==0)n--; if(n<5)n=5;                    // Number of bins
  ArrayResize(XHist,n);
  ArrayResize(YHist,n);
  ArrayInitialize(YHist,0.0);
  a=MathAbs(TSort[0]-Mean); b=MathAbs(TSort[NumTS-1]-Mean);
  if(a<b)a=b; v=Mean-a;
  delta=2.0*a/n;
  for(i=0;i<n;i++)XHist[i]=(v+(i+0.5)*delta-Mean)/StDev; // Histogram. X-axis
  for(i=0;i<NumTS;i++)
    {
    k=(int)((TS[i]-v)/delta);
    if(k>(n-1))k=n-1;
    YHist[k]++;
    }
  for(i=0;i<n;i++)YHist[i]=YHist[i]/NumTS/delta*StDev;   // Histogram. Y-axis
  
  // Normal Probability Plot. X=Xnpp[],Y=TSort[]
  ArrayResize(Xnpp,NumTS);
  Xnpp[NumTS-1]=MathPow(0.5,1.0/NumTS);
  Xnpp[0]=1-Xnpp[NumTS-1];
  a=NumTS+0.365;
  for(i=1;i<(NumTS-1);i++)Xnpp[i]=(i+0.6825)/a;
  for(i=0;i<NumTS;i++)Xnpp[i]=ndtri(Xnpp[i]);      // Normal Probability Plot. X-axis

  // Autocorrelation function (correlogram)
  NLags=(int)MathRound(10*MathLog10(NumTS));
  if(NLags>NumTS/5)NLags=NumTS/5;
  if(NLags<3)NLags=3;                         // Number of lags for ACF and PACF Plot

  IP=NLags*5;
  if(IP>NumTS*0.7)IP=(int)MathRound(NumTS*0.7);        // Autoregressive model order
  
  ArrayResize(cor,IP);
  ArrayResize(ar,IP);
  ArrayResize(tdat,IP);
  a=0;
  for(i=0;i<NumTS;i++)a+=TSCenter[i]*TSCenter[i];    
  for(i=1;i<=IP;i++)
    {  
    c=0;
    for(k=i;k<NumTS;k++)c+=TSCenter[k]*TSCenter[k-i];
    cor[i-1]=c/a;                                      // Autocorrelation
    } 
  
  LevinsonRecursion(cor,ar,tdat);                      // Levinson-Durbin recursion
  
  ArrayResize(ACF,NLags);
  ArrayCopy(ACF,cor,0,0,NLags);                        // ACF
  ArrayResize(PACF,NLags);
  ArrayCopy(PACF,tdat,0,0,NLags);                      // PACF
  
  UPLim=1.96/MathSqrt(NumTS);                  // Upper limit (5% significance level)
  LOLim=-UPLim;                                // Lower limit (5% significance level)
  ArrayResize(CBup,NLags);
  ArrayResize(CBlo,NLags);
  a=0;
  for(i=0;i<NLags;i++)
    {
    a+=ACF[i]*ACF[i];  
    CBup[i]=1.96*MathSqrt((1+2*a)/NumTS);           // Upper limit (confidence bands)
    CBlo[i]=-CBup[i];                               // Lower limit (confidence bands)
    }
  
  // Spectrum Plot
  n=320;                                              // Number of X-points
  ArrayResize(Spect,n);
  v=M_PI/n;
  for(i=0;i<n;i++)
    {
    a=i*v; b=0;
    for(k=0;k<NLags;k++)b+=((double)NLags-k)/(NLags+1.0)*ACF[k]*MathCos(a*(k+1));
    Spect[i]=2.0*(1+2*b);                             // Spectrum Y-axis
    }

  // AR Spectral Estimates Plot (maximum entropy method)
  p=12;                                               // n = 2**p = 4096
  n=((ulong)1<<p);                                    // Number of X-points
  m=n<<1;
  ArrayResize(ARSp,n);                                // AR Spectrum. Y-axis
  ArrayResize(tdat,m);
  ArrayInitialize(tdat,0);
  tdat[0]=1;
  for(i=0;i<IP;i++)tdat[i+1]=-ar[i];
  fht(tdat,p+1);                                      // Fast Hartley transform (FHT)
  for(k=1,i=m-1;k<i;++k,--i) tdat[k]=tdat[k]*tdat[k]+tdat[i]*tdat[i];
  tdat[0]=2*tdat[0]*tdat[0];
  ArrayCopy(ARSp,tdat,0,0,n);
  c=-DBL_MAX;
  for(i=0;i<n;i++)
    {
    ARSp[i]=1/ARSp[i];
    if(c<ARSp[i])c=ARSp[i];                           // c = max(ARSp)
    } 
  for(i=0;i<n;i++)                                    // logarithmic scale
    {
    b=ARSp[i]/c;                                      // normalization
    if(b<1e-7)b=1e-7; 
    ARSp[i]=10*MathLog10(b);                          // dB
    }

  // Vizualization
  show();
  
  return;
  }
//-----------------------------------------------------------------------------------
// virtual. Show result
//-----------------------------------------------------------------------------------
void TSAnalysis::show()
  {
  int i,k,n,fhandle;
  double a,x;
  string str="",path,datafile;
  
  str=StringFormat("var NumTS=%.12g;\n",NumTS);
  str+=StringFormat("var MinTS=%.12g;\n",MinTS);
  str+=StringFormat("var MaxTS=%.12g;\n",MaxTS);
  str+=StringFormat("var Mean=%.12g;\n",Mean);
  str+=StringFormat("var maxOut=%.12g;\n",maxOut);
  str+=StringFormat("var minOut=%.12g;\n",minOut);
  
  // Text
  str+="var s=[";
  str+=StringFormat("'Mean','=% g',",Mean);
  str+=StringFormat("'Var','=% g',",Var);
  str+=StringFormat("'StDev','=% g',",StDev);
  str+=StringFormat("'uVar','=% g',",uVar);
  str+=StringFormat("'uStDev','=% g',",uStDev);
  str+=StringFormat("'Skew','=% g',",Skew);
  str+=StringFormat("'ExKurt','=% g',",ExKurt);
  str+="'','','Jarque-Bera test:','',";
  str+=StringFormat("'JB','=% g',",JBTest);
  str+=StringFormat("'p-val','=% f',",JBpVal);
  str+="'','','Adjusted J-B test:','',";
  str+=StringFormat("'AJB','=% g',",AJBTest);
  str+=StringFormat("'p-val','=% f'",AJBpVal);
  str+="];\n";
  
  // Time Series Plot
  k=0;
  str+="var TS=[";
  for(i=0;i<NumTS-1;i++)
    {
    str+=StringFormat("%.12g,",TS[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("%.12g];\n",TS[NumTS-1]);
  
  // Lag Plot
  k=0;
  str+="var LagP=[";
  for(i=0;i<NumTS-2;i++)
    {
    str+=StringFormat("[%.12g,%.12g],",TS[i],TS[i+1]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%.12g,%.12g]];\n",TS[NumTS-2],TS[NumTS-1]);
  
  // Histogram
  n=ArraySize(XHist);
  str+="var HistP=[";
  k=0;
  for(i=0;i<n-1;i++)
    {
    str+=StringFormat("[%.4f,%.4f],",XHist[i],YHist[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%.4f,%.4f]];\n",XHist[n-1],YHist[n-1]);
  a=XHist[1]-XHist[0];
  str+="var Nd=[";                                   //Normal distribution line
  n=n+4;
  k=0;
  for(i=0;i<n-1;i++)
    {
    x=a*(i-n/2);
    str+=StringFormat("[%.4f,%.4f],",x,MathExp(-x*x/2)/MathSqrt(2*M_PI));
    if(10<k++){k=0; str+="\n";}
    }
  x=a*(n-1-n/2);
  str+=StringFormat("[%.4f,%.4f]];\n",x,MathExp(-x*x/2)/MathSqrt(2*M_PI));
  
  // Normal Probability Plot
  k=0;
  str+="var Xnpp=[";
  for(i=0;i<NumTS-1;i++)
    {
    str+=StringFormat("[%.12g,%.12g],",Xnpp[i],TSort[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%.12g,%.12g]];\n",Xnpp[NumTS-1],TSort[NumTS-1]);
  str+=StringFormat("var Ndl=[[%.12g,%.12g],[%.12g,%.12g]];\n",           // Line
                Xnpp[0],Xnpp[0]*StDev+Mean,Xnpp[NumTS-1],Xnpp[NumTS-1]*StDev+Mean);
  
  // Autocorrelation function (correlogram)
  str+=StringFormat("var NLags=%i;\n",NLags);        // ACF. Number of lags
  k=0;
  str+="var ACF=[";
  for(i=0;i<NLags-1;i++)
    {
    str+=StringFormat("[%i,%.4f],",i+1,ACF[i]);     // ACF
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%i,%.4f]];\n",NLags,ACF[NLags-1]);
  str+=StringFormat("var UPLim=%.4f;\n",UPLim);
  str+=StringFormat("var CBup=[[0,%.4f],",UPLim);
  k=0;
  for(i=0;i<NLags;i++)
    {
    str+=StringFormat("[%i,%.4f],",i+1,CBup[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%i,%.4f]];\n",NLags+1,CBup[NLags-1]);
  str+=StringFormat("var CBlo=[[0,%.4f],",LOLim);
  k=0;
  for(i=0;i<NLags;i++)
    {
    str+=StringFormat("[%i,%.4f],",i+1,CBlo[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%i,%.4f]];\n",NLags+1,CBlo[NLags-1]);
  
  // Spectrum
  k=0;
  n=ArraySize(Spect);
  str+="var Spec=[";
  for(i=0;i<n-1;i++)
    {
    str+=StringFormat("%.4f,",Spect[i]);
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("%.4f];\n",Spect[n-1]);
  
  // Partial autocorrelation function
  k=0;
  str+="var PACF=[";
  for(i=0;i<NLags-1;i++)
    {
    str+=StringFormat("[%i,%.4f],",i+1,PACF[i]);    // PACF
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("[%i,%.4f]];\n",NLags,PACF[NLags-1]);
  
  // Autoregressive Spectral Estimates
  str+=StringFormat("var IP=%i;\n",IP);              // AR model order
  k=0;
  n=ArraySize(ARSp);
  str+="var ARSp=[";
  for(i=0;i<n-1;i++)
    {
    str+=StringFormat("%.3f,",ARSp[i]);             // AR Spectrum. Y-axis
    if(10<k++){k=0; str+="\n";}
    }
  str+=StringFormat("%.3f];\n",ARSp[n-1]);          // AR Spectrum. Y-axis
  
  // Save
  fhandle=FileOpen("TSDat.txt",FILE_WRITE|FILE_TXT|FILE_ANSI);
  if(fhandle==INVALID_HANDLE)return;
  FileWriteString(fhandle,str);
  FileClose(fhandle);
  
  //Move data file and Execute
  datafile=TerminalInfoString(TERMINAL_DATA_PATH)+"\\MQL5\\Files\\TSDat.txt";
  path=TerminalInfoString(TERMINAL_DATA_PATH)+"\\MQL5\\Scripts\\TSAnalysis\\";
  DeleteFileW(path+"TSDat.txt");
  MoveFileW(datafile,path+"TSDat.txt");
  ShellExecuteW(NULL,"open",path+"TSA.htm",NULL,NULL,1);

  return;
  }
//-----------------------------------------------------------------------------------
// Inverse of Normal distribution function
// Prototype:
// Cephes Math Library Release 2.8: June, 2000
// Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
//-----------------------------------------------------------------------------------
double TSAnalysis::ndtri(double y0)
  {
  static double s2pi =2.50662827463100050242E0; // sqrt(2pi)
  static double P0[5]={-5.99633501014107895267E1,  9.80010754185999661536E1,
                       -5.66762857469070293439E1,  1.39312609387279679503E1,
                       -1.23916583867381258016E0};
  static double Q0[8]={ 1.95448858338141759834E0,  4.67627912898881538453E0,
                        8.63602421390890590575E1, -2.25462687854119370527E2,
                        2.00260212380060660359E2, -8.20372256168333339912E1,
                        1.59056225126211695515E1, -1.18331621121330003142E0};
  static double P1[9]={ 4.05544892305962419923E0,  3.15251094599893866154E1,
                        5.71628192246421288162E1,  4.40805073893200834700E1,
                        1.46849561928858024014E1,  2.18663306850790267539E0,
                       -1.40256079171354495875E-1,-3.50424626827848203418E-2,
                       -8.57456785154685413611E-4};
  static double Q1[8]={ 1.57799883256466749731E1,  4.53907635128879210584E1,
                        4.13172038254672030440E1,  1.50425385692907503408E1,
                        2.50464946208309415979E0, -1.42182922854787788574E-1,
                       -3.80806407691578277194E-2,-9.33259480895457427372E-4};
  static double P2[9]={ 3.23774891776946035970E0,  6.91522889068984211695E0,
                        3.93881025292474443415E0,  1.33303460815807542389E0,
                        2.01485389549179081538E-1, 1.23716634817820021358E-2,
                        3.01581553508235416007E-4, 2.65806974686737550832E-6,
                        6.23974539184983293730E-9};
  static double Q2[8]={ 6.02427039364742014255E0,  3.67983563856160859403E0,
                        1.37702099489081330271E0,  2.16236993594496635890E-1,
                        1.34204006088543189037E-2, 3.28014464682127739104E-4,
                        2.89247864745380683936E-6, 6.79019408009981274425E-9};
  double x,y,z,y2,x0,x1,a,b;
  int i,code;

  if(y0<=0.0){Print("Function ndtri() error!"); return(-DBL_MAX);}
  if(y0>=1.0){Print("Function ndtri() error!"); return(DBL_MAX);}

  code=1; y=y0;
  if(y>(1.0-0.13533528323661269189)){y=1.0-y; code=0;}  // 0.135... = exp(-2)
  if(y>0.13533528323661269189)                          // 0.135... = exp(-2)
    {
    y=y-0.5; y2=y*y;
    a=P0[0]; for(i=1;i<5;i++)a=a*y2+P0[i];
    b=y2+Q0[0]; for(i=1;i<8;i++)b=b*y2+Q0[i];
    x=y+y*(y2*a/b);
    x=x*s2pi; 
    return(x);
    }
  x=MathSqrt(-2.0*MathLog(y));
  x0=x-MathLog(x)/x;
  z=1.0/x;
  if(x<8.0)                               // y > exp(-32) = 1.2664165549e-14
    {
    a=P1[0]; for(i=1;i<9;i++)a=a*z+P1[i];
    b=z+Q1[0]; for(i=1;i<8;i++)b=b*z+Q1[i];
    x1=z*a/b;
    }
  else
    {
    a=P2[0]; for(i=1;i<9;i++)a=a*z+P2[i];
    b=z+Q2[0]; for(i=1;i<8;i++)b=b*z+Q2[i];
    x1=z*a/b;
    }
  x=x0-x1;
  if(code!=0)x=-x;

  return(x);
  }
//-----------------------------------------------------------------------------------
// Calculate the Levinson-Durbin recursion for the autocorrelation sequence R[]
// and return the autoregression coefficients A[] and partial autocorrelation
// coefficients K[]
//-----------------------------------------------------------------------------------
void TSAnalysis::LevinsonRecursion(const double &R[],double &A[],double &K[])
  {
  int p,i,m;
  double km,Em,Am1[],err;

  p=ArraySize(R);
  ArrayResize(Am1,p);
  ArrayInitialize(Am1,0);
  ArrayInitialize(A,0);
  ArrayInitialize(K,0);
  km=0;
  Em=1;
  for(m=0;m<p;m++)
    {
    err=0;
    for(i=0;i<m;i++)err+=Am1[i]*R[m-i-1];
    km=(R[m]-err)/Em;
    K[m]=km; A[m]=km;
    for(i=0;i<m;i++)A[i]=(Am1[i]-km*Am1[m-i-1]);
    Em=(1-km*km)*Em;
    ArrayCopy(Am1,A);
    }
  return;
  }
//-----------------------------------------------------------------------------------
// Radix-2 decimation in frequency (DIF) fast Hartley transform (FHT).
// Length is N = 2 ** ldn
//-----------------------------------------------------------------------------------
void TSAnalysis::fht(double &f[], ulong ldn)
  {
  const ulong n = ((ulong)1<<ldn);
  for (ulong ldm=ldn; ldm>=1; --ldm)
    {
    const ulong m = ((ulong)1<<ldm);
    const ulong mh = (m>>1);
    const ulong m4 = (mh>>1);
    const double phi0 = M_PI / (double)mh;
    for (ulong r=0; r<n; r+=m)
      {
      for (ulong j=0; j<mh; ++j)
        {
        ulong t1 = r+j;
        ulong t2 = t1+mh;
        double u = f[t1];
        double v = f[t2];
        f[t1] = u + v;
        f[t2] = u - v;
        }
      double ph = 0.0;
      for (ulong j=1; j<m4; ++j)
        {
        ulong k = mh-j;
        ph += phi0;
        double s=MathSin(ph);
        double c=MathCos(ph);
        ulong t1 = r+mh+j;
        ulong t2 = r+mh+k;
        double pj = f[t1];
        double pk = f[t2];
        f[t1] = pj * c + pk * s;
        f[t2] = pj * s - pk * c;
        }
      }
    }
  if(n>2)
    {
    ulong r = 0;
    for (ulong i=1; i<n; i++)
      {
      ulong k = n;
      do {k = k>>1; r = r^k;} while ((r & k)==0);
      if (r>i) {double tmp = f[i]; f[i] = f[r]; f[r] = tmp;}
      }
    }
  }
//-----------------------------------------------------------------------------------