cf_fraction v. 1.0 algorithm
Method by Dr. Peterson, implementation by Michel Vuijlsteke
Here's the method used in <cf_fraction>, with a partial explanation for how it works:
We want to approximate a value m (given as a decimal) between 0 and 1,
by a fraction Y/X. Think of fractions as vectors (denominator,
numerator), so that the slope of the vector is the value of the
fraction. We are then looking for a lattice vector (X, Y) whose slope
is as close as possible to m. This picture illustrates the goal, and
shows that, given two vectors A and B on opposite sides of the desired
slope, their sum A + B = C is a new vector whose slope is between the
two, allowing us to narrow our search:
num
^

+ + + + + + + + + + +

+ + + + + + + + + + +
 slope m=0.7
+ + + + + + + + + + + /
 /
+ + + + + + + + + + D < solution
 /
+ + + + + + + + + /+ +
 /
+ + + + + + + C/ + + +
 /
+ + + + + + /+ + + + +
 /
+ + + + B/ + + + + + +
 /
+ + + /A + + + + + + +
 /
+ +/ + + + + + + + + +
 /
+++++++++++> denom
Here we start knowing the goal is between A = (3,2) and B = (4,3), and
formed a new vector C = A + B. We test the slope of C and find that
the desired slope m is between A and C, so we continue the search
between A and C. We add A and C to get a new vector D = A + 2*B, which
in this case is exactly right and gives us the answer.
Given the vectors A and B, with slope(A) < m < slope(B),
we can find consecutive integers M and N such that
slope(A + M*B) < x < slope(A + N*B) in this way:
If A = (b, a) and B = (d, c), with a/b < m < c/d, solve
a + x*c  = m 

b + x*d 
to give
x =  b*m  a 

c  d*m 
If this is an integer (or close enough to an integer to consider it
so), then A + x*B is our answer. Otherwise, we round it down and up to
get integer multipliers M and N respectively, from which new lower and
upper bounds A' = A + M*B and B' = A + N*B can be obtained. Repeat the
process until the slopes of the two vectors are close enough for the
desired accuracy. The process can be started with vectors (0,1), with
slope 0, and (1,1), with slope 1. Surprisingly, this process produces
exactly what continued fractions produce, and therefore it will
terminate at the desired fraction (in lowest terms, as far as I can
tell) if there is one, or when it is correct within the accuracy of
the original data.
For example, for the slope 0.7 shown in the picture above, we get
these approximations:
Step 1: A = 0/1, B = 1/1 (a = 0, b = 1, c = 1, d = 1)
x = 
1 * 0.7  0 0.7 
= 
0.7 
= 2.3333 


1  1 * 0.7 
0.3 
M = 2: lower bound A' = (0 + 2*1) / (1 + 2*1) = 2 / 3
N = 3: upper bound B' = (0 + 3*1) / (1 + 3*1) = 3 / 4
Step 2: A = 2/3, B = 3/4 (a = 2, b = 3, c = 3, d = 4)
x = 
3 * 0.7  2 
= 
0.1 
= 0.5 


3  4 * 0.7 
0.2 
M = 0: lower bound A' = (2 + 0*3) / (3 + 0*4) = 2 / 3
N = 1: upper bound B' = (2 + 1*3) / (3 + 1*4) = 5 / 7
Step 3: A = 2/3, B = 5/7 (a = 2, b = 3, c = 5, d = 7)
x = 
3 * 0.7  2 
= 
0.1 
= 1 


5  7 * 0.7 
0.1 
N = 1: exact value A' = B' = (2 + 1*5) / (3 + 1*7) = 7 / 10
which of course is obviously right.
In most cases you will never get an exact integer, because of rounding
errors, but can stop when one of the two fractions is equal to the
goal to the given accuracy.
