/* This program implements a binary search to find the root of a cubic
   polynomial.  The primary function is called recursively to do this.
   A binary search is a search in which you narrow down the possibilities
   by half each time.  In this case, the binary search is done in terms
   of an interval [x_min,x_max] in which there is a root of the cubic
   polynomial x*x*x + a*x*x + b*x + c.  By checking the sign of the 
   polynomial at the midpoint x_mid of the interval, we can check 
   whether the root is in the first half of the interval [x_min,x_mid]
   or in the second half [x_mid,x_max].  */

#include<iostream.h>

double a,b,c, tolerance = 1e-10; /* The coefficients of the polynomial 
                                    a,b,c are global variables.  */

int sign(double);                /* returns the sign of a number */
double root(double, double);     /* This is the primary function, which
                                    is called recursively. The values 
                                    passed into it are the limits
                                    x_min, x_max of the interval in
                                    which we'll find the root, value
                                    returned is the root.
				 */
double cubic(double);            /* calculates the polynomial */
bool isroot(double);             /* checks if a given value is a root */

void main()
{
  double min_range=-10., max_range=10., r;

  cout << "Enter the coefficients a, b, c of a cubic polynomial of the form\n"
    "x*x*x + a*x*x + b*x + c:  ";
  cin >> a >> b >> c;

  /* The following loop expands the range until the signs of the polynomial
     evaluated at max_range and min_range are different.  This is always
     possible for a CUBIC polynomial (or indeed any polynomial of odd
     degree). */

  while (sign(cubic(max_range)) == sign(cubic(min_range))) {
    max_range += 10.;
    min_range -= 10.;
  }

  if (isroot(min_range)) {
    /*    cout << "min_range\n";  */
    r=min_range;
  }
  else if (isroot(max_range)) {
    /*    cout << "max_range\n";  */
    r=max_range;
  }
  else 
    r=root(min_range,max_range);

  cout << "x = " << r << " is a root of the polynomial\n"
    "x*x*x + " << a << "*x*x + " << b << "*x + " << c << ".\n";
}


int sign(double z)  /* returns the sign of z, which is -1 if z is
                       negative, 0 if z is 0, and 1 if z is positive */
{
  if (z<-tolerance) return -1;
  if (z<tolerance) return 0;
  return 1;
}

double root(double x_min, double x_max)
{
  double x_mid;

  if (x_max - x_min < tolerance) /* Since x_max > x_min always, this
                                    checks to see if they're essentially
                                    equal.  This will be true once the 
                                    function has been called recursively
                                    many times. */
    return x_min;
  
  x_mid = (x_min + x_max)/2.;

  if (isroot(x_mid)) 
    return x_mid;
  
  if (sign(cubic(x_mid))!=sign(cubic(x_min)))
    return root(x_min,x_mid);
  else
    return root(x_mid,x_max);
}

double cubic(double x)
{
  return x*x*x + a*x*x + b*x + c;
}

bool isroot(double x)
{ 
  return (sign(cubic(x))==0);
}