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