/* chebyshev.c   two ways to get coefficients for e^x            */
/* compute x(i) and w(i)=2/Pi  i=1,n   ordinates and weights     */
/* on interval  -1.0 to 1.0 (length is 2.0)                      */
/* use ordinates and weights for Gauss Chebyshev integration     */

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

static double Pi=3.141592653589793;
static double TP[18][18]; /* TP[n][i] */
static double Texp[18][18];  /* for evaluating exp(x) in -1 to 1 */
static double TTexp[18][18]; /* for evaluating exp(x) in -1 to 1 */ 
static double expp[18];
static double P[18];         /* temp polynomial */

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

static double eval(int n, double P[], double x)
{ /* evaluate polynomial P at point x */
  int i;
  double y = P[n]*x;
  for(i=n-1; i>0; i--) y = (y+P[i])*x;
  return y+P[0];
} /* end eval */

static double nfct(int n)
{ /* compute n factorial as a double */
  if(n<=1) return 1.0;
  return (double)n*nfct(n-1);
} /* end nfct */

static void make_expp(int n)
{ /* generate Taylor series polynomial for exp(x) */
  int i;
  for(i=0; i<n; i++) expp[i] = 1.0/nfct(i);
  for(i=0; i<n; i++) printf("expp[%d]=%g \n", i, expp[i]);
} /* end expp */

static void make_TP(int n)
{ /* build coefficients of Chebyshev polynomials Tn(x) */
  int i, j;
  TP[0][0] =  1.0;
  TP[1][0] =  0.0;
  TP[1][1] =  1.0;
  TP[2][0] = -1.0;
  TP[2][1] =  0.0;
  TP[2][2] =  2.0;
  TP[3][0] =  0.0;
  TP[3][1] = -3.0;
  TP[3][2] =  0.0;
  TP[3][3] =  4.0;
  /* Tn(x) = 2xTn-1(x)-Tn-2(x) */
  for(i=4; i<n; i++)
  {
    for(j=0; j<i-1; j++) TP[i][j] = -TP[i-2][j];
    for(j=0; j<i-2; j++) TP[i][j+1] = TP[i][j+1] + 2.0*TP[i-1][j];
    TP[i][i-1] = 2.0*TP[i-1][i-2];
    TP[i][i] = 2.0*TP[i-1][i-1];
  }
  printf("Tn(x) as TP[n][i]*x^i \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<=i; j++) printf("TP[%d][%d]=%g\n", i, j, TP[i][j]);
  }
} /* end make_TP */

static void make_cheb(int n) /* telescoping exp(x) */
{ /* existing polynomial, sum down = cheb sum down
    c0 1   = c0 T0
    c1 x   = c1 T1
    c2 x^2 = c2 1/2(T2 + T0)
    c3 x^3 = c3 1/4(T3 + 3T1)
    c4 x^4 = c4 1/8(T4 + 4T2 + 3T0)
    c5 x^5 = c5 1/16(T5 + 5T3 + 10T1)
    c6 x^6 = c6 1/32(T6 + 6T4 + 15T2 + 10T0)
    c7 x^7 = c7 1/64(T7 + 7T5 + 21T3 + 35T1)
    c8 x^8 = c8 1/128(T8 + 8T6 + 28T4 + 56T2 + 35T0)
    c9 x^9 = c9 1/256(T9 + 9T7 + 36T5 + 84T3 + 126T1)
  c10 X^10 = c10 1/2^-9(T10 + 10T8 + 36+9 + 
              T0    T1    T2    T3    T4    T5   T6   T7   T8   T9 */
  double c[10][10]=
           {{1.0,  0.0,  0.0,  0.0,  0.0,  0.0, 0.0, 0.0, 0.0, 0.0},/*0*/
            {0.0,  1.0,  0.0,  0.0,  0.0,  0.0, 0.0, 0.0, 0.0, 0.0},/*1*/
            {1.0,  0.0,  1.0,  0.0,  0.0,  0.0, 0.0, 0.0, 0.0, 0.0},/*2*/
            {0.0,  3.0,  0.0,  1.0,  0.0,  0.0, 0.0, 0.0, 0.0, 0.0},/*3*/
            {3.0,  0.0,  4.0,  0.0,  1.0,  0.0, 0.0, 0.0, 0.0, 0.0},/*4*/
            {0.0, 10.0,  0.0,  5.0,  0.0,  1.0, 0.0, 0.0, 0.0, 0.0},/*5*/
            {10.0,  0.0, 15.0,  0.0,  6.0,  0.0, 1.0, 0.0, 0.0, 0.0},/*6*/
            { 0.0, 35.0,  0.0, 21.0,  0.0,  7.0, 0.0, 1.0, 0.0, 0.0},/*7*/
            {35.0,  0.0, 56.0,  0.0, 28.0,  0.0, 8.0, 0.0, 1.0, 0.0},/*8*/
            { 0.0,126.0,  0.0, 84.0,  0.0, 36.0, 0.0, 9.0, 0.0, 1.0}};/*9*/  
  int i, j, k;
  double d=1.0;
  
  printf("Powers of Tn(x) for exp(x) \n");
  for(k=0; k<10; k++) /* power */
  {
    for(j=0; j<=k; j++) TTexp[k][j] = 0.0;
    for(i=0; i<=k; i++) /* loop through c[k][i] for T_i */
    {
      for(j=0; j<=i; j++)
      {
        if(j<=1) d=1.0;
        TTexp[k][j] = TTexp[k][j] + expp[k]*c[i][j]*TP[i][j]/d;
        d = d*2.0;
      }
    }
  }
  printf("telescoping Chebyshev for exp(x) series\n");
  for(k=0; k<10; k++)
  {
    for(j=0; j<=k; j++) printf("TTexp[%d][%d]=%g \n", k, j, TTexp[k][j]);
  }
  printf("\n");
} /* end make_cheb */

