home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1994 #1 / monster.zip / monster / PROG_C / SNPD9404.ZIP / FRACTION.C < prev    next >
C/C++ Source or Header  |  1994-04-03  |  3KB  |  91 lines
  1. .I 0 5
  2. /*
  3. **  FRACTION.C - Compute continued fraction series
  4. **
  5. **  cfrac() donated to the public domain by the author, Thad Smith
  6. **  original Fraction.C, public domain by Bob Stout, modified to use cfrac()
  7. .D 1 4
  8. .I 6 55
  9. #include <math.h>
  10. #include <float.h>
  11.  
  12. #define MAX_LENGTH 100
  13.  
  14. long double cfrac(long double x, long double *p, long double *q, int bits)
  15. {
  16.       double v;                                 /* integer in series    */
  17.       long double del;                          /* x - v                */
  18.       long double z;    /* approximated value from truncated series     */
  19.       long double t;                            /* temp                 */
  20.       long double p0 = 0.0, q0 = 0.0;           /* last p, q            */
  21.       long double imax;                         /* max for p, q         */
  22.       static long double cf[MAX_LENGTH]; /* continued fraction integers */
  23.       int i, j, ntimes = MAX_LENGTH;;
  24.  
  25.       if (x < 0)
  26.             x = -x;
  27.       imax = floor(pow(2.0, bits)) - 1.0;
  28.       for (i = 0; i < ntimes; i++)
  29.       {
  30.             v = floor((double)x);
  31.             cf[i] = v;
  32.             z = cf[i];
  33.             *p = z; *q = 1;
  34.             for (j = i; j--; )
  35.             {
  36.                   z = cf[j] + 1.0/z;
  37.                   t = *p;
  38.                   *p = cf[j] * (*p) + (*q);
  39.                   *q = t;
  40.             }
  41.             del = x-v;
  42.             if (del < DBL_EPSILON)
  43.                   break;
  44.             if ((*p > imax) || (*q > imax))
  45.             {
  46.                   *p = p0;
  47.                   *q = q0;
  48.                   break;
  49.             }
  50.             else
  51.             {
  52.                   p0 = *p;
  53.                   q0 = *q;
  54.             }
  55.             x = 1.0 / del;
  56.       }
  57.       return (*p)/(*q);
  58. }
  59.  
  60. /*
  61. **  Remove everything below this to use cfrac() as a stand-alone function
  62. */
  63.  
  64. .I 8 25
  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. .D 9 70
  91.