home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Media Share 13
/
mediashare_13.zip
/
mediashare_13
/
ZIPPED
/
PROGRAM
/
SNIP9404.ZIP
/
FRACTION.C
< prev
next >
Wrap
C/C++ Source or Header
|
1994-04-03
|
3KB
|
92 lines
/*
** 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;
}