Louis is a member of the Department of Industrial and Manufacturing Systems Engineering at Lehigh University. He can be contacted at [email protected].
In integer arithmetic, a common practice for calculating a fractional part of a value is to multiply by the numerator of some fraction and divide by the denominator. This is such a common operation for Forth programmers (who frequently work on small, integer-only machines) that one of the most basic Forth words is */, which does a multiplication followed by a division with a double-precision intermediate result. But even in the current age of sophisticated, single-chip computers with floating-point instructions, the need to apply the old chestnut can arise. Thus, to calculate 0.6 of some quantity, you would multiply by 3 and divide by 5. In this case, the choice of a common fraction is trivial, but it isn't always: Accuracy and magnitude may be of concern.
Suppose you need to multiply by the fractional part of 
3=1.7320508+. A simplistic approach would be to multiply by 732 and divide by 1000 for an accuracy of 5E-5, but this might not be the best solution. Unless a double-precision intermediate result is provided, overflow can occur. For example, if a multiplication routine with an unsigned 16-bit product was preferred, overflow would occur for any number greater than 216/732=89. Slightly better accuracy with less overflow complication can be achieved by multiplying by 112/153=0.732026+, which has an approximation error of 2.5E-5. To restrict the numerator to two or even one decimal digit, the best common fractions would be 71/97=0.7319587+ for an error of 9.2E5 or 8/11=0.72727+ for an error of 4.8E3. If the denominator were restricted to an 8-bit binary number so that an 8-bit divisor routine could be used, then the best common fraction would be 153/209=0.7320574+ for an error of 6.61E-6.
Even if multiple-precision multiplication and division are available, it is often necessary to obtain a better approximation. For instance, to multiply by the decimal part of 
=3.141592653+ with an accuracy of 1E-6, the simple approach would be to multiply by 141,593 and divide by 100,000, which would give an accuracy of 3.54E-7. However, the fraction 16/113 approximates the decimal part of to an error of 2.67E-7: the same order of magnitude of error without the need for multiple-precision arithmetic. If you don't mind using 16-bit divisors, 4703/33215=0.1415926539+ is accurate to within 3.32E-10. Most programmers prefer the latter two approximations.
The question, What is the best common fraction to use for a particular application? is most likely to take one of three slightly different forms:
- What is the common fraction with the smallest numerator and denominator that will approximate a given real number to within a specified tolerance?
- What is the common fraction with a numerator less than N that best approximates a given real number?
- What is the common fraction with a denominator less than N that best approximates a given real number?
Farey was a mineralogist who, in 1816, made observations about common fractions and published them without proof. The mathematician Cauchy published the proofs for Farey's observations shortly thereafter. Farey wrote down the sequence of all reduced fractions between 0 and 1 whose denominators are limited by a number N, the "order" of the Farey Sequence. Figure 1 shows the Farey Sequences of orders 1, 2, 3, 4, and 5.
Farey observed (and Cauchy proved) that for two consecutive fractions a/b<c/d of the Farey Sequence of order N, ad-bc=-1. This theorem is proved by induction. By inspection, it can be seen that it is true for orders 1, 2, 3, 4, and 5. Assuming that it is true for a series of order N, it can be proved that it is true for a series of order N+1, thus completing the inductive proof.
The "mediants" of a Farey Sequence are also interesting. The fraction (a+c)/(b+d) is called the mediant of two consecutive terms of a sequence a/b and c/d. For instance, a/b<(a+c)/(b+d)<c/d since a(b+d)<b(a+c)
ad<bca/b<c/d and (a+c)d<(b+d)cad<bca/b<c/d. Assume that a/b and c/d are consecutive terms of a Farey Sequence of order N. Now consider a Farey Sequence of order N+1 with a subsequence {a/b,h/k,c/d} contained in the sequence of order N+1 but not of order N. If this is true, then ak-bh=-1, hd-kc=-1. Solving these equations for h and k yields h=(a+c)/(bc-ad) and k=(b+d)/(bc-ad). But since a/b and c/d are consecutive terms of a Farey Sequence of order N, ab-bc=1, which reduces the expressions for h and k to h=a+c and k=c+d. This illustrates a second theorem associated with Farey Sequences: Fractions that belong to the Farey Sequence of order N+1 but not of order N are mediants of the Farey Sequence of order N.
This important theorem means that starting from the Farey Sequence [0/1,1/1], sequences of a higher order can be constructed by successively inserting mediants with appropriate denominators. This enables efficient search procedures to be defined for answering the three questions posed earlier. I'll describe the answer to Question 1; the procedures for the other two questions are simple modifications to the first.
Question 1 seeks a common fraction that will approximate a real number 
to within a given tolerance 
. We need an efficient way to search all possible Farey Sequences. Suppose you know that is contained in some interval [a/b,c/d] of a Farey Sequence. If neither endpoint approximates to within the desired tolerance , another approximation to can be made by calculating the mediant of this Farey interval. This will divide the interval into the two intervals [a/b,h/k] and [h/k,c/d], where h/k=(a+c)/(c+d). If |-h/k|
, h/k yields the desired common fraction, and you are finished. Otherwise, you need only determine which interval contains , replace the appropriate interval limit with h/k, and repeat the procedure.
This search procedure examines those intervals containing of successively higher-order Farey Sequences until an interval limit is found that approximates to the desired tolerance. Assuming that neither 0 or 1 are a satisfactory approximation to , the steps of the search procedure are:
- Initialize endpoints of the first Farey Sequence [a/b,c/d] to a=0,b=c,d=1.
- Calculate an approximation h/k, where h=a+c,k=b+d.
- Calculate error=(h/k)-. If |error|, stop. h/k is the desired common fraction.
- If error<0, change the lower limits to a=h,b=k. Go to step 2.
- If error>0, change the upper limits to c=h,d=k. Go to step 2.
and . When programmed on a programmable calculator, the execution time ranges from a second to a minute or more.
Answering the maximum-numerator question is nearly as straightforward. Intervals on successively higher Farey Sequences are searched until either the numerator of the new interval boundary in the next-higher sequence exceeds the limit N or the desired tolerance is achieved. If the absolute best approximation is desired with the magnitude constraint, a tolerance of 0 can be specified. The major modification is that if the search loop is exited due to the numerator magnitude constraint, both endpoints of the interval [a/b,c/d] must be checked for the best approximation. Listing Two is a C function that implements this procedure.
The procedure for the maximum-denominator problem is similar to that of the maximum numerator. Intervals on successively higher Farey Sequences are searched until either the denominator of the new interval boundary in the next-higher sequence exceeds the limit max_denominator or the desired tolerance is achieved; see Listing Three.
Reference
Rademacher, Hans. Lectures on Elementary Number Theory. New York, NY: Blaisdell Publishing Company, 1964.
Figure 1: Farey Sequences of orders 1, 2, 3, 4, and 5
Copyright © 1995, Dr. Dobb's Journal
0 1
- -
1, 1
0 1 1
- - -
1, 2, 1
0 1 1 2 1
- - - - -
1, 3, 2, 3, 1
0 1 1 1 2 3 1
- - - - - - -
1, 4, 3, 2, 3, 4, 1
0 1 1 1 2 1 3 2 3 4 1
- - - - - - - - - - -
1, 5, 4, 3, 5, 2, 5, 3, 4, 5, 1
Listing One
#include <stdio.h>
#include <math.h>
void Fraction( double realnum, double maxerror, double* num, double* denom) {
double a,b,c,d,h,k,error;
a = 0;
b = c = d = 1;
for( ;; ) {
h = a + c;
k = b + d;
error = h/k - realnum;
if ( fabs(error) <= maxerror )
break;
if ( error < 0 ) {
a = h;
b = k;
}else{
c = h;
d = k;
}
}
*num = h;
*denom = k;
}
void main() {
double realnum, maxerror, num, denom;
printf( "\nEnter real number to be approximated > ");
scanf("%lf", &realnum);
printf( "\nEnter maximum allowable error > ");
scanf("%lf",&maxerror);
Fraction( realnum, maxerror, &num, &denom);
printf( "\n real number = %e", realnum);
printf( "\n numerator = %.0f", num);
printf( "\n denominator = %.0f", denom);
printf( "\nnumerator / denominator = %e", num/denom);
printf( "\n error in approx = %e", num/denom - realnum);
}
Listing Two
#include <stdio.h>
#include <math.h>
void FractionLimitNumerator( double realnum, double maxerror,
double* num, double* denom, double maxnum) {
double a,b,c,d,h,k,error;
a = 0;
b = c = d = 1;
for( ;; ) {
h = a + c;
k = b + d;
if ( h > maxnum ) {
if ( fabs(a/b - realnum) < fabs(c/d - realnum) ) {
h = a;
k = b;
}else{
h = c;
k = d;
}
break;
}
error = h/k - realnum;
if ( fabs(error) <= maxerror )
break;
if ( error < 0 ) {
a = h;
b = k;
}else{
c = h;
d = k;
}
}
*num = h;
*denom = k;
}
void main() {
double realnum, maxerror, num, denom, maxnum;
printf( "\nEnter real number to be approximated > ");
scanf("%lf", &realnum);
printf( "\nEnter maximum allowable error > ");
scanf("%lf",&maxerror);
printf( "\nEnter maximum allowable numerator > ");
scanf("%lf",&maxnum);
FractionLimitNumerator( realnum, maxerror, &num, &denom, maxnum);
printf( "\n real number = %e", realnum);
printf( "\n numerator = %.0f", num);
printf( "\n denominator = %.0f", denom);
printf( "\nnumerator / denominator = %e", num/denom);
printf( "\n error in approx = %e", num/denom - realnum);
}
Listing Three
#include <stdio.h>
#include <math.h>
void FractionLimitDenominator( double realnum, double maxerror,
double* num, double* denom, double maxdenom) {
double a,b,c,d,h,k,error;
a = 0;
b = c = d = 1;
for( ;; ) {
h = a + c;
k = b + d;
if ( k > maxdenom ) {
if ( fabs(a/b - realnum) < fabs(c/d - realnum) ) {
h = a;
k = b;
}else{
h = c;
k = d;
}
break;
}
error = h/k - realnum;
if ( fabs(error) <= maxerror )
break;
if ( error < 0 ) {
a = h;
b = k;
}else{
c = h;
d = k;
}
}
*num = h;
*denom = k;
}
void main() {
double realnum, maxerror, num, denom, maxdenom;
printf( "\nEnter real number to be approximated > ");
scanf("%lf", &realnum);
printf( "\nEnter maximum allowable error > ");
scanf("%lf",&maxerror);
printf( "\nEnter maximum allowable denominator > ");
scanf("%lf",&maxdenom);
FractionLimitDenominator( realnum, maxerror, &num, &denom, maxdenom);
printf( "\n real number = %e", realnum);
printf( "\n numerator = %.0f", num);
printf( "\n denominator = %.0f", denom);
printf( "\nnumerator / denominator = %e", num/denom);
printf( "\n error in approx = %e", num/denom - realnum);
}
