Mathematics Programming

This algorithm uses the bisection algorithm to solve the equation 3x² + 12x + 1 = 0.  The equation being solved can be changed by recoding the F function.

/**************************************************
 * File:        Bisection.c
 * Description: Simple implementation of the
 *              bisection algorithm in C
 **************************************************/

#include <math.h>
#include <stdio.h>
 
int Bisect (double, double, double, unsigned, double *, double *, unsigned *); 
double F (double);
void Input (double *, double *, double *, unsigned *); 
void Output (double, double, unsigned, int); 
void Step (double *, double *, double *, double *, double *, double *); 
 
int main () {
  double X1;
  double X2;
  double Y1;
  double Y2;
  double X;
  double Y;
  double Tol;
  unsigned Limit;
  unsigned N;
  int OK;
  
  printf ("Bisection algorithm - C++\n");

  Input (&X1, &X2, &Tol, &Limit);
  OK = Bisect (X1, X2, Tol, Limit, &X, &Y, &N);
  Output (X, Y, N, OK);
  
  return 0;
}

int Bisect (double X1, double X2, double Tol, unsigned Limit, 
            double * X, double * Y, unsigned * i) {
  int OK;
  double Y1;
  double Y2;
  
  Y1 = F (X1);
  Y2 = F (X2);
  OK = (Y1 > 0) != (Y2 > 0);
  printf ("%2d  %10f  %10f\n", 0, X1, Y1);
  printf ("%2d  %10f  %10f\n", 0, X2, Y2);
  if (OK) {
    *i = 0;
    while ((fabs (X1 - X2) > Tol) && (*i < Limit)) {
      Step (&X1, &Y1, &X2, &Y2, X, Y);
      ++(*i);
      printf ("%2d  %10f  %10f\n", *i, *X, *Y);
    }
  }
  return OK;
} 

void Input (double * X1, double * X2, double * Tol, unsigned * Limit) {
  printf ("First X value?  ");
  scanf ("%lf", X1);
  printf ("Second X value?  ");
  scanf ("%lf", X2);
  printf ("Tolerance?  ");
  scanf ("%lf", Tol);
  printf ("Maximum iterations?  ");
  scanf ("%d", Limit);
}
 
void Output (double X, double Y, unsigned N, int OK) {
  if (OK) {
    printf ("Algorithm succeeded after %i steps\n", N);
  } else {
    printf ("Algorithm failed after %i steps\n", N);
  }
  printf ("X = %f\nY = %f\n", X, Y);
}

void Step (double * X1, double * Y1, double * X2, double * Y2, 
  double * X, double * Y) {
  *X = (*X1 + *X2) / 2.0;
  *Y = F (*X);
  if ((*Y1 < 0) && (*Y < 0)) {
    *X1 = *X;
    *Y1 = *Y;
  } else {
    *X2 = *X;
    *Y2 = *Y;
  }
}

double F (double X) {
  double Y;
  Y = 1.0 + 3.0 * X * (X + 4.0); 
  return Y;
}