/* The following values are used to test for too small
   and too large values of v for an integer fraction to
   represent.

   These values are currently set for the amdal 5890.
*/

#define MAX 4.0e+10
#define MIN 3.0e-10

double
frac(v, n, d, error)
double v, error;
int *n, *d;
{
	/*
	  given a number, v, this function outputs two integers,
	  d and n, such that

		v = n / d

	  to accuracy
		epsilon = | (v - n/d) / v | <= error

	  input:  v = decimal number you want replaced by fraction.
		  error = accuracy to which the fraction should
			  represent the number v.

	  output:  n = numerator of representing fraction.
		   d = denominator of representing fraction.

	  return value:  -1.0 if (v < MIN || v > MAX || error < 0.0)
			 | (v - n/d) / v | otherwise.

	  Note:  This program only works for positive numbers, v.

	  reference:  Jerome Spanier and Keith B. Oldham, "An Atlas
		of Functions," Springer-Verlag, 1987, pp. 665-7.
	*/

	int D, N, t;
	double epsilon, r, m, fabs();


	if (v < MIN || v > MAX || error < 0.0)
		return(-1.0);
	*d = D = 1;
	*n = (int)v;
	N = (*n) + 1;
	goto three;

one:	if (r > 1.0)
		goto two;
	r = 1.0/r;
two:	N += (*n)*(int)r;
	D += (*d)*(int)r;
	(*n) += N;
	(*d) += D;
three:	r = 0.0;
	if (v*(*d) == (double)(*n))
		goto four;
	r = (N - v*D)/(v*(*d) - (*n));
	if (r > 1.0)
		goto four;
	t = N;
	N = (*n);
	*n = t;
	t = D;
	D = (*d);
	*d = t;
four:	printf("%d/%d", *n, *d);
	epsilon = fabs(1.0 - (*n)/(v*(*d)));
	if (epsilon <= error)
		goto six;
	m = 1.0;
	do {
		m *= 10.0;
	} while (m*epsilon < 1.0);
	epsilon = 1.0/m * ((int)(0.5 + m*epsilon));
six:	printf("	epsilon = %e\n", epsilon);
	if (epsilon <= error)
		return(epsilon);
	if (r != 0.0)
		goto one;
}

--------------

The following program calls frac.c and determines several rational
approximations for pi.

--------------

main()
{
	double frac(), error = 1.0e-10;
	int n, d;

	frac(3.1415926535, &n, &d, error);
}

--------------

Typical output would be:

hrdcpy cc testfrac.c frac.c -lm
+ /etc/preroot /native /bin/cc testfrac.c frac.c -lm 
testfrac.c:
frac.c:
hrdcpy a.out
3/1	epsilon = 5.000000e-02
22/7	epsilon = 4.000000e-04
355/113	epsilon = 8.000000e-08
104348/33215	epsilon = 1.000000e-10
208341/66317	epsilon = 1.036518e-11

If you have any other questions, just send me mail.


			Robert J. Craig
			AT&T Bell Laboratories
			1200 East Naperville Road
			Naperville, IL 60566-7045