/************Warren D. Smith Mar 2002***************************
If you want to find all real roots, in sorted order,
  of a polynomial P(x) of low degree n (say n<10), then do this:

1. find roots in sorted order of P'(x)  (derivative).
   as well as a strict lowerbound and upperbound on roots of P(x).
2. If these values are x0,x1,...,xn, in sorted order, then consider
each interval [x[i], x[i+1]).  Such an interval
must contain either exactly 0 or exactly 1 root of P(x)
(can tell which by finding sign(P(x[i]) and sign(P(x[i+1])))
which may be found via regula-falsi or related interval-reducing
methods.  Note this method always converges to all roots
and generates them in sorted order. Also note it should
exhibit excellent numerical behavior. For large n it
will be slow, but if n<10 or so, it should be good.
It will also find as a side effect all real roots of all 
derivatives of P(x)...

To compile:  gcc polynom.c -O6 -lm -lc -o polynom
****************************************************/

#define TOOBIGDEG 10

/* a[0] + x*a[1] + ... + x^(n-1) * a[n-1] + x^n * a[n]. */
double EvalPoly(double a[], int n, double x){
  int h;  double y;
  h = n;
  y = a[h];
  for(h--; h>=0; h--){ y *= x; y += a[h]; }
  return y;
}

/* b[0] = a[1], b[1] = 2*a[2], ... b[n-1] = n*a[n]. */
void TakeDeriv(double a[], int n, double b[]){
  int h;  
  for(h=n; h>0; h--){ b[h-1] = a[h]*h; }
}

/* 1 + (|a[0]|+|a[1]|+...+|a[n-1]|)/|a[n]|. */
double StrictUpperBoundOnRootMagnitude(double a[], int n){
  int h;  double y;
  y = 0.0;
  for(h=n-1; h>=0; h--){ y += fabs(a[h]); }
  y /= fabs(a[n]);
  y += 1.0;
  return y;   
}

/* if poly P(x) assumes values flo, fhi of opposite signs at x=xlo, xhi, 
 * then finds root of P(x) in interval lo<x<hi. */
double LocateRoot(double a[], int n, double xlo, double xhi, 
                                     double flo, double fhi){
  double xinterp, xmid, finterp, fmid;
  for(;;){
    xinterp = (xlo*fhi - xhi*flo)/(fhi-flo);
    if(xinterp<=xlo || xhi<=xinterp) return(xinterp);
    finterp = EvalPoly( a, n, xinterp );
    if( finterp>0.0 ){
      if( flo<0.0 ){ xhi = xinterp; fhi = finterp; }
      if( fhi<0.0 ){ xlo = xinterp; flo = finterp; }
    }else if( finterp<0.0 ){
      if( flo>0.0 ){ xhi = xinterp; fhi = finterp; }
      if( fhi>0.0 ){ xlo = xinterp; flo = finterp; }
    }else return(xinterp);
    xmid = (xlo+xhi)*0.5;
    if(xmid<=xlo || xhi<=xmid) return(xmid);
    fmid = EvalPoly( a, n, xmid );
    if( fmid>0.0 ){
      if( flo<0.0 ){ xhi = xmid; fhi = fmid; }
      if( fhi<0.0 ){ xlo = xmid; flo = fmid; }
    }else if( fmid<0.0 ){
      if( flo>0.0 ){ xhi = xmid; fhi = fmid; }
      if( fhi>0.0 ){ xlo = xmid; flo = fmid; }
    }else return(xmid);
  }
}

