Algorithmus-Challenge: Generiere fortgesetzte Brüche für einen Float

Lesezeit: 7 Minuten

Benutzer-Avatar
P ich

(BEARBEITEN: Als Antwort auf mürrische Kommentare: Nein, es sind keine Hausaufgaben. Ich arbeite an der Tonhöhenerkennung, nehme eine Reihe potenzieller harmonischer Spitzen und versuche, Kandidaten für die Grundfrequenz zu konstruieren. Es ist also eigentlich eine sehr praktische Frage.)

Betrachten Sie die besten gebrochenen Näherungen für (z. B.) pi, geordnet nach aufsteigendem Nenner: 3/1, 22/7, 355/113, …

Herausforderung: Erstellen Sie eine ordentlich C-Algorithmus, der die n-te Quotientennäherung a/b für einen gegebenen Float generiert und auch die Diskrepanz zurückgibt.

calcBestFrac(float frac, int n, int * a, int * b, float * err) {…}

Die beste Technik, die ich glaube, ist fortgesetzte Brüche

Nimm den Bruchteil von Pi weg und du bekommst 3
Jetzt ist der Rest 0,14159 … = 1/7,06251 …

Die nächstbeste Ratio ist also 3 + 1/7 = 22/7
Nehmen Sie die 7 von 7,06251 weg und Sie erhalten 0,06251.. Ungefähr 1/15,99659..

Nennen Sie es 16, dann ist die nächstbeste Annäherung
3 + 1/(7 + 1/16) = 355/113

Dies ist jedoch alles andere als trivial in sauberen C-Code umzuwandeln. Ich werde posten, wenn ich etwas ordentlich bekomme. In der Zwischenzeit mag es jemand als Denksportaufgabe genießen.

  • Wenn dies wirklich eine Denksportaufgabe ist, ist es kein Thema. Wenn es tatsächlich Ihre Hausaufgabe ist, könnte es themenbezogen sein, aber Sie müssen uns zeigen, was Sie bisher versucht haben und welches spezifische Problem Sie haben. Wir erledigen nicht nur Ihre Arbeit für Sie.

    – Philip Potter

    9. Januar 2011 um 6:28 Uhr

  • @Philip Potter, ich habe es eindeutig als Herausforderung bezeichnet. Es ist für Leute gedacht, die Code-Herausforderungen mögen. Wie ich. Manche von uns sind so seltsam! Ich suche nach Methoden, um eine Grundfrequenz aus potenziellen Oberschwingungen abzuschätzen. Da es sich um ein ordentliches Problem handelt, sollte es eine ordentliche Codelösung haben. Ich werde posten, wenn ich einen finde. Ich denke, es ist eine gute Frage. Das ist die Art von Frage, die ich sowieso beantworten würde. Es hat praktische Anwendungen und auch einen gewissen ästhetischen Reiz. Das gewinnt jeden Tag über Sudoku.

    – P ich

    9. Januar 2011 um 8:01 Uhr

  • @GWW, diese Kritik sollte SO zu Füßen gelegt werden, weil es dem Fragesteller nicht erlaubt ist, mehrere Antworten zu akzeptieren. Wenn Sie genau hinsehen, werden Sie feststellen, dass jede einzelne Frage, auf die ich mich geweigert habe, eine Antwort zu akzeptieren, keinen klaren Gewinner enthält. Daher stimme ich einfach den Antworten zu, die mir gefallen, da es aktiv irreführend wäre, eine Antwort einer anderen vorzuziehen. Der Zweck dieser Funktion besteht darin, die Aufmerksamkeit auf eine bestimmte Antwort zu lenken. Manchmal wirkt sich das gegen die Klarheit der Ressource selbst aus.

    – P ich

    9. Januar 2011 um 8:05 Uhr

  • Ihre Frage ist unklar … Beachten Sie, dass 13/4 eine bessere Annäherung ist als 3/1. Warum ist das nicht in deiner Liste?

    – Chris Hopmann

    9. Januar 2011 um 9:50 Uhr

  • @Chris (und @Ohmu): Die Konvergenten p[k]/q[k] von Kettenbrüchen sind immer die besten rationalen Annäherungen, aber sie geben nicht alle die besten Annäherungen. Sie müssen auch die Halbkonvergenten/Medianten/was auch immer: der Form (S[k]+n*p[k+1])/(q[k]+n*q[k+1]) für eine ganze Zahl n≥1 (so dass der Nenner kleiner als q ist[k+1], Natürlich). Insbesondere können Sie 13/4 als Mediant von 1/0 (dem herkömmlichen „nullten“ Konvergenten) und 3/1 erhalten: Nehmen Sie n=4, sodass (1+4*3)/(0+4*1 ) = 13/4. Ein weiterer Algorithmus ist hier: de.wikipedia.org/wiki/…

    – ShreevatsaR

    9. Januar 2011 um 11:36 Uhr

