/* Return the value of the Chebychev polynomial Tn of degree n at x
   Degree must not be negative
   x must be between -1 and 1
   if an error is detected in the input 2.0 is returned
   if no memory is available 3.0 is returned

   Tn(x) = cos(n*arc cos(x))

   E. Lenz
   Johann-Fichte-Strasse 11
   8 Munich 40
   Germany                                                        */


extern void *calloc();
extern double fabs();

float Chebyshev(n,x)
int n;
float x;
{ int i,j;
  float t;
  float *T,*T1,*T2; /* the coefficient space is allocated dynamically */

  if (fabs(x)>1.0) return(2.0);     /* value out of range */
  if (n<0) return(2.0);             /* negative index not allowed */
  if (n==0) return(1.0);            /* P0(x)=1 */
  if (n==1) return(x);              /* P1(x)=x */
  T=calloc(3*(n+1),sizeof(float));  /* allocate an set to zero */
  if (T==0L) return(3.0);           /* No memory available */
  T1=T+n+1;
  T2=T1+n+1;                        /* set addresses */
  T[0]=T1[1]=1.0;                   /* P0(x)=1 P1(x)=x */
  for (i=1;i<n;i++)
     for (j=i+1;j>-1;j--)
        { T2[j]=2*T1[j-1]-T[j];     /* next generation */
          T[j]=T1[j];               /* make old generation new */
          T1[j]=T2[j];
        }
#ifdef DEBUG
  printf("Chebyshev Tn of the first kind\n");  /* define DEBUG to see the coefficients */
  for (i=n;i>=0;i--) printf("x%d = %f  ",i,T2[i]);
  printf("\n");
#endif
  for (i=n,t=0;i>=0;i--) /* calculate the polynomial */
      {  t*=x;           /* by means of the Horner scheme */
         t+=T2[i];
      }
  free(T);
  return(t);
}