/* Return the value of the Hermite polynomial of degree n at x
   Degree must not be negative
   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 Hermite(n,x)
int n;
float x;
{ int i,j;
  float t;
  float *H,*H1,*H2; /* the coefficient space is allocated dynamically */

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