// test_poly.c   polyprt  polyval  polyderiv  polyintegrate  polyfit
//               polynomial stored  p[0] + p[1]*x + ...p[pwr]x^pwr

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "poly.h"

int main(int argc, char * argv[])
{
  int pwr = 3;
  int n = 10;
  double p[] = {1.2, 2.4, 3.6, 4.8}; // changed for later test
  double pck[10];
  double pqck[10];
  double pre[10];
  double pderiv[3];
  double pintgr[5];
  double x[10];
  double y[10];
  double yc[10];
  double err;
  double pa[] = {3.0, 5.0, 8.0, 12.0};
  int na = 3;
  double pb[] = {1.0, 2.0, 3.0, 4.0};
  int nb = 3;
  double xr[] = {1.0, 2.0, 3.0, 4.0};
  int nrt = 4; 
  double result[10];
  int nr[1]; // returned pwr
  int nq[1];
  int npre[1];

  printf("test_poly.c running \n");
  printf("test polyval\n");
  for(int i=0; i<n; i++)
  {
    x[i] = (double)(i+1)*0.25;
  }
  for(int i=0; i<n; i++)
  {
    y[i] = polyval(pwr, p, x[i]);
    yc[i] = y[i];
    printf("i=%d, x[i]=%f, y[i]=%f \n", i, x[i], y[i]);
  }

  printf("computed values\n");
  for(int j=0; j<=pwr; j++) p[j] = -1.0;
  polyfit(pwr, n, x, y, p);
  for(int j=0; j<=pwr; j++) printf("p[%d]=%f \n", j, p[j]);

  for(int i=0; i<n; i++)
  {
    y[i] = polyval(pwr, p, x[i]);
    err = abs(y[i]-yc[i]);
    printf("i=%d, x[i]=%f, y[i]=%f, err=%f \n", i, x[i], y[i], err);
  }

  printf("\n");
  printf("test polycopy\n");
  polyprt(nb, pb);
  polycopy(nb, pb, result);
  printf("copied = \n");
  polyprt(nb, result);

  printf("\n");
  printf("test polyadd\n");
  polyprt(na, pa);
  printf("+\n");
  polyprt(nb, pb);
  printf("sum =\n");
  polyadd(na, pa, nb, pb, nr, result);
  polyprt(nr[0], result);

  printf("\n");
  printf("test polysub\n");
  polyprt(na, pa);
  printf("-\n");
  polyprt(nb, pb);
  printf("difference =\n");
  polysub(na, pa, nb, pb, nr, result);
  polyprt(nr[0], result);

  printf("\n");
  printf("test polymul\n");
  polyprt(na, pa);
  printf("*\n");
  polyprt(nb, pb);
  polymul(na, pa, nb, pb, nr, result);
  printf("product =\n");
  polyprt(nr[0], result);

  printf("\n");
  printf("test polydiv\n");
  polyprt(nr[0], result);
  printf("/\n");
  polyprt(nb, pb);
  polydiv(nr[0], result, nb, pb, nq, pqck, nr, pck);
  printf("quotient =\n");
  polyprt(nq[0], pqck);
  printf("remainder =\n");
  polyprt(nr[0], pck);

  printf("\n");
  printf("test polyderiv\n");
  p[0] = 4.0;
  p[1] = 3.0;
  p[2] = 2.0;
  p[3] = 1.0;

  polyprt(3, p);
  polyderiv(3, p, pderiv);
  printf("derivative\n");
  polyprt(2, pderiv);
  printf("integrate\n");
  polyintegrate(2, pderiv, pck);
  polyprt(3, pck);

  printf("\n");
  printf("test polyintegrate\n");
  polyprt(3, p);
  polyintegrate(3, p, pintgr);
  printf("integral\n");
  polyprt(4, pintgr);
  printf("derivative\n");
  polyderiv(4, pintgr, pck);
  polyprt(3, pck);

  printf("\n");
  printf("test polymake from roots xr\n");
  for(int i=0; i<1; i++) // nrt=1
    printf("xr[%d]=%f \n", i, xr[i]);
  polymake(1, xr, nr, result);
  printf("polynomial nr=%d  p =\n", nr[0]);
  polyprt(nr[0], result);
  printf("check 0 %f \n",polyval(nr[0], result, xr[0]));

  printf("test polymake from roots xr\n");
  for(int i=0; i<2; i++)
    printf("xr[%d]=%f \n", i, xr[i]);
  polymake(2, xr, nr, result);
  printf("polynomial nr=%d p=\n", nr[0]);
  polyprt(nr[0], result);
  printf("check 0 \n", polyval(nr[0], result, xr[0]));
  printf("check 0 \n", polyval(nr[0], result, xr[1]));

  printf("test polymake from roots xr\n");
  for(int i=0; i<3; i++)
    printf("xr[%d]=%f \n", i, xr[i]);
  polymake(3, xr, nr, result);
  printf("polynomial nr=%d p=\n", nr[0]);
  polyprt(nr[0], result);
  printf("check 0 %f \n", polyval(nr[0], result, xr[0]));
  printf("check 0 %f \n", polyval(nr[0], result, xr[1]));
  printf("check 0 %f \n", polyval(nr[0], result, xr[2]));

  printf("test polymake from roots xr\n");
  for(int i=0; i<nrt; i++)
    printf("xr[%d]=%f \n", i, xr[i]);
  polymake(nrt, xr, nr, result);
  printf("polynomial =\n");
  polyprt(nr[0], result);
  printf("check 0 %f\n", polyval(nr[0], result, xr[0]));
  printf("check 0 %f\n", polyval(nr[0], result, xr[1]));
  printf("check 0 %f\n", polyval(nr[0], result, xr[2]));
  printf("check 0 %f\n", polyval(nr[0], result, xr[3]));
  polycopy(nr[0], result, pck); // made from roots 1.0, 2.0, 3.0, 4.0

  printf("\n");
  printf("test polyremovearoot 4.0 from pre \n");
  polycopy(nr[0], result, pre); // made from roots 1.0, 2.0, 3.0, 4.0
  npre[0] = nr[0];
  x[0] = 4.0;
  polyprt(npre[0], pre);
  polyremovearoot(x[0], npre, pre);
  printf("with root 4.0 removed = \n");
  polyprt(npre[0], pre); // npre reduced

  printf("\n");
  printf("test polyfindaroot \n");
  polyprt(nr[0], pck);
  x[0] = polyfindaroot (nr[0], pck);
  printf("x[0] = %f\n",x[0]);
  printf("remove that root \n");
  polyremovearoot(x[0], nr, pck);

  printf("test polyfindaroot \n");
  polyprt(nr[0], pck);
  x[1] = polyfindaroot (nr[0], pck);
  printf("x[1] = %f\n", x[1]);
  printf("remove that root \n");
  polyremovearoot(x[1], nr, pck);

  printf("test polyfindaroot \n");
  polyprt(nr[0], pck);
  x[2] = polyfindaroot (nr[0], pck);
  printf("x[2] = %f\n",x[2]);
  printf("remove that root \n");
  polyremovearoot(x[2], nr, pck);

  printf("test polyfindaroot \n");
  polyprt(nr[0], pck);
  x[3] = polyfindaroot (nr[0], pck);
  printf("x[3] =%f\n ",x[3]);

  printf("test_poly.c finished \n");
  return 0;
} // end test_poly.c