#include<stdio.h>
#include<stdlib.h>
#include<math.h>



/* The following are parameters.  They can be changed.  
   h and eta are set as follows:

   h= hmult* 10^{ hpower}
   eta = etamult* 10^{etapower}
*/

#define kfinal 200
#define k0 2
#define hpower -8
#define hmult 2
#define etapower -15
#define etamult 5



/* The following are the variables for the numerical solution 
   of the system of differential equations.  They match the 
   notation from the paper.  

   X[i] = \tilde{x}_i(t)
   g[k] = g(k) = transition parameter for phase k
   Xp[i] = holds the next value of X[i] during the
           calculations
   k = current Phase
*/

double X[kfinal+1],Xp[kfinal+1],g[kfinal+1],t,h;
int i,k;



/* The following variables are used in bounding the errors.
   Again, they match the notation from the paper.

   E[i] = upper bound on the absolute value of
          e_i(t) = x_i(t) - \tilde{x}_i(t) 
   Ep[i] = holds the next value of E[i] during the
           calculations
   f = upper bound on the absolute value of 
       1 - \sum_{i=1}^k i \tilde{x}_i(t) 
             = 1  - \sum_{i=1}^k i X[i]
   maxe = maximum of E[i] over i =1, ..., k
            (this is diagnostic only. it is output
             but not used anywhere in the error bounds)
   eta = maximum error in a floating point operation
*/

double f, E[kfinal+1], Ep[kfinal+1], maxe, eta;




/*  These functions due the main calculations

    dxi returns F_i(\tilde{ \bf x})
    dx0 returns F_0(\tilde{ \bf x}) 
    errxi returns upper bound on the absolute value of
          e_i(t+h), following the methods discussed in 
          the paper
    errx0 return an upper bound on the absolute value
          of e_0(t+h).
*/

double dxi(double *yy, long k, long i);
double dx0(double *yy, long k);
double errxi(double *ee, double *yy, long k, long i, double h, double eta, double f);
double errx0(double *ee, double *yy, long k, double h, double eta, double f);




int main(void)
{

  /* Initialization of everything */

  h = hmult*exp( log(10.0)* hpower);
  eta = etamult*exp( log(10.0)* etapower);

  k = k0;
  t = 0;
  f = 0;

  X[0]=0;
  X[1]=1;
  for (i=2;i<=kfinal;++i) X[i] = 0.0;

  for (i=0;i<=kfinal;++i) E[i] = 0.0;

  for (i=1;i<kfinal;++i) g[i] = 0.5 + sqrt(1.0/(2*i));
  g[kfinal] = 0.5;



  while(k <= kfinal)    /* Phase loop */
    {

      while(  (X[0] >=  t*g[k]) && (X[0]-E[0] >= t*0.5) )   /* time loop */
	{

          Xp[0] = X[0] + h*dx0(X,k);
          for(i=1; i<= k; ++i) Xp[i] = X[i] + h*dxi(X,k,i);

          Ep[0] = errx0(E,X,k,h,eta,f);
          for(i=1; i<= k; ++i) Ep[i] = errxi(E,X,k,i,h,eta,f);

          for(i=0; i<= k; ++i) 
	    {
              X[i] = Xp[i];
              E[i] = Ep[i];  
	    }


	  /*  Calculation of f for use in next step */

          f=0;
          for(i=1; i<= k; ++i)  f = f + i*X[i];
          f = fabs( f - 1);
          f = f + 2*k*eta;



	  t=t+h;

	}   /* end of time loop */


      /* Output of diagnostic information */

      maxe=0;
      for(i=1; i <=k; ++i) if( E[i]>= maxe ) maxe=E[i];
      printf("k=%i  t=%8f  ratio=%8f  maxe=%le  f= %le  E[0]= %le  "
                                         ,k, t, X[0]/t, maxe, f, E[0]);
      printf("\n");


      k=k+1;

   }    /* end of Phase loop */ 

}   



