/* aquad_chebyshev.c  compute fit using adaptive integration
 *
 *  T0(x)=1  T1(x)=x   Tn+1(x)=2*x*Tn(x) - Tn-1(x)
 *  T2(x)=2 x^2 - 1   T3(x)=4 x^3 - 3 x  T4(x)=8 x^4 -8 x^2 + 1
 *  c_i = 2/Pi int -1 to 1 of f(x)*T_i(x)dx/sqrt(1-x^2)
 *  F(x) = c_0/2 sum i=1 to n c_i*T_i(x)
 */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(a)  ((a)<0.0?(-(a)):(a))

static double aquad3(double f(double x), double xmin, double xmax, double eps);
static double peval(int n, double x, double c[]); /* evaluate a polynomial */
static double teval(int n, double x); /* recursive evaluation of Tn(x) */

static double Pi=3.141592653589793;

#define N 9              /* change to max polynomial size */
static double T[N][N+1]; /* coefficients of Tn(x) */
static double C[N];   /* c_i for function */
static double Cm[N];   /* maple values */
static double rmserr, maxerr;

/* function to fit in range -1 to 1 */
static double f(double x)
{
  return exp(x);
}

static int iorder; /* set below for g(x) */
/* create function g(x) = f(x)T_i(x)/sqrt(1.0-x*x) */
static double g(double x)
{
  double a, xx;
  a = peval(iorder, x, &T[iorder][0]);
  xx = x * x;
  return (f(x)*a)/sqrt(1.0-xx);
}

