(EDIT: In response to grumpy comments, No it isn't homework. I am working on pitch detection, taking an array of potential harmonic peaks, and attempting to construct candidates for fundamental frequency. So, it is actually a very practical question.)
Consider the best fractional approximations for (eg) pi, ordered by increasing denominator: 3/1, 22/7, 355/113, ...
Challenge: create a tidy C algorithm that will generate the n'th quotient approximation a/b for a given float, returning also the discrepancy.
calcBestFrac(float frac, int n, int * a, int * b, float * err) {...}
The best technique I believe is continued fractions
Take away the fractional part of pi, and you get 3
Now, the remainder is 0.14159... = 1/7.06251..
So the next best rational is 3 + 1/7 = 22/7
Take away the 7 from 7.06251 and you get 0.06251.. Roughly 1/15.99659..
Call it 16, then the next best approximation is
3 + 1/(7 + 1/16) = 355/113
However, this is far from trivial to convert into clean C code. I will post if I get something tidy. In the meanwhile, someone may enjoy it as a brainteaser.
[Since you asked for this as an answer rather than a comment.]
For any real number, the convergents p[k]/q[k] of its continued fraction are always best rational approximations, but they aren't all the best rational approximations. To get all of them, you also have to take the semi-convergents/mediants — fractions of the form (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) for some integer n≥1. Taking n=a[k+2] gives p[k+2]/q[k+2], and the integers n to take are those from either floor(a[k+2]/2) or ceiling(a[k+2]/2), to a[k+2]. This is also mentioned on Wikipedia.
The continued fraction for π is [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…] (sequence A001203 in OEIS), the sequence of convergents is 3/1, 22/7, 333/106, 355/113, 103993/33102… (A002485/A002486), and the sequence of best approximations is 3/1, 13/4, 16/5, 19/6, 22/7, 179/57… (A063674/A063673).
So the algorithm says that the best approximations of π = [3; 7, 15, 1, 292, 1, 1,…] are
3/1 = [3]
13/4 = [3; 4]
16/5 = [3; 5]
19/6 = [3; 6]
22/7 = [3; 7]
179/57 = [3; 7, 8]
201/64 = [3; 7, 9]
223/71 = [3; 7, 10]
245/78 = [3; 7, 11]
267/85 = [3; 7, 12]
289/92 = [3; 7, 13]
311/99 = [3; 7, 14]
333/106 = [3; 7, 15]
355/113 = [3; 7, 15, 1]
52163/16604 = [3; 7, 15, 1, 146]
52518/16717 = [3; 7, 15, 1, 147]
… (all the fractions from [3; 7, 15, 1, 148] to [3; 7, 15, 1, 291])…
103993/33102 = [3; 7, 15, 1, 292]
104348/33215 = [3; 7, 15, 1, 292, 1]
...
Here's a C program that given a positive real number, generates its continued fraction, its convergents, and the sequence of best rational approximations. The function find_cf finds the continued fraction (putting the terms in a[] and the convergents in p[] and q[] — excuse the global variables), and the function all_best prints all the best rational approximations.
#include <math.h>
#include <stdio.h>
#include <assert.h>
// number of terms in continued fraction.
// 15 is the max without precision errors for M_PI
#define MAX 15
#define eps 1e-9
long p[MAX], q[MAX], a[MAX], len;
void find_cf(double x) {
int i;
//The first two convergents are 0/1 and 1/0
p[0] = 0; q[0] = 1;
p[1] = 1; q[1] = 0;
//The rest of the convergents (and continued fraction)
for(i=2; i<MAX; ++i) {
a[i] = lrint(floor(x));
p[i] = a[i]*p[i-1] + p[i-2];
q[i] = a[i]*q[i-1] + q[i-2];
printf("%ld: %ld/%ld\n", a[i], p[i], q[i]);
len = i;
if(fabs(x-a[i])<eps) return;
x = 1.0/(x - a[i]);
}
}
void all_best(double x) {
find_cf(x); printf("\n");
int i, n; long cp, cq;
for(i=2; i<len; ++i) {
//Test n = a[i+1]/2. Enough to test only when a[i+1] is even, actually...
n = a[i+1]/2; cp = n*p[i]+p[i-1]; cq = n*q[i]+q[i-1];
if(fabs(x-(double)cp/cq) < fabs(x-(double)p[i]/q[i]))
printf("%ld/%ld, ", cp, cq);
//And print all the rest, no need to test
for(n = (a[i+1]+2)/2; n<=a[i+1]; ++n) {
printf("%ld/%ld, ", n*p[i]+p[i-1], n*q[i]+q[i-1]);
}
}
}
int main(int argc, char **argv) {
double x;
if(argc==1) { x = M_PI; } else { sscanf(argv[1], "%lf", &x); }
assert(x>0); printf("%.15lf\n\n", x);
all_best(x); printf("\n");
return 0;
}
For π, here's the output of this program, in about 0.003 seconds (i.e., it's truly better than looping through all possible denominators!), line-wrapped for readability:
% ./a.out
3.141592653589793
3: 3/1
7: 22/7
15: 333/106
1: 355/113
292: 103993/33102
1: 104348/33215
1: 208341/66317
1: 312689/99532
2: 833719/265381
1: 1146408/364913
3: 4272943/1360120
1: 5419351/1725033
14: 80143857/25510582
13/4, 16/5, 19/6, 22/7, 179/57, 201/64, 223/71, 245/78, 267/85, 289/92, 311/99,
333/106, 355/113, 52163/16604, 52518/16717, 52873/16830, 53228/16943, 53583/17056,
53938/17169, 54293/17282, 54648/17395, 55003/17508, 55358/17621, 55713/17734,
56068/17847, 56423/17960, 56778/18073, 57133/18186, 57488/18299, 57843/18412,
58198/18525, 58553/18638, 58908/18751, 59263/18864, 59618/18977, 59973/19090,
60328/19203, 60683/19316, 61038/19429, 61393/19542, 61748/19655, 62103/19768,
62458/19881, 62813/19994, 63168/20107, 63523/20220, 63878/20333, 64233/20446,
64588/20559, 64943/20672, 65298/20785, 65653/20898, 66008/21011, 66363/21124,
66718/21237, 67073/21350, 67428/21463, 67783/21576, 68138/21689, 68493/21802,
68848/21915, 69203/22028, 69558/22141, 69913/22254, 70268/22367, 70623/22480,
70978/22593, 71333/22706, 71688/22819, 72043/22932, 72398/23045, 72753/23158,
73108/23271, 73463/23384, 73818/23497, 74173/23610, 74528/23723, 74883/23836,
75238/23949, 75593/24062, 75948/24175, 76303/24288, 76658/24401, 77013/24514,
77368/24627, 77723/24740, 78078/24853, 78433/24966, 78788/25079, 79143/25192,
79498/25305, 79853/25418, 80208/25531, 80563/25644, 80918/25757, 81273/25870,
81628/25983, 81983/26096, 82338/26209, 82693/26322, 83048/26435, 83403/26548,
83758/26661, 84113/26774, 84468/26887, 84823/27000, 85178/27113, 85533/27226,
85888/27339, 86243/27452, 86598/27565, 86953/27678, 87308/27791, 87663/27904,
88018/28017, 88373/28130, 88728/28243, 89083/28356, 89438/28469, 89793/28582,
90148/28695, 90503/28808, 90858/28921, 91213/29034, 91568/29147, 91923/29260,
92278/29373, 92633/29486, 92988/29599, 93343/29712, 93698/29825, 94053/29938,
94408/30051, 94763/30164, 95118/30277, 95473/30390, 95828/30503, 96183/30616,
96538/30729, 96893/30842, 97248/30955, 97603/31068, 97958/31181, 98313/31294,
98668/31407, 99023/31520, 99378/31633, 99733/31746, 100088/31859, 100443/31972,
100798/32085, 101153/32198, 101508/32311, 101863/32424, 102218/32537, 102573/32650,
102928/32763, 103283/32876, 103638/32989, 103993/33102, 104348/33215, 208341/66317,
312689/99532, 833719/265381, 1146408/364913, 3126535/995207,
4272943/1360120, 5419351/1725033, 42208400/13435351, 47627751/15160384,
53047102/16885417, 58466453/18610450, 63885804/20335483, 69305155/22060516,
74724506/23785549, 80143857/25510582,
All these terms are correct, though if you increase MAX you start getting errors because of precision. I'm myself impressed with how many terms you get with only 13 convergents. (As you can see, there's a small bug where it sometimes doesn't print the very first "n/1" approximation, or prints it incorrectly — I leave it for you to fix!)
You can try with √2, whose continued fraction is [1; 2, 2, 2, 2…]:
% ./a.out 1.41421356237309504880
1.414213562373095
1: 1/1
2: 3/2
2: 7/5
2: 17/12
2: 41/29
2: 99/70
2: 239/169
2: 577/408
2: 1393/985
2: 3363/2378
2: 8119/5741
2: 19601/13860
2: 47321/33461
3/2, 4/3, 7/5, 17/12, 24/17, 41/29, 99/70, 140/99, 239/169, 577/408, 816/577, 1393/985, 3363/2378, 4756/3363, 8119/5741, 19601/13860, 47321/33461,
Or the golden ratio φ = (1+√5)/2 whose continued fraction is [1; 1, 1, 1, …]:
% ./a.out 1.61803398874989484820
1.618033988749895
1: 1/1
1: 2/1
1: 3/2
1: 5/3
1: 8/5
1: 13/8
1: 21/13
1: 34/21
1: 55/34
1: 89/55
1: 144/89
1: 233/144
1: 377/233
2/1, 3/2, 5/3, 8/5, 13/8, 21/13, 34/21, 55/34, 89/55, 144/89, 233/144, 377/233,
(See the Fibonacci numbers? Here the convergents are all the approximants.)
Or with rational numbers like 4/3 = [1; 3]:
% ./a.out 1.33333333333333333333
1.333333333333333
1: 1/1
3: 4/3
3/2, 4/3,
or 14/11 = [1; 3, 1, 2]:
% ./a.out 1.27272727272727272727
1.272727272727273
1: 1/1
3: 4/3
1: 5/4
2: 14/11
3/2, 4/3, 5/4, 9/7, 14/11,
Enjoy!