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;
}
|
