Listing 2 Illustrates use of machine epsilon in a root-finding algorithm

/* root.c:
 *
 *   To use a different precision, change ftype
 *   and the suffix of the floating-point constants
 */

#include <stdio.h>
#include <float.h>
#include <math.h>
#include <assert.h>

#define sign(x) ((x < 0.0) ? -1 : (x > 0.0) ? 1 : 0)

#define PREC DBL_DIG
#define EPS DBL_EPSILON

typedef double ftype;

ftype root(ftype a, ftype b, ftype (*f)(ftype))
{

   ftype fofa = f(a);
   ftype fofb = f(b);
   assert(a < b);
   assert(fofa * fofb < 0.0);

   /* Close-in on root via bisection */
   while (fabs(b - a) > EPS*fabs(a))
   {
      ftype x = a + (b-a)/2.0;
      ftype fofx = f(x);
      if (x <= a || x >= b || fofx == 0.0)
         return x;
      if (sign(fofx) == sign(fofa))
      {
         a = x;
         fofa = fofx;
      }
      else
      {
         b = x;
         fofb = fofx;
      }
   }
   return a;
}

main()
{
   extern ftype f(ftype);
   printf("root == %.*f\n",PREC,root(-1.0,1.0,f));
   return 0;
}

ftype f(ftype x)
{
   return x*x + x - 1.0;
}

/* Output:
root == 0.618033988749895
*/
/* End of File */