#include <u.h>
#include <libc.h>

/*
 *		Robert J. Craig
 * 		AT&T Bell Laboratories
 *		1200 East Naperville Road
 *		Naperville, IL 60566-7045
 *
 * 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.
 */

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

double
frac(double v, int *n, int *d, double error)
{

	int D, N, t;
	double epsilon, r, m;


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

	do{
		r = 0.0;
		if(v*(*d) != (double)*n){
			r = (N - v*D)/(v * *d - *n);
			if(r <= 1.0){
				t = N;
				N = *n;
				*n = t;
				t = D;
				D = *d;
				*d = t;
			}
		}
		print("%6d/%-6d", *n, *d);
		epsilon = fabs(1.0 - *n/(v * *d));
		if(epsilon > error){
			m = 1.0;
			do{
				m *= 10.0;
			}while (m*epsilon < 1.0);
			epsilon = 1.0/m * (int)(0.5 + m*epsilon);
		}
	
		print("	ε = %e\n", epsilon);
		if(epsilon <= error)
			break;

		if(r <= 1.0)
			r = 1.0/r;

		N += *n *(int)r;
		D += *d *(int)r;
		(*n) += N;
		(*d) += D;

	}while(r != 0.0);
	return epsilon;
}

void
main(int argc, char *argv[])
{
	int n, d;
	double error = 1.0e-10;

	if(argc != 2){
		fprint(2, "usage: fract: n\n");
		exits("usage");
	}

	frac(atof(argv[1]), &n, &d, error);
}