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