home *** CD-ROM | disk | FTP | other *** search
/ Point Programming 1 / PPROG1.ISO / c / snippets / fraction.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-04-03  |  2.8 KB  |  92 lines
  1. /*
  2. **  FRACTION.C - Compute continued fraction series
  3. **
  4. **  cfrac() donated to the public domain by the author, Thad Smith
  5. **  original Fraction.C, public domain by Bob Stout, modified to use cfrac()
  6. */
  7.  
  8. #include <math.h>
  9. #include <float.h>
  10.  
  11. #define MAX_LENGTH 100
  12.  
  13. long double cfrac(long double x, long double *p, long double *q, int bits)
  14. {
  15.       double v;                                 /* integer in series    */
  16.       long double del;                          /* x - v                */
  17.       long double z;    /* approximated value from truncated series     */
  18.       long double t;                            /* temp                 */
  19.       long double p0 = 0.0, q0 = 0.0;           /* last p, q            */
  20.       long double imax;                         /* max for p, q         */
  21.       static long double cf[MAX_LENGTH]; /* continued fraction integers */
  22.       int i, j, ntimes = MAX_LENGTH;;
  23.  
  24.       if (x < 0)
  25.             x = -x;
  26.       imax = floor(pow(2.0, bits)) - 1.0;
  27.       for (i = 0; i < ntimes; i++)
  28.       {
  29.             v = floor((double)x);
  30.             cf[i] = v;
  31.             z = cf[i];
  32.             *p = z; *q = 1;
  33.             for (j = i; j--; )
  34.             {
  35.                   z = cf[j] + 1.0/z;
  36.                   t = *p;
  37.                   *p = cf[j] * (*p) + (*q);
  38.                   *q = t;
  39.             }
  40.             del = x-v;
  41.             if (del < DBL_EPSILON)
  42.                   break;
  43.             if ((*p > imax) || (*q > imax))
  44.             {
  45.                   *p = p0;
  46.                   *q = q0;
  47.                   break;
  48.             }
  49.             else
  50.             {
  51.                   p0 = *p;
  52.                   q0 = *q;
  53.             }
  54.             x = 1.0 / del;
  55.       }
  56.       return (*p)/(*q);
  57. }
  58.  
  59. /*
  60. **  Remove everything below this to use cfrac() as a stand-alone function
  61. */
  62.  
  63. #include <stdio.h>
  64. #include <stdlib.h>
  65.  
  66. main (int argc, char *argv[])
  67. {
  68.       long double x;      /* value to be approximated */
  69.       long double r,p,q;  /* approx ratio r = p/q     */
  70.       int bits;           /* bits of precision        */
  71.  
  72.       if (argc < 2 || argc > 3)
  73.       {
  74.             puts ("Use: FRACTION value [precision]");
  75.             puts ("where value = floating point value to generate "
  76.                   "continued fraction");
  77.             puts ("      precision (optional) = bits in "
  78.                   "numerator/denominator");
  79.             return 1;
  80.       }
  81.       sscanf (argv[1], "%Lf", &x);
  82.       if (argc == 3)
  83.             bits = atoi(argv[2]);
  84.       else  bits = 32;
  85.  
  86.       cfrac(x, &p, &q, bits);
  87.       printf("\n[%.20Lf]\n%.0Lf/%.0Lf = %lXh/%lXh = %.20Lf\n",
  88.             x, p, q, (long)p, (long)q, r = p/q);
  89.       printf("Error = %.10Lg, (%.10Lf%%)\n", r - x, 100. * (r - x) / x);
  90.       return EXIT_SUCCESS;
  91. }
  92.