算法挑战：为一个浮点数生成连分数

Question

(编辑: 回应一些不高兴的评论，这不是作业。我正在研究音高检测，需要处理一个潜在的谐波峰数组，并尝试构建基频候选解，所以这实际上是一个非常实际的问题。)

考虑对于（例如）π 的最佳分数逼近，按递增的分母排序：3/1，22/7，355/113，......

挑战：创建一个“整洁”的C算法，将给定浮点数的第n个商式逼近a/b生成，同时返回误差。

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

我认为最好的技术是连分数

去掉pi的小数部分，你得到3
现在，余数是0.14159... = 1/7.06251...

因此下一个最佳的有理数为3 + 1/7 = 22/7
从7.06251中减去7，你得到0.06251...大约是1/15.99659..

将其称为16，那么下一个最佳逼近就是
3 + 1/(7 + 1/16) = 355/113

但是，将其转换为简洁的C代码并不容易。如果我得到了整洁的代码，我会发布一下。与此同时，有人可能会将其视为一个谜题。

- P i

---

如果这确实是一道脑筋急转弯，那么它不属于主题。如果这是你的家庭作业，那么它可能属于主题，但你必须展示给我们你已经尝试了什么以及你遇到了哪些具体问题。我们不会替你完成作业。 - Philip Potter

@Philip Potter，我明确将其标记为挑战。它是为喜欢编程挑战的人而设计的，就像我一样。我们中有些人就是这样奇怪！我正在研究从可能的谐波中估计基本频率的方法。由于这是一个简洁的问题，应该有一个简洁的代码解决方案。如果我找到了一个，我会发布出来。我认为这是一个好问题。这是我无论如何都会回答的问题。它具有实际应用价值，也具有一定的美学吸引力。这比数独更有吸引力。 - P i

@GWW，这个批评应该归咎于SO，因为它不允许提问者接受多个答案。如果你仔细看，你会发现我没能接受答案的每一个问题都没有明显的优胜者。因此，我只是给我喜欢的答案点赞，因为有利偏袒一个答案而不是另一个答案是会导致误导的。这个功能的目的是吸引注意力到一个特定的答案上。有时候这会对资源本身的清晰度产生负面影响。 - P i

你的问题不够清楚明确...请注意，13/4是比3/1更好的近似值，为什么你的列表中没有它？ - Chris Hopman

@Chris（和@Ohmu）：连分数的收敛分数 p[k]/q[k] 总是最佳有理逼近，但它们并不提供全部最佳逼近。您还需要考虑半收敛分数/中间值/等等：形式为 (p[k]+np[k+1])/(q[k]+nq[k+1]) 的分数，其中 n≥1 是某个整数（当然是使得分母小于 q[k+1] 的整数）。特别地，您可以将 13/4 视为 1/0（常规的“零级”收敛分数）和 3/1 的中间值：取 n=4，这样 (1+43)/(0+41) = 13/4。另一个算法在这里：http://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations - ShreevatsaR

Answer 1

[由于您希望得到回答而非评论。]

对于任何实数，其连分数的收敛项 p[k]/q[k] 总是最佳有理逼近，但它们并不是 所有 最佳有理逼近。要获得所有这样的项，还需要取半收敛项/中值——形如 (p[k]+n*p[k+1])/(q[k]+n*q[k+1]) 的分数，其中n≥1为整数。取 n=a[k+2] 得到 p[k+2]/q[k+2]，可取的整数 n 为从 floor(a[k+2]/2) 或 ceiling(a[k+2]/2) 到 a[k+2]。这也在维基百科上被提及。

近似π

π 的连分数为 [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…]（OEIS 序列A001203），其收敛项序列为 3/1, 22/7, 333/106, 355/113, 103993/33102… （A002485/A002486），而最佳逼近项序列为 3/1, 13/4, 16/5, 19/6, 22/7, 179/57… （A063674/A063673）。

因此，算法表明 π = [3; 7, 15, 1, 292, 1, 1,…] 的最佳逼近项为

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

程序

这是一段 C 程序，给定一个正实数，可以生成它的连分数、收敛值以及最佳有理逼近的序列。函数 find_cf 可以找到连分数（将项保存在 a[] 中，将收敛值保存在 p[] 和 q[] 中 - 请忽略全局变量），而函数 all_best 则输出所有最佳有理逼近。

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

示例

对于圆周率π，这个程序的输出如下，在约0.003秒内完成（即，它真正比循环遍历所有可能的分母更好！），为了易读性而换行：

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

所有这些术语都是正确的，但如果您增加MAX，则会因精度问题而出现错误。我对只有13个收敛项就能得到如此多的术语印象深刻。（正如您所看到的，有一个小错误，有时不会打印第一个“n/1”近似值，或者打印不正确 - 我让您来修复！）

您可以尝试√2，其连分数为[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,

或者黄金比例φ = (1+√5)/2，其连分数为[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, 

看到斐波那契数列了吗？这里的收敛数都是近似值。

或者对于有理数，例如4/3 = [1; 3]：

% ./a.out 1.33333333333333333333
1.333333333333333

1:  1/1
3:  4/3

3/2, 4/3, 

或者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, 

祝您愉快！

- ShreevatsaR

---

太棒了！这远远超出了我的期望。这也意味着我的音高检测方法还有一些希望。非常感谢你！ - P i

Answer 2

这个C程序还可以，除了不能信任余数检查，因为计算x*p-q也可以看出这一点：

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

- user2039073