summaryrefslogtreecommitdiff
path: root/reference/C/CONTRIB/SNIP/fraction.c
blob: 35a4b4c78b10db8faea34c0062bc5d77f103489d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
/*
**  FRACTION.C - Compute continued fraction series
**
**  cfrac() donated to the public domain by the author, Thad Smith
**  original Fraction.C, public domain by Bob Stout, modified to use cfrac()
*/

#include <math.h>
#include <float.h>

#define MAX_LENGTH 100

long double cfrac(long double x, long double *p, long double *q, int bits)
{
      double v;                                 /* integer in series    */
      long double del;                          /* x - v                */
      long double z;    /* approximated value from truncated series     */
      long double t;                            /* temp                 */
      long double p0 = 0.0, q0 = 0.0;           /* last p, q            */
      long double imax;                         /* max for p, q         */
      static long double cf[MAX_LENGTH]; /* continued fraction integers */
      int i, j, ntimes = MAX_LENGTH;;

      if (x < 0)
            x = -x;
      imax = floor(pow(2.0, bits)) - 1.0;
      for (i = 0; i < ntimes; i++)
      {
            v = floor((double)x);
            cf[i] = v;
            z = cf[i];
            *p = z; *q = 1;
            for (j = i; j--; )
            {
                  z = cf[j] + 1.0/z;
                  t = *p;
                  *p = cf[j] * (*p) + (*q);
                  *q = t;
            }
            del = x-v;
            if (del < DBL_EPSILON)
                  break;
            if ((*p > imax) || (*q > imax))
            {
                  *p = p0;
                  *q = q0;
                  break;
            }
            else
            {
                  p0 = *p;
                  q0 = *q;
            }
            x = 1.0 / del;
      }
      return (*p)/(*q);
}

/*
**  Remove everything below this to use cfrac() as a stand-alone function
*/

#include <stdio.h>
#include <stdlib.h>

main (int argc, char *argv[])
{
      long double x;      /* value to be approximated */
      long double r,p,q;  /* approx ratio r = p/q     */
      int bits;           /* bits of precision        */

      if (argc < 2 || argc > 3)
      {
            puts ("Use: FRACTION value [precision]");
            puts ("where value = floating point value to generate "
                  "continued fraction");
            puts ("      precision (optional) = bits in "
                  "numerator/denominator");
            return 1;
      }
      sscanf (argv[1], "%Lf", &x);
      if (argc == 3)
            bits = atoi(argv[2]);
      else  bits = 32;

      cfrac(x, &p, &q, bits);
      printf("\n[%.20Lf]\n%.0Lf/%.0Lf = %lXh/%lXh = %.20Lf\n",
            x, p, q, (long)p, (long)q, r = p/q);
      printf("Error = %.10Lg, (%.10Lf%%)\n", r - x, 100. * (r - x) / x);
      return EXIT_SUCCESS;
}