/****************************************************************
*   Program to demonstrate Bisection & Quartile subroutines     *
* ------------------------------------------------------------- *
* References:                                                   *
*                                                               *
*  Bisection: BASIC Scientific Subroutines, Vol. II             *
*             By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981.     *
*                                                               *
*  Quartile:  After an algorithm provided By Namir Shammas.     *
*                                                               *
*                             C++ Release By J-P Moreau, Paris. *
*                                    (www.jpmoreau.fr)          *
* ------------------------------------------------------------- *
* Example:  Find a real root of f(x)=(x+1)^5                    *
*                                                               *
* Sample run:                                                   *
*                                                               *
* What is the initial range (X0,X1):                            *
*    X0 = -5                                                    *
*    X1 = 0                                                     *
* Convergence criterion: 1e-6                                   *
*                                                               *
*                                                               *
* Bisection: The calculated zero is X = -9.999996e-001          *
*                                                               *
* The associated Y value is Y =  5.850012e-033                  *
*                                                               *
* The number of steps was: 23                                   *
*                                                               *
*                                                               *
* Quartile: The calculated zero is X = -1.000000e+000           *
*                                                               *
* The associated Y value is Y = -5.712903e-036                  *
*                                                               *
* The number of steps was: 12                                   *
*                                                               *
*****************************************************************
Note: The quartile method is slightly faster and more accurate
      than the classical Bisection Method.
---------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

  double a,b,root, e,x,x0,x1;
  int    m;


/**********************************************/
double Y(double x) {
  return (x+1)*(x+1)*(x+1)*(x+1)*(x+1);
}
/**********************************************/


/**********************************************
*        Bisection method subroutine          *
* ------------------------------------------- *
* This routine iteratively seeks the zero of  *
* function Y(x) using the method of interval  *
* halving until the interval is less than e   *
* in width. It is assumed that the function   *
* Y(x) is available from a function routine.  *
* ------------------------------------------- *
* Input values: range (x0,x1), and e.         *
* Output values: root x, Y(x) and number of   *
* steps m.                                    *
**********************************************/
void Bisect(void) {
//Labels: e10
double y0, yy;

  m=0;
  if (Y(x0) * Y(x1) > 0.0) {
    printf(" Bisection: No root found in given interval.\n");
    return;
  }
  else {
e10: y0=Y(x0);
    x=(x0+x1)/2; yy=Y(x);
    m++;
    if (yy*y0==0.0) return;
    if (yy*y0<0.0)  x1=x;
    if (yy*y0>0.0)  x0=x;
    if (fabs(x1-x0)>e) goto e10;
  }
}

/***********************************************
*         Quartile method subroutine           *
* -------------------------------------------- *
* This routine iteratively seeks the zero of   *
* function Y(x) using the method of interval   *
* quarting until the interval is less than tol *
* in width. It is assumed that the function    *
* Y(x) is available from a function routine.   *
* -------------------------------------------- *
* Input values: range (a, b), and tol.         *
* Output values: found root and number of used *
* steps.                                       *
***********************************************/
void Quartile(double a,double b,double tol, double *root, int *step) {
double m, co;

  co=0.25; //can be 0.3
  *step=0;
  if (Y(a) * Y(b) > 0.0) {
    printf(" Quartile: No root found in given interval.\n");
    return;
  }
  else {
	while (fabs(a-b) > tol) {
      if (fabs(Y(a)) < fabs(Y(b)))
        m = a + co*(b-a);
      else
        m = a + (1.0-co)*(b-a);
      if (Y(m)*Y(a) > 0.0)  a= m; else b=m;
      (*step)=(*step)+1;
    }
    *root = (a+b)/2.0;
  }
}


void main() {

  printf("\n What is the initial range (X0,X1):\n\n");
  printf("    X0 = "); scanf("%lf", &x0);
  printf("    X1 = "); scanf("%lf", &x1);
  printf("\n Convergence criterion: "); scanf("%lf", &e);
  printf("\n");

  a=x0; b=x1;

  Bisect();   // Call Bisection routine

  if (m!=0) {
    printf("\n Bisection: The calculated zero is X = %e\n\n", x);
    printf(" The associated Y value is Y = %e\n\n", Y(x));
    printf(" The number of steps was: %d\n\n", m);
  }

  Quartile(a, b, e, &root, &m); //Call Quartile routine

  if (m!=0) {
    printf("\n Quartile: The calculated zero is X = %e\n\n", root);
    printf(" The associated Y value is Y = %e\n\n", Y(root));
    printf(" The number of steps was: %d\n\n", m);
  }

}

//End of file tquart.cpp