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