#include <stdio.h>
#include <sys/file.h>
#include <math.h>

#define F (float)
#define D (double)
#define MAXORD 13
#define MAXITER 200
#define EPS (double) 1.0e-20
#define forever for(;;)
int MINDETREND=100;
#define close_enough 0.95

double coef[MAXORD+1];

double evalpoly(), evalpoly2(), evalpoly3();
int newton(), newton2(), polyfit();

main(ac,av)
int ac;
char **av;
{
   int fdi, fdo, nx, nt, xint, first, k, knx, order=5, where, flag, end;
   char ifile[80], ofile[80];
   float *fts;
   double *dts, last, l;
   int check, exval = 0;

   setpar(ac,av);
   if ( ! getpar("ntr","d",&nx)) getpar("nx","d",&nx);
   if (getpar("if","s",ifile)) {
      if ((fdi = open(ifile,O_RDONLY,0644)) <= 2) {
         fprintf(stderr,"%s:  Cannot open %s\n",av[0],ifile);
         exit(-1);
      }
   }
   else fdi = 0;
   if (getpar("of","s",ofile)) {
      if ((fdo = open(ofile,O_WRONLY|O_CREAT|O_TRUNC,0644)) <= 2) {
         fprintf(stderr,"%s:  Cannot open %s\n",av[0],ofile);
         exit(-1);
      }
   }
   else fdo = 1;
   mstpar("nt","d",&nt);
   getpar("order","d",&order);
   getpar("min","d",&MINDETREND);
   endpar();

   if (order > MAXORD) {
      fprintf(stderr,"%s:  using maximum order (%d)\n",av[0],MAXORD);
      order = MAXORD;
   }

   if ((fts = (float *) malloc(sizeof(float)*nt)) == NULL) {
      fprintf(stderr,"%s:  cannot allocate %d bytes\n",av[0],sizeof(float)*nt);
      exit(-1);
   }

   if ((dts = (double *) malloc(sizeof(double)*nt)) == NULL) {
      fprintf(stderr,"%s:  cannot allocate %d bytes\n",av[0],sizeof(double)*nt);
      exit(-1);
   }

   for (knx = 0; knx < nx; knx++) {
      flag = 0;
      if (doread(fdi,fts,sizeof(float)*nt) != sizeof(float)*nt) {
         fprintf(stderr,"%s:  read error\n",av[0]);
         exit(-1);
      }
      for (k = 0; k < nt; k++) dts[k] = D fts[k];
      last = dts[nt - 1];
      for (k = nt - 1; k > 0; k--)
	 if (dts[k] * last <= D 0.0) break;
      if (nt - k < (int) (MINDETREND * close_enough)) {
	 fprintf(stderr,"min flag, trace %d\n",knx);
	 flag = 1;
         last = dts[nt - MINDETREND - 1];
         for (k = nt - MINDETREND - 1; k > 0; k--)
	    if (dts[k] * last <= D 0.0) break;
      }
      first = k + 1;
      if (polyfit(dts+first,order,nt-first) < 0) {
	 fprintf(stderr,"%s:  warning - ill-conditioned, trace %d\n",av[0],knx);
      }
/* old method: */
#ifdef DUM
      where = end = nt-first;
      if (flag) where = end = nt-first-MINDETREND;
      for (k=0; k<MAXITER; k++) {
         if (newton(&xint,order,where) < 0) {
	    /* let x-intercept be turning point */
	    if (newton2(&xint,order,where) >= 0) {
	       if (xint >= end) continue;
	       fprintf(stderr,"%s:  using turning point, trace %d\n",av[0],knx);
	       /* set coef[0] so poly intersects time axis at turning point */
	       /* this does leave some drift */
	       coef[0] = D 0.0;
	       coef[0] = -evalpoly( D xint,order);
            }
         }
         if (xint < end) break;
         where -= 10;
      }
      if (k < MAXITER) {
         xint += first;
         if ((xint < 0) || (xint > nt))
            fprintf(stderr,"%s:  warning - bad fit1, trace %d, int=%d\n",av[0],knx,xint);
         else
            for (k = xint + 1; k < nt; k++) {
	       l = D (k - first);
               dts[k] -= evalpoly(l,order);;
	    }
      }
      else fprintf(stderr,"%s:  warning - bad fit2, trace %d, int=%d\n",av[0],knx,xint);
#endif
/* better method: */
      check = newton(&xint,order,0);
      if ((check < 0) || (xint + first < 0) || (xint + first > nt)) {
         /* let x-intercept be turning point */
         if (newton2(&xint,order,0) >= 0) {
            fprintf(stderr,"%s:  using turning point, trace %d\n",av[0],knx);
            /* set coef[0] so poly intersects time axis at turning point */
            /* this does leave some drift */
            coef[0] = D 0.0;
            coef[0] = -evalpoly( D xint,order);
         }
      }
      xint += first;
      if ((xint < 0) || (xint > nt)) {
         fprintf(stderr,"%s:  warning - bad fit1, trace %d, int=%d\n",av[0],knx,xint);
	 exval = 2;
      }
      else
         for (k = xint + 1; k < nt; k++) {
            l = D (k - first);
            dts[k] -= evalpoly(l,order);;
         }
      for (k = 0; k < nt; k++) fts[k] = F dts[k];
      if (dowrite(fdo,fts,4*nt) != 4*nt) {
         fprintf(stderr,"%s:  write error\n",av[0]);
         exit(-1);
      }
   }
   exit(exval);
}

