/* Return the value of the Gegenbauer polynomial 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
   E. Lenz
   Johann-Fichte-Strasse 11
   8 Munich 40
   Germany                                                        */

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

float Gegenbauer(n,x,a)
int n;
float x,a;
{ int i,j;
  float t;
  float *G,*G1,*G2; /* 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)
     if (a==0)
       return(2*x);
     else
       return(2*a*x);               /* P1(x)=2*a*x */
  G=calloc(3*(n+1),sizeof(float));  /* allocate an set to zero */
  if (G==0L) return(3.0);           /* No memory available */
  G1=G+n+1;
  G2=G1+n+1;                        /* set addresses */
  G[0]=1.0;                         /* P0(x)=1 */
  if (a==0)
      G1[1]=2;
  else
      G1[1]=2*a;                    /* P1(x)=2*a*x */
  for (i=1;i<n;i++)
     for (j=i+1;j>-1;j--)
        { G2[j]=(2*(i+a)*G1[j-1]-(i+2*a-1)*G[j])/(i+1);  /* next generation */
          G[j]=G1[j];               /* make old generation new */
          G1[j]=G2[j];
        }
#ifdef DEBUG
  printf("Gegenbauer(%f)\n",a);  /* define DEBUG to see the coefficients */
  for (i=n;i>=0;i--) printf("x%d = %f  ",i,G2[i]);
  printf("\n");
#endif
  for (i=n,t=0;i>=0;i--) /* calculate the polynomial */
      {  t*=x;           /* by means of the Horner scheme */
         t+=G2[i];
      }
  free(G);
  return(t);
}