When doing frequency synthesis with fractional-N PLLs, one often needs to find a rational approximation of a floating point number with the constraint that the numerator must not be larger than a certain number. The more exact the approximation is, the closer the actual frequency will be to the desired one.

The integer part is obviously easy, but the fractional part requires a more sophisticated algorithm. One such algorithm is based on Farey sequences and the formula for finding the next fraction between two Farey neighbors.

The Farey sequence of order N consists of all completely reduced fractions between 0 and 1. So e.g. the Farey sequence of order 3 consists of {0/1, 1/3, 1/2, 2/3, 1/1}. A sequence of a higher order contains all the terms of all lower orders and then some more. So, with a fractional-N PLL where the denominator is limited to some value D, the possible fractional parts of the N in the PLL is precisely the fractions present in the Farey sequency of order D. And our goal is to find the best one to approximate the fractional part of the desired N, i.e. a number between 0 and 1.

There is a neat and useful formula for finding the next fraction that will appear between two Farey neighbors as the order of the sequence is increased. If a/b < c/d are neighbors in some Farey sequence, the next term to appear between them is the *mediant* (a+c)/(b+d). So an algorithm to find better and better rational approximations to a number *x* is to

- Start with the Farey sequence of order 1, i.e. {0/1, 1/1}, where a/b = 0/1 and c/d = 1/1 are neighbors.
- Calculate the mediant (a+c)/(b+d) = (0+1)/(1+1) = 1/2.
- Figure out which of a/b and c/d are further from
*x*than the mediant and replace it with the mediant. - Go to 2 if the denominator is not yet larger than what is allowed.
- Use the best result previous to the denominator becoming too large.

This algorithm is described on this archived web page.

This seems fine, until one considers some special cases. Say e.g. the target number is 0.000 001 and the maximum allowed denominator is 2 000 000. The algorithm would narrow down the top end of the interval from 1/1 to 1/2 to 1/3 to 1/4… in each successive iteration. It would thus take a million steps before the perfect approximation of 1/1 000 000 is reached. This is hardly efficient in this case, even though it converges quickly in many (most) other cases.

A way to speed up these degenerate cases is to not just take a single step in each iteration, but instead figure out how many times in a row either a/b or c/d will be discarded in favor of the mediant. It turns out that this is not too hard to do.

If e.g. a/b is to be discarded K times in a row, the resulting number will be (a + K*c)/(b + K*d). So how many times K will a/b be discarded until it is time to narrow down the interval from the other end? Set (a + K*c)/(b + K*d) = x and solve for K, which gives K = (x*b – a)/(c – x*d). K is often not an integer, so then select the biggest integer smaller than K.

If it is instead c/d that is to be discarded, the formula for K is (c – x*d)/(x*b – a).

With this improvement, it takes just one step to get to 1/1 000 000 instead of a million. Quite an improvement.

I wrote some C-code to implement this improved algorithm. Here it is:

