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