home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_13_03 / allison / root.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-03  |  1.2 KB  |  65 lines
  1. LISTING 2 - Illustrates use of Machine Epsilon in a root-finding
  2. algorithm
  3.  
  4. /* root.c:
  5.  *
  6.  *   To use a different precision, change ftype
  7.  *   and the suffix of the floating-point constants
  8.  */
  9.  
  10. #include <stdio.h>
  11. #include <float.h>
  12. #include <math.h>
  13. #include <assert.h>
  14.  
  15. #define sign(x) ((x < 0.0) ? -1 : (x > 0.0) ? 1 : 0)
  16.  
  17. #define PREC DBL_DIG
  18. #define EPS  DBL_EPSILON
  19.  
  20. typedef double ftype;
  21.  
  22. ftype root(ftype a, ftype b, ftype (*f)(ftype))
  23. {
  24.     ftype fofa = f(a);
  25.     ftype fofb = f(b);
  26.     assert(a < b);
  27.     assert(fofa * fofb < 0.0);
  28.  
  29.     /* Close-in on root via bisection */
  30.     while (fabs(b - a) > EPS*fabs(a))
  31.     {
  32.         ftype x = a + (b-a)/2.0;
  33.         ftype fofx = f(x);
  34.         if (x <= a || x >= b || fofx == 0.0)
  35.             return x;
  36.         if (sign(fofx) == sign(fofa))
  37.         {
  38.             a = x;
  39.             fofa = fofx;
  40.         }
  41.         else
  42.         {
  43.             b = x;
  44.             fofb = fofx;
  45.         }
  46.     }
  47.     return a;
  48. }
  49.  
  50. main()
  51. {
  52.     extern ftype f(ftype);
  53.     printf("root == %.*f\n",PREC,root(-1.0,1.0,f));
  54.     return 0;
  55. }
  56.  
  57. ftype f(ftype x)
  58. {
  59.     return x*x + x - 1.0;
  60. }
  61.  
  62. /* Output:
  63. root == 0.618033988749895
  64. */
  65.