```
// Farey sequence-based rational approximation of numbers.
// Per Magnusson, 2024
// MIT licence, http://www.opensource.org/licenses/mit-license.php
#include <cstdint>
#include <cmath>
typedef struct {
uint32_t numerator;
uint32_t denominator;
} rational_t;
// Find the best rational approximation to a number between 0 and 1.
//
// target - a number between 0 and 1 (inclusive)
// maxdenom - the maximum allowed denominator
//
// The algorithm is based on Farey sequences/fractions. See
// https://web.archive.org/web/20181119092100/https://nrich.maths.org/6596
// a, b, c, d notation from
// https://en.wikipedia.org/wiki/Farey_sequence is used here (not
// from the above reference). I.e. narrow the interval between a/b
// and c/d by splitting it using the mediant (a+c)/(b+d) until we are
// close enough with either endpoint, or we have a denominator that is
// bigger than what is allowed.
// Start with the interval 0 to 1 (i.e. 0/1 to 1/1).
// A simple implementation of just calculating the mediant (a+c)/(b+d) and
// iterating with the mediant replacing the worst value of a/b and c/d is very
// inefficient in cases where the target is close to a rational number
// with a small denominator, like e.g. when approximating 10^-6.
// The straightforward algorithm would need about 10^6 iterations as it
// would try all of 1/1, 1/2, 1/3, 1/4, 1/5 etc. To resolve this slow
// convergence, at each step, it is calculated how many times the
// interval will need to be narrowed from the same side and all those
// steps are taken at once.
rational_t rational_approximation(double target, uint32_t maxdenom)
{
rational_t retval;
double mediant; // float does not have enough resolution
// to deal with single-digit differences
// between numbers above 10^8.
double N, Ndenom, Ndenom_min;
uint32_t a = 0, b = 1, c = 1, d = 1, ac, bd, Nint;
const int maxIter = 1000;
if(target > 1) {
// Invalid
retval.numerator = 1;
retval.denominator = 1;
return retval;
}
if(target < 0) {
// Invalid
retval.numerator = 0;
retval.denominator = 1;
return retval;
}
if(maxdenom < 1) {
maxdenom = 1;
}
mediant = 0;
Ndenom_min = 1/((double) 10*maxdenom);
int ii = 0;
// Farey approximation loop
while(1) {
ac = a+c;
bd = b+d;
if(bd > maxdenom || ii > maxIter) {
// The denominator has become too big, or too many iterations.
// Select the best of a/b and c/d.
if(target - a/(double)b < c/(double)d - target) {
ac = a;
bd = b;
} else {
ac = c;
bd = d;
}
break;
}
mediant = ac/(double)bd;
if(target < mediant) {
// Discard c/d as the mediant is closer to the target.
// How many times in a row should we do that?
// N = (c - target*d)/(target*b - a), but need to check for division by zero
Ndenom = target * (double)b - (double)a;
if(Ndenom < Ndenom_min) {
// Division by zero, or close to it!
// This means that a/b is a very good approximation
// as we would need to update the c/d side a
// very large number of times to get closer.
// Use a/b and exit the loop.
ac = a;
bd = b;
break;
}
N = (c - target * (double)d)/Ndenom;
Nint = floor(N);
// Check if the denominator will become too large
if(d + Nint*b > maxdenom) {
// Limit N, as the denominator would otherwise become too large
N = (maxdenom - d)/(double)b;
Nint = floor(N);
}
// Fast forward to a good c/d.
c = c + Nint*a;
d = d + Nint*b;
} else {
// Discard a/b as the mediant is closer to the target.
// How many times in a row should we do that?
// N = (target*b - a)/(c - target*d), but need to check for division by zero
Ndenom = (double)c - target * (double)d;
if(Ndenom < Ndenom_min) {
// Division by zero, or close to it!
// This means that c/d is a very good approximation
// as we would need to update the a/b side a
// very large number of times to get closer.
// Use c/d and exit the loop.
ac = c;
bd = d;
break;
}
N = (target * (double)b - a)/Ndenom;
Nint = floor(N);
// Check if the denominator will become too large
if(b + Nint*d > maxdenom) {
// Limit N, as the denominator would otherwise become too large
N = (maxdenom - b)/(double)d;
Nint = floor(N);
}
// Fast forward to a good a/b.
a = a + Nint*c;
b = b + Nint*d;
}
ii++;
}
retval.numerator = ac;
retval.denominator = bd;
return retval;
}
```

To test the algorithm on an Arduino (actually a Pi Pico), I wrote the following code:

```
#include <arduino.h>
typedef struct {
double target;
uint32_t maxdenom;
uint32_t expected_numerator;
uint32_t expected_denominator;
} rational_test_case_t;
void test_rational_approx()
{
rational_t result;
rational_test_case_t test[] = {
{0, 3000, 0, 1},
{1, 3000, 1, 1},
{0.5, 3000, 1, 2},
{0.5+1/3001.0, 3000, 751, 1501},
{1/3001.0, 2500, 1, 2500},
{1/3001.0, 1500, 0, 1},
{1/3001.0, 3001, 1, 3001},
{0.472757439, 1816, 564, 1193},
{0.472757439, 1817, 859, 1817},
};
uint32_t n_tests = sizeof(test)/sizeof(test[0]);
for(uint32_t ii = 0; ii < n_tests; ii++) {
result = rational_approximation(test[ii].target, test[ii].maxdenom);
Serial.print("target = ");
Serial.print(test[ii].target, 7);
Serial.print(", maxdenom = ");
Serial.print(test[ii].maxdenom);
Serial.print(", approx = ");
Serial.print(result.numerator);
Serial.print("/");
Serial.print(result.denominator);
if(result.numerator == test[ii].expected_numerator && result.denominator == test[ii].expected_denominator) {
Serial.println(" OK");
} else {
Serial.print(" Expected: ");
Serial.print(test[ii].expected_numerator);
Serial.print("/");
Serial.println(test[ii].expected_denominator);
}
}
}
```

The printout from the code is:

```
target = 0.0000000, maxdenom = 3000, approx = 0/1 OK
target = 1.0000000, maxdenom = 3000, approx = 1/1 OK
target = 0.5000000, maxdenom = 3000, approx = 1/2 OK
target = 0.5003332, maxdenom = 3000, approx = 751/1501 OK
target = 0.0003332, maxdenom = 2500, approx = 1/2500 OK
target = 0.0003332, maxdenom = 1500, approx = 0/1 OK
target = 0.0003332, maxdenom = 3001, approx = 1/3001 OK
target = 0.4727574, maxdenom = 1816, approx = 564/1193 OK
target = 0.4727574, maxdenom = 1817, approx = 859/1817 OK
```

I also implemented the algorithm in an Excel sheet:

The excel document is available here: