/***************************************************
*   Program to demonstrate Bisection subroutine    *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*                C++ version 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                      *
*                                                  *
* The calculated zero is X = -1.000000             *
* The associated Y value is Y = -0.000010          *
* The number of steps was: 23                      *
*                                                  *
***************************************************/
#include <stdio.h>
#include <math.h>

double  e,x,x0,x1;
int     m;


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


/**********************************************
*        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()  {
  double y0,yy;
  m=0;
e100: y0=Y(x0);
  x=(x0+x1)/2; yy=Y(x);
  m++;
  if (yy*y0==0) return;
  if (yy*y0<0) x1=x;
  if (yy*y0>0) x0=x;
  if (fabs(x1-x0)>e) goto e100;
}


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

  Bisect();   // Call bisection routine

  printf("\n\n The calculated zero is X = %f\n",x);
  printf("\n The associated Y value is Y = %f\n",Y(x));
  printf("\n The number of steps was: %d\n\n",m);
}


// End of file Bisect.cpp