Benutzer-Avatar
ShreevatsaR

[Since you asked for this as an answer rather than a comment.]

Für jede reelle Zahl sind die Konvergenten p[k]/q[k] seines fortgesetzten Bruchs sind immer die besten rationalen Annäherungen, aber das sind sie nicht alle die besten rationalen Annäherungen. Um alle zu erhalten, müssen Sie auch die Halbkonvergenten/Medianten nehmen – Bruchteile der Form (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) für eine ganze Zahl n≥1. Nimm n=a[k+2] gibt p[k+2]/q[k+2]und die zu nehmenden ganzen Zahlen n sind die von beiden Stockwerken (a[k+2]/2) oder Decke(a[k+2]/2), zu a[k+2]. Dies wird auch erwähnt auf Wikipedia.

Annäherung an π

Der fortgesetzte Bruch für π ist [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…] (Reihenfolge A001203 im OEIS), ist die Konvergenzfolge 3/1, 22/7, 333/106, 355/113, 103993/33102… (A002485/A002486), und die Reihenfolge der besten Annäherungen ist 3/1, 13/4, 16/5, 19/6, 22/7, 179/57… (A063674/A063673).

Der Algorithmus sagt also, dass die besten Näherungen von π = [3; 7, 15, 1, 292, 1, 1,…] sind

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

Programm

Hier ist ein C-Programm, das bei einer gegebenen positiven reellen Zahl ihren fortgesetzten Bruch, ihre Konvergenten und die Folge der besten rationalen Näherungen erzeugt. Die Funktion find_cf findet den fortgesetzten Bruch (setzt die Terme in a[] und die Konvergenzen in p[] und q[] — entschuldigen Sie die globalen Variablen) und die Funktion all_best druckt die besten rationalen Annäherungen.

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

Beispiele

Für π ist hier die Ausgabe dieses Programms in etwa 0,003 Sekunden (dh es ist wirklich besser, als alle möglichen Nenner zu durchlaufen!), Zeilenumbruch für die Lesbarkeit:

% ./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, 

Alle diese Begriffe sind korrekt, aber wenn Sie MAX erhöhen, erhalten Sie aufgrund der Genauigkeit Fehler. Ich bin selbst beeindruckt, wie viele Terme man mit nur 13 Konvergenten bekommt. (Wie Sie sehen können, gibt es einen kleinen Fehler, bei dem manchmal die allererste „n/1“-Annäherung nicht oder falsch gedruckt wird – die Behebung überlasse ich Ihnen!)

Du kannst es mit √2 versuchen, dessen Kettenbruch ist [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,

Oder der Goldene Schnitt φ = (1+√5)/2 dessen fortgesetzter Bruch ist [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, 

(Siehe die Fibonacci-Zahlen? Hier sind die Konvergenten alle Näherungen.)

Oder mit rationalen Zahlen wie 4/3 = [1; 3]:

% ./a.out 1.33333333333333333333
1.333333333333333

1:  1/1
3:  4/3

3/2, 4/3, 

oder 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, 

Genießen!

  • Was für eine fantastische Antwort! Das geht weit über alles hinaus, was ich mir erhofft hatte. Es bedeutet auch, dass es etwas Hoffnung für meinen Tonhöhenerkennungsansatz gibt. Ich danke dir sehr!

    – P ich

    10. Januar 2011 um 0:32 Uhr

Das C-Programm ist in Ordnung, abgesehen davon, dass Sie der Überprüfung des Rests nicht trauen können, wie auch aus der Berechnung von x*pq ersichtlich ist:

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

1057570cookie-checkAlgorithmus-Challenge: Generiere fortgesetzte Brüche für einen Float

This website is using cookies to improve the user-friendliness. You agree by using the website further.

Privacy policy