/* returns true degree of P(x) */
int FindAllRealRoots(double a[], int n, 
        int NumRoots[TOOBIGDEG],
        double root[TOOBIGDEG][TOOBIGDEG]){
  int deg, k, h, ct;
  double lo, hi, flo, fhi, dis, u, coeff[TOOBIGDEG][TOOBIGDEG];
  /* find true degree, which could be <n: */
  for(deg=n; deg>=0; deg--){ if(a[deg] != 0.0) break; }
  /* find coeffs of poly and all derivatives: */
  for(h=deg; h>=0; h--){ coeff[deg][h] = a[h]; }
  for(h=deg; h>=0; h--){ TakeDeriv(coeff[h], h, coeff[h-1]); }
  /* find roots of linear poly, i.e. of (deg-1)th deriv: */
  NumRoots[0]=0; NumRoots[1] = 1;
  root[1][0] = -coeff[1][0] / coeff[1][1];
  /* find roots of quadratic poly, i.e. of (deg-2)th deriv: */
  dis = coeff[2][1]; dis *= dis; dis -= 4.0*coeff[2][0]*coeff[2][2];
  if(dis < 0.0){ NumRoots[2] = 0; }
  else if(dis==0.0){ 
    NumRoots[2] = 1; 
    dis = coeff[2][1];
    if(dis != 0.0){ root[2][0] = -0.5*dis/coeff[2][0]; }
    else{ root[2][0] = 0.0; }
  }else{ 
    NumRoots[2] = 2; 
    if( coeff[2][1] > 0.0 ){ /*avoid cancellation error*/
      root[2][0] = -0.5*(sqrt(dis) + coeff[2][1])/coeff[2][2];
    }else{
      root[2][0] =  0.5*(sqrt(dis) - coeff[2][1])/coeff[2][2];
    }
    root[2][1] = coeff[2][0] / (coeff[2][2]*root[2][0]);
    if( root[2][0] > root[2][1] ){ /*sort*/
      dis = root[2][0]; root[2][0] = root[2][1]; root[2][1] = dis;
    }
  }
  /* find roots of successive polys of degrees 3,4,...,deg: */
  for(h=3; h<=deg; h++){
     u = StrictUpperBoundOnRootMagnitude(coeff[h], h);
     hi = -u;
     fhi = EvalPoly( coeff[h], h, hi );
     ct = 0;
     for(k=0; k<NumRoots[h-1]; k++){
       lo = hi; 
       flo = fhi;
       hi = root[h-1][k];
       fhi = EvalPoly( coeff[h], h, hi );
       if( (flo<0.0 && fhi>0.0) || (flo>0.0 && fhi<0.0) ){
	 root[h][ct] = LocateRoot(coeff[h], h, lo, hi, flo, fhi); ct++;
       }else if( flo==0.0 ){ root[h][ct] = lo; ct++; }
     }
     lo = hi; 
     flo = fhi;
     hi = u;
     fhi = EvalPoly( coeff[h], h, hi );
     if( (flo<0.0 && fhi>0.0) || (flo>0.0 && fhi<0.0) ){
       root[h][ct] = LocateRoot(coeff[h], h, lo, hi, flo, fhi); ct++;
     }else if( flo==0.0 ){ root[h][ct] = lo; ct++; }
     NumRoots[h] = ct;
  }  
  return(deg);
}

/* Finds the degree-n polynomial P(x) so that P(x[i])=a[i] for i=0..n.
 * a[0..n] is overwritten by the coefficients of P(x). */
void PolyInterpolate(double a[], int n, double x[]){
  int i,k;
  for(k=0; k<n; k++){
    for(i=n; i>k; i--){ a[i] = (a[i]-a[i-1])/(x[i]-x[i-k-1]); }
  }
  for(k=n-1; k>=0; k--){
    for(i=k; i<n; i++){ a[i] -= a[i+1]*x[k]; }
  }
}

void main(){ /*test driver*/
  int deg,i,j;
  char which[TOOBIGDEG];
  double a[TOOBIGDEG], x[TOOBIGDEG];
  int NumRoots[TOOBIGDEG];
  double root[TOOBIGDEG][TOOBIGDEG];
  printf("degree of your poly?\n");
  scanf("%u", &deg);
  if(deg >= TOOBIGDEG){ printf("too large\n"); exit(1); }
  if(deg <= 0){ printf("too small\n"); exit(1); }
  printf("Polynomial from: interpolate data (i) or give explicit coeffs (c)?\n");
  scanf("%s", which);
  if(which[0]=='c'){
    printf(
     "input %d coeffs a[0], a[1]..., a[%d] of poly a[0]+a[1]*x+...+a[%d]*x^%d:\n",
	   deg+1, deg, deg, deg);
    for(i=0; i<=deg; i++){     scanf("%lf", &(a[i]) );    }
  }else if(which[0]=='i'){
    printf("input %d  x,P(x) pairs\n", deg+1);
    for(i=0; i<=deg; i++){     scanf("%lf %lf", &(x[i]), &(a[i]) );    }
    printf("Your x,P(x) pairs were:\n");
    for(j=0; j<=deg; j++){  printf("%g %g\n", x[j], a[j] ); }
    PolyInterpolate(a, deg, x);
  }else{
    printf("must choose i or c\n"); exit(1);
  }
  printf("Your coeffs are:\n");
  for(j=0; j<=deg; j++){  printf("%g ", a[j] ); }
  printf("\nHere are the real roots of P(x) and all derivatives:\n");
  deg = FindAllRealRoots(a, deg, NumRoots, root);
  for(i=0; i<=deg; i++){
     printf("degree=%d  NumRoots=%d  roots=", i, NumRoots[i]);
     for(j=0; j<NumRoots[i]; j++){  printf("%g ", root[i][j]);  }
     printf("\n");
  }
}