double dxi(double *yy, long k, long i)
{
	long j,ell;
	double sum;

	sum=0;

	for(j=1;j<i;++j)
	{
		ell=i-j;
		sum+= j*yy[j]*ell*yy[ell];
	}
	for(j=1;j<=k-i;++j)
	{
		sum = sum -  2 * j*yy[j]*i*yy[i];
	}
	return sum;
}




double dx0(double *yy,long k)
{
	long j,i;
	double sum;
	sum=0;
	for(j=1;j<k;++j)
		for(i=1;i+j<=k;++i)
			sum+=i*yy[i]*j*yy[j];
	return sum;
}



double errxi(double *ee, double *yy, long k, long i, double h, double eta, double f)
{
	long j,ell;
	double alpha,beta,error,sum;


        /* We begin with equation (25) from the paper */

        sum=0.0;
        for(j=1; j<=k-i; ++j)
	     {
		sum+= yy[j]*j  ;
	     }
        error = ee[i] * (1 - 2*i*h*sum);


   
        /* Then add in equation (27) */

        sum=0.0;
        for(j=1;j<i;++j)
	     {
		ell=i-j;
		sum+= ee[j]*j*ee[ell]*ell;
	     }
        error = error + sum*h;
        sum=0.0;
        for(j=1; j<=k-i; ++j)
	     {
		sum += ee[j]*j;
	     }
        error = error + 2*h*ee[i]*i*sum;



        /* alpha is the quantity from equation (28)   */

        alpha = 0.0;
        if( k <= 2*i -1)
	  {
              for(j=1; j<=k-i; ++j)
		{
		  ell=i-j;
		  alpha += j*ee[j]* fabs(ell*yy[ell] - i* yy[i]);
		}
              for(j=k-i+1; j<=i-1; ++j)
		{
		  ell=i-j;
		  alpha += j*ee[j]* ell*yy[ell];
		}
	  }
        else
	  {
              for(j=1; j<=i-1; ++j)
		{
		  ell=i-j;
		  alpha += j*ee[j]* fabs(ell*yy[ell] - i* yy[i]);
		}
              sum=0.0;
	      for(j=i; j<=k-i; ++j);
	        {
		  sum += j *ee[j];
		}
	      alpha += sum* i *yy[i];
	  }
        


        /*beta ia the quantity from equation (29) */

        beta = 0.0;
	for(j=1;j<i;++j)
	     {
		ell=i-j;
		beta+= yy[j]*j*ee[ell]*ell;
	     }
        sum= 0.0;
        for(j=k-i+1; j<=k; ++j)
	     {
		sum = sum +  ee[j]*j ;
	     }
        beta = beta + i*yy[i]*(sum + f); 



        /* We then add 2h times the smaller of alpha and beta */

        if( alpha <= beta ) 
	  {
	    error = error + 2*h* alpha;
	  }
        else 
	  {
            error = error + 2*h* beta;
	  }


       
        /* We finish by adding in the truncation and rounding errors 
           We use 13*eta to account for potetial rounding errors in
           both the calaculation of X[i] and E[i].
        */

        error = error + 4*h*h*k + 13*eta;
	return error;
}


double errx0(double *ee, double *yy, long k, double h, double eta, double f )
{
	long j,i;
	double sum, alpha, beta;

	sum=0.0;
	for(j=1;j<k;++j)
		for(i=1;i+j<=k;++i)
			sum+= i*j*ee[i]*ee[j];
        sum = sum*h;

	for(j=1;j<k;++j)
	  {
	    alpha=0;
 	    for(i=1;i<=k-j;++i)  alpha+= i*ee[i];

            beta=0;
 	    for(i=k-j+1;i<=k;++i)  beta+= i*ee[i];
            beta += f;


            if( alpha <= beta ) 
	       {
	        sum = sum + 2*h*j*yy[j]*alpha;
	       }
            else 
	       {
                sum = sum + 2*h*j*yy[j]*beta;
	       }
	  }

        sum = ee[0] + sum +  4*h*h*k + 9*eta;
	return sum;
}