Algorithm Challenge: Generate Continued Fractions for a float

Question

(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.

Warily

If this is truly a brain teaser, it's off topic. If it's actually your homework, it could be on-topic, but you have to show us what you've tried so far and what specific problem you're having. We aren't just going to do your work for you. – Wiggs

@Ohmu, is this homework related? – Ligneous

I looks like you want the next fraction with better accuracy, but you don't state that. If you don't actually want that restriction, then shouldn't it be more like 3/1, 6/2, 9/3, 13/4, 16/5, 19/6, 22/7 and so on? – Somite

@Philip Potter, I clearly labelled it as a challenge. It is intended for people that like code challenges. Like me. Some of us are strange like that! I am looking at methods to estimate a fundamental frequency from potential harmonics. As it is a neat problem, it should have a neat code solution. I will post if I find one. I think it's a good question. It's the sort of question I would answer anyway. It has practical applications, and also some aesthetic appeal. That wins over Sudoku any day. – Warily

@GWW, this criticism should be laid at the feet of SO, for not allowing the asker to accept multiple answers. If you look carefully, you'll see that every single question I've refused to accept an answer on contains no clear winner. Hence I simply upvote the answers I like, as it would be actively misleading to favour one answer over another. The purpose of this feature is to draw attention to one particular answer. Sometimes that works against the clarity of the resource itself. – Warily

Your question is unclear... Note that 13/4 is a better approximation than 3/1. Why isn't that in your list? – Differentiation

Chris, good observation. thank you! I am going to have to re-strategise! – Warily

@Chris (and @Ohmu): The convergents p[k]/q[k] of continued fractions are always best rational approximations, but they don't give all the best approximations. You also have to take the semi-convergents/mediants/whatever: of the form (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) for some integer n≥1 (such that the denominator is less than q[k+1], of course). In particular, you can get 13/4 as a mediant of 1/0 (the conventional "zeroth" convergent) and 3/1: take n=4, so that (1+4*3)/(0+4*1) = 13/4. Another algorithm is here: en.wikipedia.org/wiki/… – Hammad

Great response -- could you possibly put this as an answer? It really belongs there rather than as a comment! – Warily

@Ohmu: Done. I felt like writing a program, so I went ahead and wrote one as well. – Hammad

Answer 1

[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.

Approximating π

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]
...

Program

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;
}

Examples

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!

Hammad

What a fantastic answer! That is far above and beyond anything I had hoped for. It also means that there is some hope for my pitch detection approach. Thank you so much! – Warily

Answer 2

The C program is fine, apart from the fact that you cannot trust the check on the remainder, as can be seen from computing x*p-q as well:

Iteration #1: 3:  3/1 - delta: 0.141592653589793116, rem: 0.141592653589793116
Iteration #2: 7:  22/7 - delta: -0.008851424871448188, rem: 0.062513305931051878
Iteration #3: 15:  333/106 - delta: 0.008821280518070296, rem: 0.996594406684156776
Iteration #4: 1:  355/113 - delta: -0.000030144353377892, rem: 0.003417231014946418
Iteration #5: 292:  103993/33102 - delta: 0.000019129331725765, rem: 0.634590879621879211
Iteration #6: 1:  104348/33215 - delta: -0.000011015021655680, rem: 0.575818424298580694
Iteration #7: 1:  208341/66317 - delta: 0.000008114310077190, rem: 0.736658567704091524
Iteration #8: 1:  312689/99532 - delta: -0.000002900711592702, rem: 0.357480987585133375
Iteration #9: 2:  833719/265381 - delta: 0.000002312886920208, rem: 0.797351564778957706
Iteration #10: 1:  1146408/364913 - delta: -0.000000587824615650, rem: 0.254151925163927682
Iteration #11: 3:  4272943/1360120 - delta: 0.000000549413016415, rem: 0.934654436927838420
Iteration #12: 1:  5419351/1725033 - delta: -0.000000038411599235, rem: 0.069914142051204637
Iteration #13: 14:  80143857/25510582 - delta: 0.000000011648808140, rem: 0.303257833981669641
Iteration #14: 3:  245850922/78256779 - delta: -0.000000003463355824, rem: 0.297524047014214316
Iteration #15: 3:  817696623/260280919 - delta: 0.000000001280568540, rem: 0.361072861287829440
Iteration #16: 2:  1881244168/598818617 - delta: -0.000000000931322575, rem: 0.769524124392304913
Iteration #17: 1:  2698940791/859099536 - delta: 0.000000000232830644, rem: 0.299504418772708979
Iteration #18: 3:  9978066541/3176117225 - delta: 0.000000000000000000, rem: 0.338848902789946401 ******* 'true' deviation below epsilon threshold

Halfmoon