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

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

[因为您要求的是答案而不是评论。]

对于任何实数，其连续分数的收敛点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是下限（a [k + 2] / 2）或上限（a [ k + 2] / 2）到a [k +
2]。维基百科上也提到了这一点。

约π

π的连续分数为[3; 7，15，1，292，1，1，1，2，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-9long 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/%ldn", 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("%.15lfnn", x);  all_best(x); printf("n");  return 0;}

例子

对于π，这是该程序的输出，大约需要0.003秒（即，这比循环遍历所有可能的分母要好！），并进行换行以提高可读性：

% ./a.out3.1415926535897933:  3/17:  22/715:  333/1061:  355/113292:  103993/331021:  104348/332151:  208341/663171:  312689/995322:  833719/2653811:  1146408/3649133:  4272943/13601201:  5419351/172503314:  80143857/2551058213/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.414213562373095048801.4142135623730951:  1/12:  3/22:  7/52:  17/122:  41/292:  99/702:  239/1692:  577/4082:  1393/9852:  3363/23782:  8119/57412:  19601/138602:  47321/334613/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.618033988749894848201.6180339887498951:  1/11:  2/11:  3/21:  5/31:  8/51:  13/81:  21/131:  34/211:  55/341:  89/551:  144/891:  233/1441:  377/2332/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.333333333333333333331.3333333333333331:  1/13:  4/33/2, 4/3,

或14/11 = [1; 3，1，2]：

% ./a.out 1.272727272727272727271.2727272727272731:  1/13:  4/31:  5/42:  14/113/2, 4/3, 5/4, 9/7, 14/11,

请享用！