/* poly_arith_test.c  test poly_arith.h .c functions */
/* poly_arith.c  polynomial arithmetic, include file or use
 * poly_arith.h
 *
 * a polynomial is stored with the constant coefficient in p[0]
 * and the n th order constant coefficient in p[np]
 *
 * example 1 + 2*x + 4*x^2 + 3*x^3 - x^4  has  np = 4 and
 * p[0] = 1;  p[1] = 2;  p[2] = 4; p[3] = 3; p[4] = -1; 
 */

#include <stdio.h>
/* #include "poly_arith.h" */
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

double peval(int n, double x, double c[]);
void add(int np, double p[], int nq, double q[], int *nr, double r[]);
void sub(int np, double p[], int nq, double q[], int *nr, double r[]);
void mul(int np, double p[], int nq, double q[], int *nr, double r[]);
void div(int np, double p[], int nd, double d[],
         int *nq, double q[], int *nr, double r[]);
void norm(int *np, double p[]);
void diff(int np, double p[], int *nr, double r[]);
void integ(int np, double p[], int *nr, double r[]);

int main(int argc, char *argv[])
{
  double p[10];
  double q[10];
  double r[10];
  double d[10];
  double x, val;
  int np, nq, nr, nd;
  int i;

  printf("poly_arith_test.c running \n");
  p[0] = 2.0;
  p[1] = 3.0;
  p[2] = 5.0;
  p[3] = 4.0;
  np = 3;
  q[0] = 1.0;
  q[1] = 2.0;
  q[2] = 4.0;
  q[3] = 99.0;
  nq = 2;

  x = 2.0;
  val = peval(np, x, p);
  printf("test peval(%d, %g, p)=%g \n", np, x, val );
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n\n",
         p[0], p[1], p[2], p[3]);

  add(np, p, nq, q, &nr, r);
  printf("test add(np=%d, p, nq=%d, q, nr=%d, r) \n", np, nq, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("q[0]=%g, q[1]=%g, q[2]=%g \n",
         q[0], q[1], q[2]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g \n\n",
         r[0], r[1], r[2], r[3]);

  sub(np, p, nq, q, &nr, r);
  printf("test sub(np=%d, p, nq=%d, q, nr=%d, r) \n", np, nq, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("q[0]=%g, q[1]=%g, q[2]=%g \n",
         q[0], q[1], q[2]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g \n\n",
         r[0], r[1], r[2], r[3]);

  sub(nq, q, np, p, &nr, r);
  printf("test sub(nq=%d, q, np=%d, p, nr=%d, r) \n", np, nq, nr);
  printf("q[0]=%g, q[1]=%g, q[2]=%g \n",
         q[0], q[1], q[2]);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g \n\n",
         r[0], r[1], r[2], r[3]);

  mul(np, p, nq, q, &nr, r);
  printf("test mul(np=%d, p, nq=%d, q, nr=%d, r) \n", np, nq, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("q[0]=%g, q[1]=%g, q[2]=%g \n",
         q[0], q[1], q[2]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g, r[4]=%g, r[5]=%g \n\n",
         r[0], r[1], r[2], r[3], r[4], r[5]);

  diff(np, p, &nr, r);
  printf("test diff(np=%d, p, nr=%d, r) \n", np, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g \n\n",
         r[0], r[1], r[2]);

  integ(np, p, &nr, r);
  printf("test integ(np=%d, p, nr=%d, r) \n", np, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g \n",
         p[0], p[1], p[2], p[3]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g, r[4]=%g \n\n",
         r[0], r[1], r[2], r[3], r[4]);

  np = 4;
  p[4] = 8.0;
  nd = 2;
  d[0] = 4.0;
  d[1] = 3.0;
  d[2] = 2.0;
  div(np, p, nd, d, &nq, q, &nr, r);
  printf("test div(np=%d, p, nd=%d, d, nq=%d, q, nr=%d, r) \n",
          np, nd, nq, nr);
  printf("p[0]=%g, p[1]=%g, p[2]=%g, p[3]=%g, p[4]=%g \n",
         p[0], p[1], p[2], p[3], p[4]);
  printf("d[0]=%g, d[1]=%g, d[2]=%g \n",
         d[0], d[1], d[2]);
  printf("q[0]=%g, q[1]=%g, q[2]=%g \n",
         q[0], q[1], q[2]);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g, r[4]=%g \n\n",
         r[0], r[1], r[2], r[3], r[4]);

  printf("test norm nr=%d \n", nr);
  norm(&nr, r);
  printf("r[0]=%g, r[1]=%g, r[2]=%g, r[3]=%g, r[4]=%g \n",
         r[0], r[1], r[2], r[3], r[4]);
  printf("nr=%d \n\n", nr);

  printf("poly_arith_test.c finished \n");
  return 0;
}

/* peval  Horners method for evaluating a polynomial */
double peval(int n, double x, double c[])
{ /* 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 */

/* add  polynomial p to polynomial q and place result in polynomial r */
void add(int np, double p[], int nq, double q[], int *nr, double r[])
{
  int i;
  *nr = max(np,nq);
  for(i=0; i<=*nr; i++)
    if(i>np)      r[i] = q[i];
    else if(i>nq) r[i] = p[i];
    else          r[i] = p[i] + q[i];
}

/* subtract polynomial q from polynomial p and place result in polynomial r */
void sub(int np, double p[], int nq, double q[], int *nr, double r[])
{
  int i;
  *nr = max(np,nq);
  for(i=0; i<=*nr; i++)
    if(i>np)      r[i] = -q[i];
    else if(i>nq) r[i] = p[i];
    else          r[i] = p[i] - q[i];
}

/* multiply polynomial p by polynomial q and place result in polynomial r */
void mul(int np, double p[], int nq, double q[], int *nr, double r[])
{
  int i, j;
  *nr = max(np,nq);
  for(i=0; i<=*nr; i++)
    if(i>np)      r[i] = -q[i];
    else if(i>nq) r[i] = p[i];
    else          r[i] = p[i] - q[i];
  *nr = np+nq;
  for(i=0; i<=*nr; i++) r[i] = 0.0;
  for(i=0; i<=np; i++)
    for(j=0; j<=nq; j++) r[i+j] += p[i]*q[j];
}

/* divide polynomial p by polynomial d (np>nd) and place 
 * quotient in polynomial q and remainder in polynomial r */
void div(int np, double p[], int nd, double d[],
         int *nq, double q[], int *nr, double r[])
{
  int i, j, k;
  *nr = np;
  *nq = np-nd;
  k = np;
  for(j=0; j<=np; j++) r[j] = p[j];
  for(i=*nq; i>=0; i--)
  {
    q[i] = r[k]/d[nd];
    for(j=nd; j>=0; j--)
      r[k-nd+j] = r[k-nd+j] - q[i]*d[j];
    k--;
  }
}

/* normalize size, np, such that p[np] not zero */
void norm(int *np, double p[])
{
  while(p[*np] == 0.0 && *np > 0) *np = *np -1;
}

/* differentiate polynomial p and place result in polynomial r */
void diff(int np, double p[], int *nr, double r[])
{
  int i;
  *nr = np-1;
  for(i=0; i<=*nr; i++) r[i] = (double)(i+1)*p[i+1];
}

/* integrate polynomial p and place result in polynomial r */
void integ(int np, double p[], int *nr, double r[])
{
  int i;
  *nr = np+1;
  r[0] = 0.0; /* a reasonable choice */
  for(i=1; i<=*nr; i++) r[i] = p[i-1]/(double)i;
}