/* peval.c Horners method for evaluating a polynomial */
#include <stdio.h>
#include "udrnrt.h"

double peval(int n, double x, double c[])
{ /* Horners method for evaluating a polynomial   */
  /* an nth order polynomial has n+1 coefficients */
  /* stored in c[0] through c[n]                  */
  /* y = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n  */
  int i;
  double y = c[n]*x;
  if(n<=0) return c[0];
  for(i=n-1; i>0; i--) y = (c[i]+y)*x;
  return y+c[0];
} /* end peval */

int main(int argc, char *argv[])
{
  double y, c[11];
  double r[11];
  int i, j, n;

  printf("peval.c running\n");
  /* test peval by generation roots of a polynomial, */
  /* then generating the polynomial, then evaluating */
  for(n=3; n<11; n++) /* n is degree of polynomial */
  {
    printf("\ny = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n, for n=%d\n", n);
    printf("roots are:\n");
    r[0] = 0.1;
    printf("r[0]=%g\n", r[0]);
    r[1] = 2.5;
    printf("r[1]=%g\n", r[1]);
    for(i=2; i<n; i++)
    {
      r[i] = r[i-1]+r[i-2];
      printf("r[%d]=%g\n", i, r[i]);
    }
    /* generate polynomial coefficients */
    c[0] = r[0]*r[1];
    c[1] = -(r[0]+r[1]);
    c[2] = 1.0;
    for(j=2; j<n; j++)
    {
      for(i=j+1; i>0; i--) c[i] = c[i-1];
      c[0] = 0.0;
      for(i=0; i<=j; i++) c[i] = c[i] - r[j]*c[i+1];
    }
    printf("polynomial coefficients are:\n");
    for(i=0; i<=n; i++) printf("c[%d]=%g\n", i, c[i]);
    printf("polynomial evaluated at each root:\n");
    for(i=0; i<n; i++)
    {
      y = peval(n, r[i], c);
      printf("r[%d]=%g, peval(%d,r[%d],c)=%g\n",
             i, r[i], n, i, y);
    }
  }
  printf("\ntest with random numbers 0 < r < 1.0 \n");
  n = 10;
  for(i=0; i<n; i++) r[i] = udrnrt();
  /* generate polynomial coefficients */
  c[0] = r[0]*r[1];
  c[1] = -(r[0]+r[1]);
  c[2] = 1.0;
  for(j=2; j<n; j++)
  {
    for(i=j+1; i>0; i--) c[i] = c[i-1];
    c[0] = 0.0;
    for(i=0; i<=j; i++) c[i] = c[i] - r[j]*c[i+1];
  }
  printf("random polynomial coefficients are:\n");
  for(i=0; i<=n; i++) printf("c[%d]=%g\n", i, c[i]);
  printf("random polynomial evaluated at eachroot:\n");
  for(i=0; i<n; i++)
  {
    y = peval(n, r[i], c);
    printf("r[%d]=%g, peval(%d,r[%d],c)=%g\n",
           i, r[i], n, i, y);
  }
  printf("thus we see the size of the largest root affects error\n");
  return 0;
} /* end peval.c */