/* for integrating exp(x)*Tn(x)/sqrt(1-x^2) using GauLeg */
static void gaulegf(double x1, double x2, double x[], double w[], int n)
{
  int i, j, m;
  double eps = 3.0E-15;
  double p1, p2, p3, pp, xl, xm, z, z1;
  
  m = (n+1)/2;
  xm = 0.5*(x2+x1);
  xl = 0.5*(x2-x1);
  for(i=1; i<=m; i++)
  {
    z = cos(3.141592654*((double)i-0.25)/((double)n+0.5));
    while(1)
    {
      p1 = 1.0;
      p2 = 0.0;
      for(j=1; j<=n; j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*(double)j-1.0)*z*p2-((double)j-1.0)*p3)/
              (double)j;
      }
      pp = (double)n*(z*p1-p2)/(z*z-1.0);
      z1 = z;
      z = z1 - p1/pp;
      if(abs(z-z1) <= eps) break;
    }
    x[i] = xm - xl*z;
    x[n+1-i] = xm + xl*z;
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
} /* end gaulegf */


int main(int argc, char *argv[])
{
  int i, j;
  double sum, a, b, xx, h;
  double x[100], w[100];

  printf("test chebyshev.c on interval -1.0 to 1.0  ordinates, weights\n");
  make_expp(17);
  h = 2.0/10.0;
  for(j=2; j<17; j++)
  {
    printf("evaluate expp(x) n=%d \n", j);
    for(i=0; i<11; i++)
    {
      xx = -1.0+h*(double)i;
      sum = eval(j, expp, xx);
      printf("expp[%d](%g)=%g, err=%g\n", j, xx, sum, sum-exp(xx));  
    }
  }
  printf("\n");
  
  make_TP(17);
  printf("Check roots of Tn(x)\n");
  for(j=2; j<17; j++)
  {
    for(i=0; i<=j; i++) P[i] = TP[j][i];
    for(i=1; i<=j; i++)
    {
      xx = cos((double)(2*i-1)*Pi/(double)(2*j));
      sum = eval(j, P, xx);
      printf("TP[%d][%d](%g)=%g zero?\n", j, i, xx, sum);    
    }
  }
  printf("\n");
  
  make_cheb(10);
  printf("Check telescoping at roots of TTexp\n");
  for(j=3; j<10; j++)
  {
    for(i=0; i<=j; i++) P[i] = TTexp[j][i];
    for(i=1; i<=j; i++)
    {
      xx = cos((double)(2*i-1)*Pi/(double)(2*j));
      sum = eval(j, P, xx);
      printf("TPexp[%d](%g)=%g, exp(%g)=%g \n",
             j, xx, sum, xx, exp(xx));    
    }
  }
  printf("\n");

  for(i=1; i<=17; i++)
  {
    gaulegf(-1.0, 1.0, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      printf("x[%d]=%21.13E, w[%d]=%21.13E \n", j, x[j], j, w[j]);
      sum = sum + w[j];
    }
    printf("              integral(1.0, -1.0..1.0)=%21.13E", sum);
    printf("\n");
  }

  a = 0.5;
  b = 1.0;
  for(i=2; i<=10; i++)
  {
    gaulegf(a, b, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*sin(x[j]);
    }
    printf("integral (0.5,1.0) sin(x) dx = %21.13E\n", sum);
  }
  printf("-cos(1.0)+cos(0.5) = %21.13E\n", (-cos(1.0)+cos(0.5)));
  printf("Maple says 3.372802560E-001 \n");
  printf("\n");

  a = 0.5;
  b = 5.0;
  for(i=2; i<=10; i++)
  {
    gaulegf(a, b, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*exp(x[j]);
    }
    printf("integral (0.5,5.0) exp(x) dx = %21.13E\n", sum);
  }
  printf("exp(5.0)-exp(0.5) = %21.13E\n", (exp(5.0)-exp(0.5)));
  printf("Maple says 1.467644378E+002 \n");
  printf("\n");

  a = 0.5;
  b = 5.0;
  for(i=2; i<=30; i++)
  {
    gaulegf(0.5, 5.0, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*((pow(pow(x[j],x[j]),x[j])*(x[j]*(log(x[j])+1.0)+x[j]*log(x[j]))));
    }
    printf("integral (0.5,5.0) mess(x) dx = %21.13E\n", sum);
  }
  printf("((5.0**5.0)**5.0)-(0.5**0.5)**0.5 = %21.13E\n",
         (pow(pow(5.0,5.0),5.0)-pow(pow(0.5,0.5),0.5)));
  printf("Maple says 2.980232239E+017 \n");
  printf("\nDone.\n");
  return 0;
} /* end main of chebyshev.c */