/*
*   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>

#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(FILE **, 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
   FILE *OUP[1];

   // 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(OUP, &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) fprintf(*OUP,"%3d   %15.8e   %15.7e \n",I,P,FP); 
   
       /* Check if meets the stopping criteria */
       if ((absval(FP)<ZERO) || (C<TOL)) 
       {
	    /* procedure completed successfully */
	    fprintf(*OUP,"\nApproximate solution P = %11.8e \n",P);
	    fprintf(*OUP,"with F(P) = %12.8e\n",FP);
	    fprintf(*OUP,"Number of iterations = %3d",I);
	    fprintf(*OUP,"    Tolerance = %15.8e, C = %15.8e\n",TOL, C);
	    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 */
       fprintf(*OUP,"\nIteration number %3d",NI);
       fprintf(*OUP," gave approximation %12.8f\n",P);
       fprintf(*OUP,"F(P) = %12.8f not within tolerance : %15.8e\n",FP, TOL);
   }
   fclose(*OUP);

   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(FILE **OUP, int *FLAG)
{
   char NAME[30];

   printf("Select output destination\n");
   printf("1. Screen\n");
   printf("2. Text file\n");
   printf("Enter 1 or 2: ");

   scanf("%d", FLAG);
   if (*FLAG == 2) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("A:OUTPUT.DTA\n");
      scanf("%s", NAME);
      *OUP = fopen(NAME, "w");
   }
   else *OUP = stdout;
   fprintf(*OUP,"Bisection Method\n");
   printf("Select amount of output\n");
   printf("1. Answer only\n");
   printf("2. All intermeditate approximations\n");
   printf("Enter 1 or 2: ");
   scanf("%d", FLAG);
   if (*FLAG == 2) fprintf(*OUP, "  I                 P              F(P)\n");
}   
     

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