/*
*   BISECTION ALGORITHM 2.1
*
*   To find a solution to f(x) = 0 given the continuous function
*   f on the interval [A,B], where f(A) and f(B) have
*   opposite signs:
*
* THis code solves Example 1 on Page 50 of the book.
*
*   INPUT:   endpoints A,B; tolerance TOL;
*            maximum number of iterations NI.
*
*   OUTPUT:  approximate solution p or
*            a message that the algorithm fails.
*/

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <iomanip>

#define ZERO      1.0E-20
#define true      1
#define false     0

static   double absval(double);         /* return absolute value of a number */
static   double Func(double);           /* Compute f(x) */
static   void   OUTPUT(int *); /* Set up output files for saving the answer */

using namespace std;

main()
{
   double A,FA,B,FB,C,P,FP,TOL;  // A, B are end points of the interval.
                                 // C is the length of the half interval.
                                 // P is the middle point of the interval.
                                 // TOL is the tolerance for the minimum of the length of the interval [A,B]
   int I,NI,OK,FLAG;             // I is the counter of steps.
                                 // NI is the maximum number of iterations

   // Initialize the problem
   A = 1.0;
   B = 2.0;
   FA = Func(A);
   FB = Func(B);
   TOL = 1.0e-9; // Tolerance for specifying the minimum of the length of the interval [A, B]
   NI = 50;

   OUTPUT(&FLAG);

   /* STEP 1 */
   I = 1;     
   /* STEP 2 */ 
   OK = true; 
   while ((I<=NI) && OK) 
   {
       /* STEP 3 */ 
       /* compute P at the I's step */     
       C = (B - A) / 2.0;
       P = A + C;

       /* STEP 4 */ 
       FP = Func(P);  

       /* print out all intermeditate approximations */
       if (FLAG == 2) cout<<setw(3)<<I <<"       " <<
           setiosflags(ios::fixed) << setprecision(8) << P <<"       "  <<
           setiosflags(ios::fixed) << setprecision(8) <<FP<<endl; 
   
       /* Check if meets the stopping criteria */
       if ((absval(FP)<ZERO) || (C<TOL)) 
       {
	    /* procedure completed successfully */
	    cout<<"\nApproximate solution P = " << setiosflags(ios::fixed) << setprecision(8) << P<<endl;
            cout<< "with F(P) = "<< setiosflags(ios::fixed) << setprecision(8) << FP <<endl;
            cout<<"Number of iterations = "<<I <<endl;
            cout<<" Tolerance = "<< setiosflags(ios::fixed) << setprecision(10) << TOL <<
             " C = "<< setiosflags(ios::fixed) << setprecision(10) << C <<endl;
	    OK = false;
       }
       else 
       {
	    /* STEP 5 */ 
	    I++; 
	    
	    /* STEP 6 */ 
	    /* compute A and B at the I's step */ 
	    if ((FA*FP) > 0.0) 
            {
	       A = P; FA = FP;
	    }
	    else 
            {
	       B = P; FB = FP;
	    }
       }
   }

   if (OK) 
   {
       /* STEP 7 */ 
       /* procedure completed unsuccessfully */
       cout<<"\nIteration number "<< NI <<endl;
       cout<<" gave approximation "<< setiosflags(ios::fixed) << setprecision(8) << P <<endl;
       cout<<"F(P) = "<<setiosflags(ios::fixed) << setprecision(8)<< FP <<
          " not within tolerance :"<< setiosflags(ios::fixed) << setprecision(8)<< TOL<<endl;
   }

   return 0;
}

/* Change function F for a new problem */
double Func(double X)
{
   double f; 

   f = ( X + 4.0 ) * X * X - 10.0;  // Use nested formulation.
   return f;
}

void OUTPUT(int *FLAG)
{
   cout<<"Bisection Method"<<endl;
   cout<<"Select amount of output\n";
   cout<<"1. Answer only\n";
   cout<<"2. All intermeditate approximations\n";
   cout<<"Enter 1 or 2: ";
   scanf("%d", FLAG);
   if (*FLAG == 2) cout<<"  I                 P              F(P)\n";
}   
     

/* Absolute Value Function */
double absval(double val)
{
   if (val >= 0) return val;
   else return -val;
}