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

double abs_double(double y)
{
  y = abs(y) ;
    return y;
}

double fx_val(double a, double b, double c, double d, double e, double x)
{
  double a_x = a * x * x * x * x ;
  double b_x = b * x * x * x ;
  double c_x = c * x * x ;
  double d_x = d * x ;
  double e_x = e ;
  double sum_x = a_x + b_x + c_x + d_x + e_x ;
    return sum_x;
}

double fx_dval(double a, double b, double c, double d, double e, double x)
{
  double a_dx = 4 * a * x * x * x ;
  double b_dx = 3 * b * x * x ;
  double c_dx = 2 * c * x ;
  double d_dx = d ;
  double e_dx = 0 ;
  double sum_dx = a_dx + b_dx + c_dx + d_dx + e_dx ;
    return sum_dx;
}

double fx_ddval(double a, double b, double c, double d, double e, double x)
{
  double a_ddx = 4 * 3 * a * x * x ;
  double b_ddx = 3 * 2 * b * x ;
  double c_ddx = 2 * 1 * c ;
  double d_ddx = 0 ;
  double e_ddx = 0 ;
  double sum_ddx = a_ddx + b_ddx + c_ddx + d_ddx + e_ddx ;
    return sum_ddx;
}


double newrfind_halley(double a, double b, double c, double d, double e, double x)
{
  double halley_counter = 0;
  double z_change;
  while (halley_counter < 10000)
    {
      double sum_h = fx_val(a,b,c,d,e,x);
      double sum_dh = fx_dval(a,b,c,d,e,x);
      double sum_ddh = fx_ddval(a,b,c,d,e,x);
      double abs_sum_dh = abs_double(sum_dh);
      z_change = (2 * sum_h * sum_dh)/((2 * abs_sum_dh * abs_sum_dh)-(sum_h * sum_ddh));
      x = x - z_change;
	if (z_change < .000001 || z_change > -.000001)
	  {
	    return x;
	  }
	else
	  {
	    halley_counter = halley_counter + 1;
	  }
    }
  return 500000;
}

int rootbound(double a, double b, double c, double d, double e, double r, double l)
{
  double v = 0;
    if (a*(4*a*l+b) < 0 )
      {
	v = v + 1 ;
      }

    if (((4*a*l+b)*(6*a*l*l+3*b*l+c))<0)
      {
	v = v + 1;
      }

    if ((6*a*l*l+3*b*l+c)*(4*a*l*l*l+3*b*l*l+2*c*l+d) < 0)
      {
	v = v + 1;
      }

    if ((4*a*l*l*l+3*b*l*l+2*c*l+d)*(a*l*l*l*l+b*l*l*l+c*l*l+d*l+e) < 0)
      {
	v= v + 1 ;
      }

    if (a*(4*a*r+b) < 0 )
      {
	v = v - 1 ;
      }

    if (((4*a*r+b)*(6*a*r*r+3*b*r+c))<0)
      {
	v = v - 1;
      }

    if ((6*a*r*r+3*b*r+c)*(4*a*r*r*r+3*b*r*r+2*c*r+d) < 0)
      {
	v = v - 1;
      }

    if ((4*a*r*r*r+3*b*r*r+2*c*r+d)*(a*r*r*r*r+b*r*r*r+c*r*r+d*r+e) < 0)
      {
	v= v - 1 ;
      }

return v;
}

int main(int argc, char **argv)
{
	double a, b, c, d, e, l, r;
	FILE *in;

	if (argv[1] == NULL) {
		printf("You need an input file.\n");
		return -1;
	}
	in = fopen(argv[1], "r");
	if (in == NULL)
		return -1;
	fscanf(in, "%lf", &a);
	fscanf(in, "%lf", &b);
	fscanf(in, "%lf", &c);
	fscanf(in, "%lf", &d);
	fscanf(in, "%lf", &e);
	fscanf(in, "%lf", &l);
	fscanf(in, "%lf", &r);

	double og_l = l;
	double og_r = r;
	double x = 0;
	double counter ;
	double rootchecker ;
	r = l + .5 ;
	l = l - .5;
	double boudan_var = rootbound(a,b,c,d,e,l,r);
	 {
	  if (boudan_var == 0)
	    printf("The polynomial has no roots in the given interval.");
	  else
	    {
	      for(counter = og_l ; counter == og_r ; counter = counter + .5)
		{
		  l = l + .5;
		  printf("%f", l);
		  x = l ;
		  rootchecker = newrfind_halley(a,b,c,d,e,x);
		  if (rootchecker != 500000)
		    {
		      printf("Root found: %f", x);
		    }
		  else
		    {
		      printf("No Root Found");
		    }
		}
	    }
	 }

    //The values are read into the variable a, b, c, d, e, l and r.
    //You have to find the bound on the number of roots using rootbound function.
    //If it is > 0 try to find the roots using newrfind function.
    //You may use the fval, fdval and fddval funtions accordingly in implementing the halleys's method.


    fclose(in);

    return 0;
}