int main(int argc, char * argv[])
{
  int i, j, k, n;
  double a=-1.0;
  double b=1.0;
  double x, exact;
  double approx=0.0;
  double err=0.0;
  double T5, T6, T7;

  printf("aquad_chebyshev.c running \n");
  /* generate family of Chebyshev polynomials Tn(x) */
  /* as T[degree][coefficients] */
  /* build coefficients of Tn(x) */
  for(n=0; n<N; n++) for(i=0; i<=N; i++) T[n][i]=0.0;
  T[0][0]=1.0;
  printf("T[%d][%d]=%g \n\n", 0, 0, T[0][0]);
  T[1][1]=1.0; /* start for recursion */
  printf("T[%d][%d]=%g \n", 1, 0, T[1][0]);
  printf("T[%d][%d]=%g \n\n", 1, 1, T[1][1]);
  for(n=2; n<N; n++)
  {
    for(i=0; i<=n; i++)
    {
      T[n][i]=T[n][i]-T[n-2][i];
      T[n][i+1]=T[n][i+1]+2.0*T[n-1][i];
      printf("T[%d][%d]=%g \n", n, i, T[n][i]);
    }
    printf("\n");
  }
  k = 16;
  for(i=0; i<=k; i++) /* sample x's */
  {
    x = (double)(i-k/2)/(double)(k/2);
    T5 = teval(5,x);
    T6 = teval(6,x);
    T7 = teval(7,x);
    printf("teval i=%d, x=%8.5f, T5=%g, T6=%g T7=%g \n",
            i, x, T5, T6, T7);
  }
  for(i=0; i<=k; i++) /* sample x's */
  {
    x = (double)(i-k/2)/(double)(k/2);
    T5 = peval(5,x,&T[5][0]);
    T6 = peval(6,x,&T[6][0]);
    T7 = peval(7,x,&T[7][0]);
    printf("peval i=%d, x=%8.5f, T5=%g, T6=%g T7=%g \n",
            i, x, T5, T6, T7);
  }
  for(i=0; i<=k; i++) /* sample x's */
  {
    x = (double)(i-k/2)/(double)(k/2);
    T5 = f(x);
    if(x<=-1.0) x=-0.99999;
    if(x>=1.0)  x= 0.99999;
    iorder = 0;
    T6 = g(x);
    iorder = 1;
    T7 = g(x);
    printf("i=%d, x=%8.5f, f(x)=%g, g0(x)=%g, g1(x)=%g \n",
            i, x, T5, T6, T7);
  }
  


  Cm[0]=1.2660658777;
  Cm[1]=1.1303182079;
  Cm[2]=0.2714953395;
  Cm[3]=0.0443368498;
  Cm[4]=0.0054742404;
  Cm[5]=0.0005429263;
  Cm[6]=0.0000449773;
  Cm[7]=0.0000031984;
  for(i=0; i<=7; i++) printf("Cm[%d]=%13.10f \n", i, Cm[i]);
  printf("\n");
  k = 16;
  rmserr = 0.0;
  maxerr = 0.0;
  for(i=0; i<=k; i++) /* sample x's */
  {
    x = (double)(i-k/2)/(double)(k/2);
    approx = Cm[0];
    for(j=1; j<=7; j++)
    {
      approx = approx + Cm[j]*peval(j,x,&T[j][0]);
    }
    err = f(x)-approx;
    err = abs(err);
    printf("maple i=%d, x=%8.5f, approx=%g, exact=%g err=%g \n",
            i, x, approx, f(x), err);
    rmserr = rmserr + err*err;
    if(err>maxerr) maxerr = err;
  }
  rmserr = sqrt(rmserr/(double)(k+1));
  printf("maple  maxerr=%g, rmserr=%g \n", maxerr, rmserr);

  /* generate C[i] for function f(x)                      */
  /*  c_i = 2/Pi int -1 to 1 of f(x)*T_i(x)dx/sqrt(1-x^2) */
  /*  create function g(x) = f(x)T_i(x)/sqrt(1.0-x*x)     */
  printf("generating coefficients \n");

  for(n=0; n<N; n++)
  {
    iorder = n; /* needed by g(x) */
    C[n] = aquad3(g, -0.999999, 0.999999, 0.0001);
    printf("n=%d, max=0.999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.9999999, 0.9999999, 0.00001);
    printf("n=%d, max=0.9999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.99999999, 0.99999999, 0.00001);
    printf("n=%d, max=0.99999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.999999999, 0.999999999, 0.000001);
    printf("n=%d, max=0.999999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.9999999999, 0.9999999999, 0.0000001);
    printf("n=%d, max=0.9999999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.99999999999, 0.99999999999, 0.00000001);
    printf("n=%d, max=0.99999999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.999999999999, 0.999999999999, 0.00000001);
    printf("n=%d, max=0.999999999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = aquad3(g, -0.9999999999999, 0.9999999999999, 0.00000001);
    printf("n=%d, max=0.9999999999999, C[%d]=%g \n", n, n, C[n]);
    C[n] = 2.0*C[n]/Pi;
    if(n==0) C[n] = C[n]/2.0;
    printf("n=%d, scaled err=%g       C[%d]=%13.10f \n",
            n, Cm[n]-C[n], n, C[n]);
  }

  /* test fit at many points for each orders -1 < x < 1 */
  printf("test fit based on integration coefficients \n");
  k = 16; /* must be even, does -1 to 1 with extra point 0.0 */
  for(n=2; n<N; n++) /* order of fit */
  {
    rmserr = 0.0;
    for(i=0; i<=k; i++) /* sample x's */
    {
      x = (double)(i-k/2)/(double)(k/2);
      approx = C[0];
      for(j=1; j<=n; j++)
      {
        approx = approx + C[j]*peval(j,x,&T[j][0]);
      }
      err = f(x)-approx;
      err = abs(err);
      printf("n=%d, i=%d, x=%8.5f, approx=%g, exact=%g err=%g \n",
              n, i, x, approx, f(x), err);
      rmserr = rmserr + err*err;
      if(i==0) maxerr = err;
      if(err>maxerr) maxerr = err;
    }
    rmserr = sqrt(rmserr/(double)(k+1));
    printf("n=%d, maxerr=%g, rmserr=%g \n", n, maxerr, rmserr);
  } /* end order loop */
  return 0; 
} /* end main of aquad_chebyshev */

/* peval.c Horners method for evaluating a polynomial */
static double peval(int n, double x, double c[])
{ /* Horners method for evaluating a polynomial   */
  /* an nth order polynomial has n+1 coefficients */
  /* stored in c[0] through c[n]                  */
  /* y = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n  */
  int i;
  double y = c[n]*x;
  if(n<=0) return c[0];
  for(i=n-1; i>0; i--) y = (c[i]+y)*x;
  return y+c[0];
} /* end peval */

/* aquad3  from book page 301, with modifications */
static double stack[500][7];   /* a place to store/retrieve */
static int top, maxtop;        /* top points to where stored */
static int hitop;              /* highest top */
static int funeval=0;          /* count function evaluations */
static void store(double s0, double s1, double s2, double s3,
                  double s4, double s5, double s6)
{
  stack[top][0] = s0;
  stack[top][1] = s1;
  stack[top][2] = s2;
  stack[top][3] = s3;
  stack[top][4] = s4;
  stack[top][5] = s5;
  stack[top][6] = s6;
} /* end store for aquad3 */

static void retrieve(double* s0, double* s1, double* s2, double* s3,
                     double* s4, double* s5, double* s6)
{
  *s0 = stack[top][0];
  *s1 = stack[top][1];
  *s2 = stack[top][2];
  *s3 = stack[top][3];
  *s4 = stack[top][4];
  *s5 = stack[top][5];
  *s6 = stack[top][6];
} /* end retrieve for aquad3 */

static double Sn(double F0, double F1, double F2, double h)
{
  return h*(F0 + 4.0*F1 + F2)/3.0; /* 2h/6 */
} /* end Sn for aquad3 */

static double RS(double F0, double F1, double F2, double F3,
                 double F4, double h) /* 4h/16 */
{
  return h*(14.0*F0 +64.0*F1 + 24.0*F2 + 64.0*F3 + 14.0*F4)/45.0; 
  /* error term  8/945  h^7 f^(8)(c) */
} /* end RS for aquad3 */

double aquad3(double f(double x), double xmin, double xmax, double eps)
{
  double a, b, c, d, e;
  double Fa, Fb, Fc, Fd, Fe;
  double h1, h2, hmin;
  double Sab, Sac, Scb, S2ab;
  double tol, reqtol;
  double val, value;

  maxtop = 499;
  hitop = 0;
  top = 0;
  value = 0.0;
  reqtol = eps;
  tol = reqtol/1.0;
  a = xmin;
  b = xmax;
  funeval = 0;
  h1 = (b-a)/2.0;
  hmin = h1;
  c = a + h1;
  Fa = f(a); funeval++;
  Fc = f(c); funeval++;
  Fb = f(b); funeval++;
  Sab = Sn(Fa, Fc, Fb, h1);
  store(a, Fa, Fc, Fb, h1, tol, Sab);
  top = 1;
  while(top > 0)
  {
    top--;
    retrieve(&a, &Fa, &Fc, &Fb, &h1, &tol, &Sab);
    c = a + h1;
    b = a + 2.0*h1;
    h2 = h1/2;
    d = a + h2;
    e = a + 3.0*h2;
    Fd = f(d); funeval++;
    Fe = f(e); funeval++;
    Sac = Sn(Fa, Fd, Fc, h2);
    Scb = Sn(Fc, Fe, Fb, h2);
    S2ab = Sac + Scb;
    /* printf("S2ab=%15.6f, Sab=%15.6f, err=%g \n", S2ab, Sab, S2ab-Sab); */
    if(abs(S2ab-Sab) < tol || h2 < 1.0e-13)
    {
      val = RS(Fa, Fd, Fc, Fe, Fb, h2);
      value = value + val;
    }
    else
    {
      h1 = h2;
      /* printf("splitting %15.6f to %15.6f to %15.6f \n", a, c, b); */
      tol = tol/2.0;
      store(a, Fa, Fd, Fc, h1, tol, Sac);
      top++;
      store(c, Fc, Fe, Fb, h1, tol, Scb);
      top++;
      if(h1 < hmin) hmin = h1;
    }
    if(top>hitop) hitop = top;
    if(top>=maxtop) break;
  } /* end while */
  printf("aquad3 hitop=%d, funeval=%d, hmin=%g \n", hitop, funeval, hmin);
  return value;
} /* end main of aquad3.c */


static double aquad3x(double f(double x), double xmin, double xmax, double eps)
{
  double a, b, c, d, e;
  double Fa, Fb, Fc, Fd, Fe;
  double h1, h2;
  double Sab, Sac, Scb, S2ab;
  double tol; /* eps */
  double val, value;

  maxtop = 99;
  top = 0;
  value = 0.0;
  tol = eps;
  a = xmin;
  b = xmax;
  h1 = (b-a)/2.0;
  c = a + h1;
  Fa = f(a);
  Fc = f(c);
  Fb = f(b);
  Sab = Sn(Fa, Fc, Fb, h1);
  store(a, Fa, Fc, Fb, h1, tol, Sab);
  top = 1;
  while(top > 0)
  {
    top--;
    retrieve(&a, &Fa, &Fc, &Fb, &h1, &tol, &Sab);
    c = a + h1;
    b = a + 2.0*h1;
    h2 = h1/2;
    d = a + h2;
    e = a + 3.0*h2;
    Fd = f(d);
    Fe = f(e);
    Sac = Sn(Fa, Fd, Fc, h2);
    Scb = Sn(Fc, Fe, Fb, h2);
    S2ab = Sac + Scb;
    if(abs(S2ab-Sab) < tol || h2 < 1.0e-13)
    {
      val = RS(Fa, Fd, Fc, Fe, Fb, h2);
      value = value + val;
    }
    else
    {
      h1 = h2;
      tol = tol/2.0;
      store(a, Fa, Fd, Fc, h1, tol, Sac);
      top++;
      store(c, Fc, Fe, Fb, h1, tol, Scb);
      top++;
    }
    if(top>=maxtop) break;
  } /* end while */
  return value;
} /* end aquad3 */

static double teval(int n, double x) /* recursive evaluation of Tn(x) */
{
  if(n<=0) return 1.0;
  if(n==1) return x;
  return 2.0*x*teval(n-1,x) - teval(n-2,x);
} /* end Teval */

/* end aquad_chebyshev */