int
polyfit(y,order,n)
double *y;
int order, n;
{
   static int first_call=1;
   double sumxiy[MAXORD+1], sumxi[2*MAXORD+1];
   double xval, yval;
   static double **a;
   static int *rows;
   int i, j, k, ret_val=0;
   double **init_pivot();
   int solve_pivot();
   void backsub(), reaugment();

   /* initialize sum values */
   for (i = 0; i <=     order; i++) sumxiy[i] = D 0.0;
   for (i = 1; i <= 2 * order; i++)  sumxi[i] = D 0.0;
   sumxi[0] = D n;

   /* fill sum values */
   for (k = 0; k < n; k++) {
      yval = y[k];
      for (i = 0; i <= order; i++) {
         sumxiy[i] += yval;
         yval *= D k;
      }
      xval = D k;
      for (i = 1; i <= 2 * order; i++) {
         sumxi[i] += xval;
         xval *= D k;
      }
   }

   /* get array to fill */
   if (first_call) { a = init_pivot(order+1,&rows); first_call = 0; }

   /* fill array */
   for (i = 0; i <= order; i++) {
      for (j = 0; j <= order; j++) a[i][j] = sumxi[i + j];
      a[i][order+1] = sumxiy[i];
   }

   /* solve for poly */
   if (solve_pivot(a,rows,order+1) < 0) ret_val = -1;
   backsub(a,order+1);
   for (i = 0; i <= order; i++) coef[i] = a[i][order+1];
   return(ret_val);
}

double
evalpoly(x,order)
double x;
int order;
{
   double val;

   val = coef[order];
   while (order) {
      val *= x;
      val += coef[--order];
   }
   return(val);
}

double
evalpoly2(x,order)
double x;
int order;
{
   double val;

   val = D order * coef[order];
   while (order > 1) {
      val *= x;
      --order;
      val += D order * coef[order];
   }
   return(val);
}

double
evalpoly3(x,order)
double x;
int order;
{
   double val;

   val = D (order * (order - 1)) * coef[order];
   while (order > 2) {
      val *= x;
      --order;
      val += D (order * (order - 1)) * coef[order];
   }
   return(val);
}

int
newton(xint,order,n)
int *xint, order, n;
{
   double yval, xval, xold;
   int i;

   i = 0;
   xval = D (n-1);
   yval = evalpoly(xval, order);
   while (fabs(yval) > EPS) {
      xold = xval;
      xval = xold - yval / evalpoly2(xold,order);
      yval = evalpoly(xval,order);
      if (++i > MAXITER) return(-1);
   }
   *xint = (int) xval;
   return(0);
}

int
newton2(xint,order,n)
int *xint, order, n;
{
   double yval, xval, xold;
   int i;

   i = 0;
   xval = D (n-1);
   yval = evalpoly2(xval, order);
   while (fabs(yval) > EPS) {
      xold = xval;
      xval = xold - yval / evalpoly3(xold,order);
      yval = evalpoly2(xval,order);
      if (++i > MAXITER) return(-1);
   }
   *xint = (int) xval;
   return(0);
}

doread(fd,buf,bytes)
int fd, bytes;
char *buf;
{
	int nread, total;

	if ((total = read(fd,buf,bytes)) <= 0) return(total);
	buf += total;
	while (bytes > total) {
		if ((nread = read(fd,buf,bytes-total)) <= 0) return(total);
		total += nread;
		buf += nread;
	}
	return(total);
}

int
dowrite(fd,buf,bytes)
int fd, bytes;
char *buf;
{
	int nwrite, total;

	if ((total = write(fd,buf,bytes)) <= 0) return(total);
	buf += total;
	while (bytes > total) {
		if ((nwrite = write(fd,buf,bytes-total)) <= 0) return(total);
		total += nwrite;
		buf += nwrite;
	}
	return(total);
}